Ignore:
Timestamp:
Jul 28, 2025, 6:44:28 PM (4 months ago)
Author:
aborella
Message:

Major modifs to treatment of contrails (from 2 classes to 2 moments) + diagnostics. Increased numerical efficiency

Location:
LMDZ6/branches/contrails
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/contrails/DefLists/context_input_lmdz.xml

    r5779 r5790  
    409409  <file id="aviation_file" name="aviation" enabled="false" mode="read" output_freq="1mo" type="one_file" time_counter_name="toto" >
    410410
    411     <field id="KMFLOWN_id" name="seg_length_m" operation="instant"  grid_ref="aviation_grid" freq_offset="1ts"  />
     411    <field id="DISTFLOWN_id" name="seg_length" operation="instant"  grid_ref="aviation_grid" freq_offset="1ts"  />
     412    <field id="FUELBURN_id" name="fuel_burn" operation="instant"  grid_ref="aviation_grid" freq_offset="1ts"  />
    412413    <field id="levaviation_id" name="pressure_Pa" axis_ref="aviation_lev" operation="instant" freq_offset="1ts" />
    413414    <field id="timeaviation_id" name="time" axis_ref="aviation_time" operation="instant" freq_offset="1ts" />
    … …  
    444445<field_definition>
    445446
    446   <field id ="KMFLOWN_read" field_ref="KMFLOWN_id"  enabled="false" read_access="true" />
    447   <field id ="KMFLOWN_interp" field_ref="KMFLOWN_read"  enabled="false" read_access="true" grid_ref="grid_from_aviation" />
     447  <field id ="DISTFLOWN_read" field_ref="DISTFLOWN_id"  enabled="false" read_access="true" />
     448  <field id ="DISTFLOWN_interp" field_ref="DISTFLOWN_read"  enabled="false" read_access="true" grid_ref="grid_from_aviation" />
     449  <field id ="FUELBURN_read" field_ref="FUELBURN_id"  enabled="false" read_access="true" />
     450  <field id ="FUELBURN_interp" field_ref="FUELBURN_read"  enabled="false" read_access="true" grid_ref="grid_from_aviation" />
    448451
    449452</field_definition>

  • LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml

    r5779 r5790  
    920920        <field id="qsati"      long_name="Saturation with respect to ice"    unit="kg/kg" />
    921921
    922         <field id="cflseri"    long_name="Linear contrail fraction" unit="-" />
    923         <field id="dcfldyn"    long_name="Dynamics linear contrail fraction tendency" unit="s-1" />
    924         <field id="cfcseri"    long_name="Contrail induced cirrus fraction" unit="-" />
    925         <field id="dcfcdyn"    long_name="Dynamics contrail induced cirrus fraction tendency" unit="s-1" />
    926         <field id="qtlseri"    long_name="Linear contrail total specific humidity" unit="kg/kg" />
    927         <field id="dqtldyn"    long_name="Dynamics linear contrail total specific humidity tendency" unit="kg/kg/s" />
    928         <field id="qtcseri"    long_name="Contrail cirrus total specific humidity" unit="kg/kg" />
    929         <field id="dqtcdyn"    long_name="Dynamics contrail cirrus total specific humidity tendency" unit="kg/kg/s" />
     922        <field id="cfcseri"    long_name="Contrail fraction" unit="-" />
     923        <field id="dcfcdyn"    long_name="Dynamics contrail fraction tendency" unit="s-1" />
     924        <field id="qtcseri"    long_name="Contrail total specific humidity" unit="kg/kg" />
     925        <field id="dqtcdyn"    long_name="Dynamics contrail total specific humidity tendency" unit="kg/kg/s" />
     926        <field id="nicseri"    long_name="Contrail ice crystals concentration" unit="#/kg" />
     927        <field id="dnicdyn"    long_name="Dynamics contrail ice crystals concentration tendency" unit="#/kg/s" />
    930928        <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
    931929        <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
    932930        <field id="potcontfraP"  long_name="Fraction with pontential persistent contrail"    unit="-" />
    933931        <field id="potcontfraNP" long_name="Fraction with potential non-persistent contrail"    unit="-" />
    934         <field id="qice_lincont" long_name="Linear contrails ice specific humidity"    unit="kg/kg" />
    935         <field id="qice_circont" long_name="Contrail induced cirrus ice specific humidity"    unit="kg/kg" />
    936         <field id="dcflini"    long_name="Initial formation linear contrail fraction tendency"    unit="s-1" />
    937         <field id="dqilini"    long_name="Initial formation linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
    938         <field id="dqtlini"    long_name="Initial formation linear contrail total specific humidity tendency"    unit="kg/kg/s" />
    939         <field id="dcflcir"    long_name="Conversion to cirrus linear contrail fraction tendency"    unit="s-1" />
    940         <field id="dqtlcir"    long_name="Conversion to cirrus linear contrail total specific humidity tendency"    unit="kg/kg/s" />
    941         <field id="dcflsub"    long_name="Sublimation linear contrail fraction tendency"    unit="s-1" />
    942         <field id="dqilsub"    long_name="Sublimation linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
    943         <field id="dqtlsub"    long_name="Sublimation linear contrail total specific humidity tendency"    unit="kg/kg/s" />
    944         <field id="dcflmix"    long_name="Mixing linear contrail fraction tendency"    unit="s-1" />
    945         <field id="dqilmix"    long_name="Mixing linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
    946         <field id="dqtlmix"    long_name="Mixing linear contrail total specific humidity tendency"    unit="kg/kg/s" />
    947         <field id="dcflsed"    long_name="Ice sedimentation linear contrail fraction tendency"    unit="s-1" />
    948         <field id="dqilsed"    long_name="Ice sedimentation linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
    949         <field id="dqtlsed"    long_name="Ice sedimentation linear contrail total specific humidity tendency"    unit="kg/kg/s" />
    950         <field id="dcflauto"   long_name="Ice autoconversion linear contrail fraction tendency"    unit="s-1" />
    951         <field id="dqilauto"   long_name="Ice autoconversion linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
    952         <field id="dqtlauto"   long_name="Ice autoconversion linear contrail total specific humidity tendency"    unit="kg/kg/s" />
    953         <field id="dcfcsub"    long_name="Sublimation contrail cirrus fraction tendency"    unit="s-1" />
    954         <field id="dqicsub"    long_name="Sublimation contrail cirrus ice specific humidity tendency"    unit="kg/kg/s" />
    955         <field id="dqtcsub"    long_name="Sublimation contrail cirrus total specific humidity tendency"    unit="kg/kg/s" />
    956         <field id="dcfcmix"    long_name="Mixing contrail cirrus fraction tendency"    unit="s-1" />
    957         <field id="dqicmix"    long_name="Mixing contrail cirrus ice specific humidity tendency"    unit="kg/kg/s" />
    958         <field id="dqtcmix"    long_name="Mixing contrail cirrus total specific humidity tendency"    unit="kg/kg/s" />
    959         <field id="dcfcsed"    long_name="Ice sedimentation contrail cirrus fraction tendency"    unit="s-1" />
    960         <field id="dqicsed"    long_name="Ice sedimentation contrail cirrus ice specific humidity tendency"    unit="kg/kg/s" />
    961         <field id="dqtcsed"    long_name="Ice sedimentation contrail cirrus total specific humidity tendency"    unit="kg/kg/s" />
    962         <field id="dcfcauto"   long_name="Ice autoconversion contrail cirrus fraction tendency"    unit="s-1" />
    963         <field id="dqicauto"   long_name="Ice autoconversion contrail cirrus ice specific humidity tendency"    unit="kg/kg/s" />
    964         <field id="dqtcauto"   long_name="Ice autoconversion contrail cirrus total specific humidity tendency"    unit="kg/kg/s" />
     932        <field id="qice_cont"  long_name="Contrails ice specific humidity"    unit="kg/kg" />
     933        <field id="dcfcini"    long_name="Initial formation contrail fraction tendency"    unit="s-1" />
     934        <field id="dqicini"    long_name="Initial formation contrail ice specific humidity tendency"    unit="kg/kg/s" />
     935        <field id="dqtcini"    long_name="Initial formation contrail total specific humidity tendency"    unit="kg/kg/s" />
     936        <field id="dnicini"    long_name="Initial formation contrail ice crystals concentration tendency"    unit="#/kg/s" />
     937        <field id="dcfcsub"    long_name="Sublimation contrail fraction tendency"    unit="s-1" />
     938        <field id="dqicsub"    long_name="Sublimation contrail ice specific humidity tendency"    unit="kg/kg/s" />
     939        <field id="dqtcsub"    long_name="Sublimation contrail total specific humidity tendency"    unit="kg/kg/s" />
     940        <field id="dnicsub"    long_name="Sublimation contrail ice crystals concentration tendency"    unit="#/kg/s" />
     941        <field id="dcfcmix"    long_name="Mixing contrail fraction tendency"    unit="s-1" />
     942        <field id="dqicmix"    long_name="Mixing contrail ice specific humidity tendency"    unit="kg/kg/s" />
     943        <field id="dqtcmix"    long_name="Mixing contrail total specific humidity tendency"    unit="kg/kg/s" />
     944        <field id="dnicmix"    long_name="Mixing contrail ice crystals concentration tendency"    unit="#/kg/s" />
     945        <field id="dnicagg"    long_name="Aggregation contrail ice crystals concentration tendency"    unit="#/kg/s" />
     946        <field id="dcfcsed"    long_name="Ice sedimentation contrail fraction tendency"    unit="s-1" />
     947        <field id="dqicsed"    long_name="Ice sedimentation contrail ice specific humidity tendency"    unit="kg/kg/s" />
     948        <field id="dqtcsed"    long_name="Ice sedimentation contrail total specific humidity tendency"    unit="kg/kg/s" />
     949        <field id="dnicsed"    long_name="Ice sedimentation contrail ice crystals concentration tendency"    unit="#/kg/s" />
     950        <field id="dcfcauto"   long_name="Ice autoconversion contrail fraction tendency"    unit="s-1" />
     951        <field id="dqicauto"   long_name="Ice autoconversion contrail ice specific humidity tendency"    unit="kg/kg/s" />
     952        <field id="dqtcauto"   long_name="Ice autoconversion contrail total specific humidity tendency"    unit="kg/kg/s" />
     953        <field id="dnicauto"   long_name="Ice autoconversion contrail ice crystals concentration tendency"    unit="#/kg/s" />
    965954        <field id="flightdist" long_name="Aviation flown distance concentration"    unit="m/s/m3" />
    966         <field id="flighth2o"  long_name="Aviation emitted H2O concentration"    unit="kg H2O/s/m3" />
    967         <field id="dqavi"      long_name="Water vapor emissions from aviation tendency"    unit="kg/kg/s" />
     955        <field id="flightfuel" long_name="Aviation fuel consumption concentration"    unit="kg/s/m3" />
     956        <field id="AEI_cont"   long_name="Apparent emission index contrails"    unit="#/kg" />
     957        <field id="AEI_surv_cont" long_name="Apparent emission index contrails after vortex loss"    unit="#/kg" />
     958        <field id="fsurv_cont"    long_name="Survival fraction after vortex loss"    unit="-" />
     959        <field id="section_cont"  long_name="Cross section of newly formed contrails"    unit="m2" />
     960        <field id="dqavi"         long_name="Water vapor emissions from aviation tendency"    unit="kg/kg/s" />
    968961        <field id="cldfra_nocont" long_name="Cloud fraction w/o contrails"    unit="-" />
    969962        <field id="cldtau_nocont" long_name="Cloud optical thickness w/o contrails"    unit="1" />

  • LMDZ6/branches/contrails/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90

    r5641 r5790  
    5050    ale_bl, ale_bl_trig, alp_bl, &
    5151    ale_wake, ale_bl_stat, AWAKE_S, &
    52     cf_ancien, qvc_ancien, qvcon, qccon, cfl_ancien, cfc_ancien, qtl_ancien, qtc_ancien
     52    cf_ancien, qvc_ancien, qvcon, qccon, cfc_ancien, qtc_ancien, nic_ancien
    5353
    5454  USE comconst_mod, ONLY: pi, dtvr
    … …  
    245245  qvcon = 0.
    246246  qccon = 0.
    247   cfl_ancien = 0.
    248247  cfc_ancien = 0.
    249   qtl_ancien = 0.
    250248  qtc_ancien = 0.
     249  nic_ancien = 0.
    251250 
    252251  z0m(:,:)=0 ! ym missing 5th subsurface initialization

  • LMDZ6/branches/contrails/libf/phylmd/infotrac_phy.F90

    r5684 r5790  
    2727   !=== FOR ISOTOPES: Specific to water
    2828   PUBLIC :: iH2O                                          !--- Value of "ixIso" for "H2O" isotopes class
    29    PUBLIC :: ivap, iliq, isol, ibs, icf, iqvc, icfl, icfc, iqtl, iqtc
     29   PUBLIC :: ivap, iliq, isol, ibs, icf, iqvc, icfc, iqtc, inic
    3030   !=== FOR ISOTOPES: Depending on the selected isotopes family
    3131   PUBLIC :: isotope                                       !--- Selected isotopes database (argument of getKey)
    … …  
    103103
    104104   !=== INDICES FOR WATER
    105    INTEGER, SAVE :: ivap, iliq, isol, ibs, icf, iqvc, icfl, icfc, iqtl, iqtc
    106 !$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, icfl, icfc, iqtl, iqtc)
     105   INTEGER, SAVE :: ivap, iliq, isol, ibs, icf, iqvc, icfc, iqtc, inic
     106!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, icfc, iqtc, inic)
    107107
    108108   !=== DIMENSIONS OF THE TRACERS TABLES AND OTHER SCALAR VARIABLES
    … …  
    365365   icf  = strIdx(tracers(:)%name, 'CLDFRA')
    366366   iqvc = strIdx(tracers(:)%name, 'CLDVAP')
    367    icfl = strIdx(tracers(:)%name, 'LINCONTFRA')
    368    icfc = strIdx(tracers(:)%name, 'CIRCONTFRA')
    369    iqtl = strIdx(tracers(:)%name, 'LINCONTWATER')
    370    iqtc = strIdx(tracers(:)%name, 'CIRCONTWATER')
     367   icfc = strIdx(tracers(:)%name, 'CONTFRA')
     368   iqtc = strIdx(tracers(:)%name, 'CONTWATER')
     369   inic = strIdx(tracers(:)%name, 'CONTNICE')
    371370   !--Two ways of declaring tracers - the way below should be deleted in the future
    372371   IF (icf.EQ.0)  icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
    373372   IF (iqvc.EQ.0) iqvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
    374    IF (icfl.EQ.0) icfl = strIdx(tracers(:)%name, addPhase('H2O', 'z'))
    375    IF (icfc.EQ.0) icfc = strIdx(tracers(:)%name, addPhase('H2O', 'y'))
    376    IF (iqtl.EQ.0) iqtl = strIdx(tracers(:)%name, addPhase('H2O', 'x'))
    377    IF (iqtc.EQ.0) iqtc = strIdx(tracers(:)%name, addPhase('H2O', 'w'))
     373   IF (icfc.EQ.0) icfc = strIdx(tracers(:)%name, addPhase('H2O', 'x'))
     374   IF (iqtc.EQ.0) iqtc = strIdx(tracers(:)%name, addPhase('H2O', 'y'))
     375   IF (inic.EQ.0) inic = strIdx(tracers(:)%name, addPhase('H2O', 'z'))
    378376
    379377   IF(CPPKEY_STRATAER .AND. type_trac == 'coag') THEN

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

    r5779 r5790  
    99! The size is (klon, nleva, 1) where
    1010! nleva            is the size of the vertical axis (read from file)
    11 ! flight_dist_read is the number of km per second
    12 ! flight_h2o_read  is the water content added to the air
     11! flight_dist_read is the number of km per second per m2
     12! flight_fuel_read is the fuel consumed per second per m2
    1313! aviation_lev     is the value of the levels
    1414INTEGER, SAVE :: nleva  ! Size of the vertical axis in the file
    1515!$OMP THREADPRIVATE(nleva)
    16 REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_dist_read ! Aviation distance flown within the mesh [m/s/mesh]
     16REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_dist_read ! Aviation distance flown within the mesh [m/s/m2/vertmesh]
    1717!$OMP THREADPRIVATE(flight_dist_read)
    18 REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_h2o_read  ! Aviation H2O emitted within the mesh [kgH2O/s/mesh]
    19 !$OMP THREADPRIVATE(flight_h2o_read)
     18REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE, PRIVATE :: flight_fuel_read ! Aviation fuel consumed within the mesh [kg/s/m2/vertmesh]
     19!$OMP THREADPRIVATE(flight_fuel_read)
    2020REAL, ALLOCATABLE, DIMENSION(:),     SAVE, PRIVATE :: aviation_lev     ! Pressure in the middle of the layers [Pa]
    2121!$OMP THREADPRIVATE(aviation_lev)
    … …  
    2525SUBROUTINE aviation_water_emissions( &
    2626      klon, klev, dtime, pplay, temp, qtot, &
    27       flight_h2o, d_q_avi &
     27      flight_fuel, d_q_avi &
    2828      )
    2929
    30 USE lmdz_lscp_ini, ONLY: RD
     30USE lmdz_lscp_ini, ONLY: RD, EI_H2O_aviation
    3131
    3232IMPLICIT NONE
    3333
    34 INTEGER,                      INTENT(IN)  :: klon, klev ! number of horizontal grid points and vertical levels
    35 REAL,                         INTENT(IN)  :: dtime      ! time step [s]
    36 REAL, DIMENSION(klon,klev),   INTENT(IN)  :: pplay      ! mid-layer pressure [Pa]
    37 REAL, DIMENSION(klon,klev),   INTENT(IN)  :: temp       ! temperature (K)
    38 REAL, DIMENSION(klon,klev),   INTENT(IN)  :: qtot       ! total specific humidity (in vapor phase) [kg/kg]
    39 REAL, DIMENSION(klon,klev),   INTENT(IN)  :: flight_h2o ! aviation emitted H2O concentration [kgH2O/s/m3]
    40 REAL, DIMENSION(klon,klev),   INTENT(OUT) :: d_q_avi    ! water vapor tendency from aviation [kg/kg]
     34INTEGER,                      INTENT(IN)  :: klon, klev  ! number of horizontal grid points and vertical levels
     35REAL,                         INTENT(IN)  :: dtime       ! time step [s]
     36REAL, DIMENSION(klon,klev),   INTENT(IN)  :: pplay       ! mid-layer pressure [Pa]
     37REAL, DIMENSION(klon,klev),   INTENT(IN)  :: temp        ! temperature (K)
     38REAL, DIMENSION(klon,klev),   INTENT(IN)  :: qtot        ! total specific humidity (in vapor phase) [kg/kg]
     39REAL, DIMENSION(klon,klev),   INTENT(IN)  :: flight_fuel ! aviation fuel consumption concentration [kg/s/m3]
     40REAL, DIMENSION(klon,klev),   INTENT(OUT) :: d_q_avi     ! water vapor tendency from aviation [kg/kg]
    4141! Local
    4242INTEGER :: i, k
    43 REAL :: rho
     43REAL :: rho, flight_h2o
    4444
    4545DO i=1, klon
    … …  
    4747    !--Dry density [kg/m3]
    4848    rho = pplay(i,k) / temp(i,k) / RD
     49    flight_h2o = flight_fuel(i,k) * EI_H2O_aviation
    4950
    5051    !--q is the specific humidity (kg/kg humid air) hence the complicated equation to update q
    … …  
    5455    !--flight_h2O is in kg H2O / s / m3
    5556   
    56     !d_q_avi(i,k) = ( M_cell * qtot(i,k) + V_cell * flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
    57     !             / ( M_cell             + V_cell * flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
     57    !d_q_avi(i,k) = ( M_cell * qtot(i,k) + V_cell * flight_h2o * dtime * ( 1. - qtot(i,k) ) ) &
     58    !             / ( M_cell             + V_cell * flight_h2o * dtime * ( 1. - qtot(i,k) ) ) &
    5859    !             - qtot(i,k)
    5960    !--NB., M_cell = V_cell * rho
    60     !d_q_avi(i,k) = ( rho * qtot(i,k) + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
    61     !             / ( rho             + flight_h2o(i,k) * dtime * ( 1. - qtot(i,k) ) ) &
     61    !d_q_avi(i,k) = ( rho * qtot(i,k) + flight_h2o * dtime * ( 1. - qtot(i,k) ) ) &
     62    !             / ( rho             + flight_h2o * dtime * ( 1. - qtot(i,k) ) ) &
    6263    !             - qtot(i,k)
    6364    !--Same formula, more computationally effective but less readable
    64     d_q_avi(i,k) = flight_h2o(i,k) * ( 1. - qtot(i,k) ) &
    65                  / ( rho / dtime / ( 1. - qtot(i,k) ) + flight_h2o(i,k) )
     65    d_q_avi(i,k) = flight_h2o * ( 1. - qtot(i,k) ) &
     66                 / ( rho / dtime / ( 1. - qtot(i,k) ) + flight_h2o )
    6667  ENDDO
    6768ENDDO
    … …  
    7273!**********************************************************************************
    7374SUBROUTINE contrails_formation( &
    74       klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, flight_dist, &
     75      klon, dtime, pplay, temp, rho, qsat, qsatl, gamma_cond, flight_dist, flight_fuel, &
    7576      clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, keepgoing, pt_pron_clds, &
    7677      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    77       dcfl_ini, dqil_ini, dqtl_ini)
     78      AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, &
     79      dcfc_ini, dqic_ini, dqtc_ini, dnic_ini)
    7880
    7981USE lmdz_lscp_ini, ONLY: RPI, RCPD, EPS_W, RTT
    8082USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation
    8183USE lmdz_lscp_ini, ONLY: eps, temp_nowater
     84USE lmdz_lscp_ini, ONLY: Nice_init_contrails
    8285
    8386USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf
    … …  
    9295REAL,     INTENT(IN) , DIMENSION(klon) :: pplay        ! layer pressure [Pa]
    9396REAL,     INTENT(IN) , DIMENSION(klon) :: temp         ! temperature [K]
     97REAL,     INTENT(IN) , DIMENSION(klon) :: rho          ! air density [kg/m3]
    9498REAL,     INTENT(IN) , DIMENSION(klon) :: qsat         ! saturation specific humidity [kg/kg]
    9599REAL,     INTENT(IN) , DIMENSION(klon) :: qsatl        ! saturation specific humidity w.r.t. liquid [kg/kg]
    96100REAL,     INTENT(IN) , DIMENSION(klon) :: gamma_cond   ! condensation threshold w.r.t. qsat [-]
    97101REAL,     INTENT(IN) , DIMENSION(klon) :: flight_dist  ! aviation distance flown concentration [m/s/m3]
     102REAL,     INTENT(IN) , DIMENSION(klon) :: flight_fuel  ! aviation fuel consumption concentration [kg/s/m3]
    98103REAL,     INTENT(IN) , DIMENSION(klon) :: clrfra       ! clear sky fraction [-]
    99104REAL,     INTENT(IN) , DIMENSION(klon) :: qclr         ! clear sky specific humidity [kg/kg]
    … …  
    110115REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
    111116REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-]
    112 REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfl_ini     ! linear contrails cloud fraction tendency bec ause of initial formation [s-1]
    113 REAL,     INTENT(INOUT), DIMENSION(klon) :: dqil_ini     ! linear contrails ice specific humidity tende ncy because of initial formation [kg/kg/s]
    114 REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtl_ini     ! linear contrails total specific humidity ten dency because of initial formation [kg/kg/s]
     117REAL,     INTENT(INOUT), DIMENSION(klon) :: AEI_contrails ! Apparent emission index contrails [#/kg]
     118REAL,     INTENT(INOUT), DIMENSION(klon) :: AEI_surv_contrails ! Apparent emission index contrails after vortex loss [#/kg]
     119REAL,     INTENT(INOUT), DIMENSION(klon) :: fsurv_contrails ! Survival fraction after vortex loss [-]
     120REAL,     INTENT(INOUT), DIMENSION(klon) :: section_contrails ! Cross section of newly formed contrails [m2]
     121REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_ini     ! contrails cloud fraction tendency bec ause of initial formation [s-1]
     122REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_ini     ! contrails ice specific humidity tende ncy because of initial formation [kg/kg/s]
     123REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_ini     ! contrails total specific humidity ten dency because of initial formation [kg/kg/s]
     124REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_ini     ! contrails ice crystals concentration ten dency because of initial formation [#/kg/s]
    115125!
    116126! Local
    … …  
    128138!
    129139! for new contrail formation
    130 REAL :: contrail_cross_section, contfra_new
     140REAL :: dist_contrails, Nice_per_m_init_contrails
     141REAL :: icesat_ratio
    131142
    132143qzero(:) = 0.
    … …  
    213224   
    214225    !--Add a source of contrails from aviation
    215     IF ( potcontfraP(i) .GT. eps ) THEN
    216       contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
    217       contfra_new = MIN(1., flight_dist(i) * dtime * contrail_cross_section)
    218       dcfl_ini(i) = potcontfraP(i) * contfra_new
    219       dqtl_ini(i) = qpotcontP * contfra_new
    220       dqil_ini(i) = dqtl_ini(i) - qsat(i) * dcfl_ini(i)
     226    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( &
     229              dtime, rho(i), flight_dist(i), flight_fuel(i))
     230      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
     234      CALL CONTRAIL_ICE_NUMBER_CONCENTRATION(temp(i), icesat_ratio, rho(i), &
     235              flight_dist(i), flight_fuel(i), Nice_per_m_init_contrails, &
     236              AEI_contrails(i), AEI_surv_contrails(i), fsurv_contrails(i))
     237
     238      !--If Nice_per_m_init_contrails .EQ. 0., all the crystals were lost in the vortex
     239      IF ( Nice_per_m_init_contrails .GT. 0. ) THEN
     240        !--dist_contrails is in meters of contrails per m3 (gridbox-average)
     241        dist_contrails = flight_dist(i) * dtime * potcontfraP(i)
     242        dcfc_ini(i) = dist_contrails * section_contrails(i)
     243        dqtc_ini(i) = icesat_ratio * qsat(i) * dcfc_ini(i)
     244        dqic_ini(i) = dqtc_ini(i) - qsat(i) * dcfc_ini(i)
     245        dnic_ini(i) = Nice_per_m_init_contrails * dist_contrails / rho(i)
     246      ENDIF
    221247    ENDIF
    222248
    … …  
    228254
    229255!**********************************************************************************
    230 FUNCTION contrail_cross_section_onera()
     256FUNCTION CONTRAIL_CROSS_SECTION_ONERA()
    231257
    232258USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
    … …  
    244270contrail_cross_section_onera = initial_width_contrails * initial_height_contrails
    245271
    246 END FUNCTION contrail_cross_section_onera
     272END FUNCTION CONTRAIL_CROSS_SECTION_ONERA
     273
     274
     275!**********************************************************************************
     276!--Based on Schumann (1998)
     277FUNCTION CONTRAIL_CROSS_SECTION_SCHUMANN(dtime, rho_air, flight_dist, flight_fuel)
     278
     279USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
     280
     281IMPLICIT NONE
     282
     283!
     284! Input
     285!
     286REAL :: dtime       ! timestep [s]
     287REAL :: rho_air     ! air density [kg/m3]
     288REAL :: flight_dist ! flown distance [m/s/m3]
     289REAL :: flight_fuel ! fuel consumed [kg/s/m3]
     290!
     291! Output
     292!
     293REAL :: contrail_cross_section_schumann ! [m2]
     294!
     295! Local
     296!
     297REAL :: dilution_factor
     298
     299!--The contrail is formed on average in the middle of the timestep
     300dilution_factor = 7000. * (dtime / 2.)**0.8
     301contrail_cross_section_schumann = flight_fuel / flight_dist * dilution_factor / rho_air
     302
     303END FUNCTION CONTRAIL_CROSS_SECTION_SCHUMANN
     304
     305
     306!**********************************************************************************
     307SUBROUTINE CONTRAIL_ICE_NUMBER_CONCENTRATION(temp, icesat_ratio, rho_air, &
     308        flight_dist, flight_fuel, Nice_per_m_init_contrails, &
     309        AEI_contrails, AEI_surv_contrails, fsurv_contrails)
     310
     311USE lmdz_lscp_ini, ONLY: RPI
     312USE lmdz_lscp_ini, ONLY: EI_soot_aviation, air_to_fuel_ratio_engine, wingspan
     313USE lmdz_lscp_ini, ONLY: Naer_amb, raer_amb_mean, raer_amb_std, r_soot_mean, r_soot_std
     314USE lmdz_lscp_ini, ONLY: circ_0_loss, plume_area_loss, nice_init_ref_loss
     315USE lmdz_lscp_ini, ONLY: N_Brunt_Vaisala_aviation, EI_H2O_aviation
     316
     317IMPLICIT NONE
     318
     319!
     320! Input
     321!
     322REAL, INTENT(IN)  :: temp         ! air temperature [K]
     323REAL, INTENT(IN)  :: icesat_ratio ! ice saturation ratio [-]
     324REAL, INTENT(IN)  :: rho_air      ! air density [kg/m3]
     325REAL, INTENT(IN)  :: flight_dist  ! flown distance [m/s/m3]
     326REAL, INTENT(IN)  :: flight_fuel  ! fuel consumed [kg/s/m3]
     327!
     328! Output
     329!
     330REAL, INTENT(OUT) :: Nice_per_m_init_contrails ! [#/m]
     331REAL, INTENT(OUT) :: AEI_contrails             ! [#/kg]
     332REAL, INTENT(OUT) :: AEI_surv_contrails        ! [#/kg]
     333REAL, INTENT(OUT) :: fsurv_contrails           ! [-]
     334!
     335! Local
     336!
     337! For initial ice nucleation
     338REAL :: log_liqsat_ratio, coef_expo, dil_coef
     339REAL :: rd_act_amb, rd_act_soot, phi_amb, phi_soot
     340REAL :: AEI_soot, AEI_amb
     341! For ice crystals loss
     342REAL :: rho_emit, temp_205, nice_init
     343REAL :: z_atm, z_emit, z_desc, z_delta
     344REAL :: Nice_per_m, fuel_per_m, frac_surv
     345
     346!------------------------------
     347!--  INITIAL ICE NUCLEATION  --
     348!------------------------------
     349!
     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
     375AEI_contrails = EI_soot_aviation
     376
     377
     378!----------------------------
     379!--    ICE CRYSTALS LOSS   --
     380!----------------------------
     381!
     382!--Based on Lottermoser and Unterstrasser (2025, LU25 hereinafter)
     383!--which is an update of Unterstrasser (2016, U16 hereinafter)
     384
     385! fuel consumption in kg/m flown
     386fuel_per_m = flight_fuel / flight_dist
     387
     388! LU25, Eq. A2
     389z_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
     393rho_emit = fuel_per_m * EI_H2O_aviation / plume_area_loss
     394! LU25, Eq. A3
     395temp_205 = temp - 205.
     396z_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
     400z_desc = SQRT(8. * circ_0_loss / RPI / N_Brunt_Vaisala_aviation )
     401
     402! initial number of ice crystals [#/m flown], LU25, Eq. 1
     403Nice_per_m = fuel_per_m * AEI_contrails
     404! ice crystals number concentration [#/m3], LU25, Eq. A1
     405nice_init = Nice_per_m / plume_area_loss
     406! LU25, Eq. 9, 13d, 13e, 13f and 13g
     407z_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
     410fsurv_contrails = 0.42 + 1.31 / RPI * ATAN(-1. + z_delta / 100.)
     411fsurv_contrails = MIN(1., MAX(0., fsurv_contrails))
     412
     413Nice_per_m_init_contrails = Nice_per_m * fsurv_contrails
     414AEI_surv_contrails = AEI_contrails * fsurv_contrails
     415
     416END SUBROUTINE CONTRAIL_ICE_NUMBER_CONCENTRATION
    247417
    248418
    … …  
    258428      CALL xios_set_file_attr("aviation_file", enabled=.TRUE.)
    259429      ! Activate aviation fields
    260       CALL xios_set_field_attr("KMFLOWN_read", enabled=.TRUE.)
    261       CALL xios_set_field_attr("KMFLOWN_interp", enabled=.TRUE.)
     430      CALL xios_set_field_attr("DISTFLOWN_read", enabled=.TRUE.)
     431      CALL xios_set_field_attr("DISTFLOWN_interp", enabled=.TRUE.)
     432      CALL xios_set_field_attr("FUELBURN_read", enabled=.TRUE.)
     433      CALL xios_set_field_attr("FUELBURN_interp", enabled=.TRUE.)
    262434    ENDIF
    263435  ENDIF
    … …  
    286458    ! Local variable
    287459    !----------------------------------------------------
    288     REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:), flight_h2o_mpi(:,:,:)
     460    REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:), flight_fuel_mpi(:,:,:)
    289461    INTEGER :: ierr
    290462    LOGICAL, SAVE :: first = .TRUE.
    … …  
    297469    REAL                  :: npole, spole
    298470    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_dist_glo2D
    299     REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_h2o_glo2D
     471    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: flight_fuel_glo2D
    300472    REAL, ALLOCATABLE, DIMENSION(:)       :: vartmp_lev
    301473    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: vartmp
    … …  
    318490      ALLOCATE(flight_dist_read(klon, nleva, 1), STAT=ierr)
    319491      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_dist',1)
    320       ALLOCATE(flight_h2o_read(klon, nleva, 1), STAT=ierr)
    321       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_h2o',1)
     492      ALLOCATE(flight_fuel_read(klon, nleva, 1), STAT=ierr)
     493      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate flight_fuel',1)
    322494      ALLOCATE(aviation_lev(nleva), STAT=ierr)
    323495      IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'problem to allocate aviation_lev',1)
    … …  
    336508        ALLOCATE(flight_dist_mpi(klon_mpi, nleva,1), STAT=ierr)
    337509        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist_mpi',1)
    338         ALLOCATE(flight_h2o_mpi(klon_mpi, nleva,1), STAT=ierr)
    339         IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_h2o_mpi',1)
    340         CALL xios_recv_field("KMFLOWN_interp", flight_dist_mpi(:,:,1))
    341         !CALL xios_recv_field("KGH2O_interp", flight_h2o_mpi(:,:,1))
    342         flight_h2o_mpi(:,:,:) = 0.
     510        ALLOCATE(flight_fuel_mpi(klon_mpi, nleva,1), STAT=ierr)
     511        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_fuel_mpi',1)
     512        CALL xios_recv_field("DISTFLOWN_interp", flight_dist_mpi(:,:,1))
     513        CALL xios_recv_field("FUELBURN_interp", flight_fuel_mpi(:,:,1))
    343514    ENDIF
    344515
    345516    ! Propagate to other OMP threads: flight_dist_mpi(klon_mpi,klev) to flight_dist(klon,klev)
    346517    CALL scatter_omp(flight_dist_mpi, flight_dist_read)
    347     CALL scatter_omp(flight_h2o_mpi, flight_h2o_read)
     518    CALL scatter_omp(flight_fuel_mpi, flight_fuel_read)
    348519
    349520ELSE
    … …  
    358529    ELSE
    359530      ALLOCATE(flight_dist_glo2D(0,0,0,0))
    360       ALLOCATE(flight_h2o_glo2D(0,0,0,0))
     531      ALLOCATE(flight_fuel_glo2D(0,0,0,0))
    361532    ENDIF
    362533
    … …  
    413584       
    414585     ! Allocate variables depending on the number of vertical levels
    415        ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_h2o_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr)
     586       ALLOCATE(flight_dist_glo2D(nbp_lon, nbp_lat, nleva, 1), flight_fuel_glo2D(nbp_lon, nbp_lat, nleva, 1), stat=ierr)
    416587       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 1',1)
    417588
    … …  
    423594
    424595       ! Get variable id
    425        CALL check_err( nf90_inq_varid(nid, 'seg_length_m', varid),"pb inq var seg_length_m" )
     596       CALL check_err( nf90_inq_varid(nid, 'seg_length', varid),"pb inq var seg_length_m" )
    426597       ! Get the variable
    427598       CALL check_err( nf90_get_var(nid, varid, flight_dist_glo2D),"pb get var seg_length_m"  )
    428599
    429 !       ! Get variable id
    430 !       CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" )
    431 !       ! Get the variable
    432 !       CALL check_err( nf90_get_var(nid, varid, flight_h2o_glo2D),"pb get var fuel_burn"  )
    433        flight_h2o_glo2D(:,:,:,:) = 0.
     600       ! Get variable id
     601       CALL check_err( nf90_inq_varid(nid, 'fuel_burn', varid),"pb inq var fuel_burn" )
     602       ! Get the variable
     603       CALL check_err( nf90_get_var(nid, varid, flight_fuel_glo2D),"pb get var fuel_burn"  )
    434604
    435605       ! Get variable id
    … …  
    462632          END DO
    463633
    464           ! Inverse vertical levels for flight_h2o_glo2D
    465           vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1)
     634          ! Inverse vertical levels for flight_fuel_glo2D
     635          vartmp(:,:,:) = flight_fuel_glo2D(:,:,:,1)
    466636          DO k=1, nleva
    467637             DO j=1, nbp_lat
    468638                DO i=1, nbp_lon
    469                    flight_h2o_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
     639                   flight_fuel_glo2D(i,j,k,1) = vartmp(i,j,nleva+1-k)
    470640                END DO
    471641             END DO
    … …  
    502672
    503673          ! Invert latitudes for the variable
    504           vartmp(:,:,:) = flight_h2o_glo2D(:,:,:,1) ! use varmth temporarly
     674          vartmp(:,:,:) = flight_fuel_glo2D(:,:,:,1) ! use varmth temporarly
    505675          DO k=1,nleva
    506676             DO j=1,nbp_lat
    507677                DO i=1,nbp_lon
    508                    flight_h2o_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
     678                   flight_fuel_glo2D(i,j,k,1) = vartmp(i,nbp_lat+1-j,k)
    509679                END DO
    510680             END DO
    … …  
    533703          spole=0.  ! South pole, j=nbp_lat       
    534704          DO i=1,nbp_lon
    535              npole = npole + flight_h2o_glo2D(i,1,k,1)
    536              spole = spole + flight_h2o_glo2D(i,nbp_lat,k,1)
     705             npole = npole + flight_fuel_glo2D(i,1,k,1)
     706             spole = spole + flight_fuel_glo2D(i,nbp_lat,k,1)
    537707          END DO
    538708          npole = npole/REAL(nbp_lon)
    539709          spole = spole/REAL(nbp_lon)
    540           flight_h2o_glo2D(:,1,    k,1) = npole
    541           flight_h2o_glo2D(:,nbp_lat,k,1) = spole
     710          flight_fuel_glo2D(:,1,    k,1) = npole
     711          flight_fuel_glo2D(:,nbp_lat,k,1) = spole
    542712       END DO
    543713     
    544        ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_h2o_mpi(klon_glo, nleva, 1), stat=ierr)
     714       ALLOCATE(flight_dist_mpi(klon_glo, nleva, 1), flight_fuel_mpi(klon_glo, nleva, 1), stat=ierr)
    545715       IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 3',1)
    546716     
    547717       ! Transform from 2D to 1D field
    548718       CALL grid2Dto1D_glo(flight_dist_glo2D, flight_dist_mpi)
    549        CALL grid2Dto1D_glo(flight_h2o_glo2D, flight_h2o_mpi)
     719       CALL grid2Dto1D_glo(flight_fuel_glo2D, flight_fuel_mpi)
    550720   
    551721    ELSE
    552         ALLOCATE(flight_dist_mpi(0,0,0), flight_h2o_mpi(0,0,0))
     722        ALLOCATE(flight_dist_mpi(0,0,0), flight_fuel_mpi(0,0,0))
    553723    END IF ! is_mpi_root .AND. is_omp_root
    554724
    … …  
    571741
    572742    ! Allocate space for output pointer variable at local process
    573     ALLOCATE(flight_dist_read(klon, nleva, 1), flight_h2o_read(klon, nleva, 1), stat=ierr)
     743    ALLOCATE(flight_dist_read(klon, nleva, 1), flight_fuel_read(klon, nleva, 1), stat=ierr)
    574744    IF (ierr /= 0) CALL abort_physic('lmdz_aviation', 'pb in allocation 5',1)
    575745
    576746    ! Scatter global field to local domain at local process
    577747    CALL scatter(flight_dist_mpi, flight_dist_read)
    578     CALL scatter(flight_h2o_mpi, flight_h2o_read)
     748    CALL scatter(flight_fuel_mpi, flight_fuel_read)
    579749
    580750ENDIF
    … …  
    582752END SUBROUTINE read_aviation_emissions
    583753
    584 SUBROUTINE vertical_interpolation_aviation(klon, klev, paprs, pplay, temp, flight_dist, flight_h2o)
     754SUBROUTINE vertical_interpolation_aviation(klon, klev, paprs, pplay, temp, flight_dist, flight_fuel)
    585755    ! This subroutine performs the vertical interpolation from the read data in aviation.nc
    586756    ! where there are nleva vertical levels described in aviation_lev to the klev levels or
    587757    ! the model.
    588758    ! flight_dist_read(klon,nleva) -> flight_dist(klon, klev)
    589     ! flight_h2o_read(klon,nleva) -> flight_h2o(klon, klev)
     759    ! flight_fuel_read(klon,nleva) -> flight_fuel(klon, klev)
    590760
    591761    USE lmdz_lscp_ini, ONLY: RD, RG
    … …  
    599769    REAL, INTENT(IN)    :: temp(klon, klev) ! temperature [K]
    600770    REAL, INTENT(OUT)   :: flight_dist(klon,klev) ! Aviation distance flown concentration [m/s/m3]
    601     REAL, INTENT(OUT)   :: flight_h2o(klon,klev)  ! Aviation H2O emitted concentration [kgH2O/s/m3]
     771    REAL, INTENT(OUT)   :: flight_fuel(klon,klev) ! Aviation fuel burn concentration [kg/s/m3]
    602772
    603773    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    613783    ! Initialisation
    614784    flight_dist(:,:) = 0.
    615     flight_h2o(:,:) = 0.
     785    flight_fuel(:,:) = 0.
    616786
    617787    ! Compute the array with the vertical interface
    … …  
    649819                 ! Vertical reprojection for each desired array
    650820                 flight_dist(i,k) = flight_dist(i,k) + zfrac * flight_dist_read(i,kori,1)
    651                  flight_h2o(i,k)  = flight_h2o(i,k) + zfrac * flight_h2o_read(i,kori,1)
     821                 flight_fuel(i,k) = flight_fuel(i,k) + zfrac * flight_fuel_read(i,kori,1)
    652822            ENDDO
    653823
    … …  
    661831            !--Normalisation with the cell thickness
    662832            flight_dist(i,k) = flight_dist(i,k) / dz
    663             flight_h2o(i,k) = flight_h2o(i,k) / dz
     833            flight_fuel(i,k) = flight_fuel(i,k) / dz
    664834           
    665835            !--Enhancement factor
    666836            flight_dist(i,k) = flight_dist(i,k) * aviation_coef
    667             flight_h2o(i,k) = flight_h2o(i,k) * aviation_coef
     837            flight_fuel(i,k) = flight_fuel(i,k) * aviation_coef
    668838        ENDDO
    669839    ENDDO

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

    r5641 r5790  
    1212    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
    1313    !--AB contrails
    14     lincontfra, circontfra, qice_lincont, qice_circont, pclc_nocont, &
     14    contfra, qice_cont, Nice_cont, pclc_nocont, &
    1515    pcltau_nocont, pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
    1616    xfiwp_nocont, xfiwc_nocont, reice_nocont)
    … …  
    9595
    9696  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
    97   REAL, INTENT(IN)  :: lincontfra(klon, klev)    ! linear contrails fraction [-]
    98   REAL, INTENT(IN)  :: circontfra(klon, klev)    ! contrail induced cirrus fraction [-]
    99   REAL, INTENT(IN)  :: qice_lincont(klon, klev)  ! linear contrails condensed water [kg/kg]
    100   REAL, INTENT(IN)  :: qice_circont(klon, klev)  ! contrail induced cirrus condensed water [kg/kg]
     97  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
     98  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
     99  REAL, INTENT(IN)  :: Nice_cont(klon, klev)     ! contrails condensed water [kg/kg]
    101100  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
    102101  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
    … …  
    139138    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
    140139    !--AB for contrails
    141     lincontfra, circontfra, qice_lincont, qice_circont, pclc_nocont, pcltau_nocont, &
     140    contfra, qice_cont, Nice_cont, pclc_nocont, pcltau_nocont, &
    142141    pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
    143142    xfiwp_nocont, xfiwc_nocont, reice_nocont)

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

    r5717 r5790  
    1111    icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, &
    1212    !--AB contrails
    13     lincontfra, circontfra, qice_lincont, qice_circont, pclc_nocont, pcltau_nocont, &
     13    contfra, qice_cont, Nice_cont, pclc_nocont, pcltau_nocont, &
    1414    pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
    1515    xfiwp_nocont, xfiwc_nocont, reice_nocont)
    … …  
    3333  USE lmdz_cloud_optics_prop_ini , ONLY : rei_coef, rei_min_temp
    3434  USE lmdz_cloud_optics_prop_ini , ONLY : zepsec, novlp, iflag_ice_thermo, ok_new_lscp
    35   USE lmdz_cloud_optics_prop_ini , ONLY : ok_plane_contrail
    36   USE lmdz_cloud_optics_prop_ini , ONLY : rei_linear_contrails, rei_cirrus_contrails
     35  USE lmdz_cloud_optics_prop_ini , ONLY : ok_plane_contrail, eff2vol_radius_contrails, rho_ice
    3736 
    3837
    … …  
    120119
    121120  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
    122   REAL, INTENT(IN)  :: lincontfra(klon, klev)    ! linear contrails fraction [-]
    123   REAL, INTENT(IN)  :: circontfra(klon, klev)    ! contrail induced cirrus fraction [-]
    124   REAL, INTENT(IN)  :: qice_lincont(klon, klev)  ! linear contrails condensed water [kg/kg]
    125   REAL, INTENT(IN)  :: qice_circont(klon, klev)  ! contrail induced cirrus condensed water [kg/kg]
     121  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
     122  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
     123  REAL, INTENT(IN)  :: Nice_cont(klon, klev)     ! contrails ice crystals concentration [#/kg]
    126124  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
    127125  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
    … …  
    174172  REAL zflwp_var, zfiwp_var
    175173  REAL d_rei_dt
    176   REAL pclc_lincont(klon,klev), pclc_circont(klon,klev)
    177   REAL rei_lincont, rei_circont
     174  REAL pclc_cont(klon,klev)
     175  REAL mice_cont, rei_cont
    178176
    179177
    … …  
    196194  xfiwc = 0.D0
    197195  !--AB
    198   IF (ok_plane_contrail) THEN
     196  IF ( ok_plane_contrail ) THEN
    199197    xfiwp_nocont = 0.D0
    200198    xfiwc_nocont = 0.D0
    … …  
    252250
    253251    !--AB
    254     IF (ok_plane_contrail) THEN
     252    IF ( ok_plane_contrail ) THEN
    255253      !--If contrails are activated, we diagnose iwc without contrails
    256254      DO k = 1, klev
    257255        DO i = 1, klon
    258           pclc_nocont(i,k) = MAX(0., pclc(i, k) - lincontfra(i, k) - circontfra(i, k))
    259           xfiwc_nocont(i, k) = MAX(0., xfiwc(i, k) - qice_lincont(i, k) - qice_circont(i, k))
     256          pclc_nocont(i,k) = MAX(0., pclc(i, k) - contfra(i, k))
     257          xfiwc_nocont(i, k) = MAX(0., xfiwc(i, k) - qice_cont(i, k))
    260258        ENDDO
    261259      ENDDO
    … …  
    575573        reice_nocont(i,k) = 0.
    576574        pclc_nocont(i,k) = 0.
    577         pclc_lincont(i,k) = 0.
    578         pclc_circont(i,k) = 0.
     575        pclc_cont(i,k) = 0.
    579576        pcltau_cont(i,k) = 0.
    580577        pclemi_cont(i,k) = 0.
    … …  
    683680          reice_nocont(i,k) = 1.
    684681          pclc_nocont(i,k) = 0.
    685           pclc(i,k) = lincontfra(i,k) + circontfra(i,k)
     682          pclc(i,k) = contfra(i,k)
    686683          pcltau_nocont(i, k) = 0.
    687684          pclemi_nocont(i, k) = 0.
    … …  
    690687
    691688
    692         IF ( lincontfra(i,k) .GT. zepsec ) THEN
    693           pclc_lincont(i,k) = lincontfra(i,k)
    694           rei_lincont = rei_linear_contrails * 1.E6
    695         ELSE
    696           pclc_lincont(i,k) = 0.
    697           rei_lincont = 1.
    698         ENDIF
    699 
    700 
    701         IF ( circontfra(i,k) .GT. zepsec ) THEN
    702           !--Everything is the same but with contrails
    703           pclc_circont(i,k) = circontfra(i,k)
    704 
    705           !tc = temp(i, k) - 273.15
    706           !rei_circont = d_rei_dt*tc + rei_max
    707           !IF (tc<=-81.4) rei_circont = rei_min
    708           rei_circont = rei_cirrus_contrails * 1.E6
    709 
    710         ELSE
    711           pclc_circont(i,k) = 0.
    712           rei_circont = 1.
    713         ENDIF
    714 
    715         IF ( MAX(lincontfra(i,k), circontfra(i,k)) .GT. zepsec ) THEN
    716 
    717           rei = ( rei_lincont * pclc_lincont(i,k) + rei_circont * pclc_circont(i,k) ) &
    718               / ( pclc_lincont(i,k) + pclc_circont(i,k) )
     689        IF ( contfra(i,k) .GT. zepsec ) THEN
     690          pclc_cont(i,k) = contfra(i,k)
     691          mice_cont = qice_cont(i,k) / MAX(1e-10, Nice_cont(i,k))
     692          !--Volumic radius (in meters)
     693          rei_cont = (mice_cont / rho_ice * 3. / 4. / RPI)**(1./3.)
     694          !--Effective radius (in microns)
     695          rei_cont = MIN(100., MAX(1., rei_cont / eff2vol_radius_contrails * 1e6))
    719696          zfiwp_var = 1000.*(xfiwc(i, k)-xfiwc_nocont(i, k))&
    720697                    / (pclc(i, k)-pclc_nocont(i, k))*rhodz(i, k)
    721698
    722           pcltau_cont(i, k) = zfiwp_var*(3.448E-03+2.431/rei)
     699          pcltau_cont(i, k) = zfiwp_var*(3.448E-03+2.431/rei_cont)
    723700          ! -- cloud infrared emissivity:
    724701          ! [the broadband infrared absorption coefficient is PARAMETERized
    725702          ! as a function of the effective cld droplet radius]
    726703          ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
    727           k_ice = k_ice0 + 1.0/rei
     704          k_ice = k_ice0 + 1.0/rei_cont
    728705          pclemi_cont(i, k) = 1.0 - exp(-df*k_ice*zfiwp_var)
    729 
    730706        ELSE
     707          pclc_cont(i,k) = 0.
     708          rei_cont = 1.
    731709          pcltau_cont(i, k) = 0.
    732710          pclemi_cont(i, k) = 0.
    733711        ENDIF
    734712
    735 
    736         rei = ( rei_lincont * pclc_lincont(i,k) + rei_circont * pclc_circont(i,k) &
    737             + reice_nocont(i, k) * pclc_nocont(i, k) ) &
    738             / ( pclc_lincont(i,k) + pclc_circont(i,k) + pclc_nocont(i,k) )
     713        rei = ( rei_cont * pclc_cont(i,k) + reice_nocont(i, k) * pclc_nocont(i, k) ) &
     714            / ( pclc_cont(i,k) + pclc_nocont(i,k) )
    739715
    740716        zflwp_var = 1000.*xflwc(i, k)/pclc(i, k)*rhodz(i, k)
    … …  
    756732      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
    757733      xfiwp_nocont(i) = xfiwp_nocont(i) + xfiwc_nocont(i, k)*rhodz(i, k)
     734
     735      !--We weight the optical properties with the cloud fractions
     736      !--This is only used for the diagnostics
     737      pcltau_nocont(i, k) = pcltau_nocont(i, k) * pclc_nocont(i,k)
     738      pclemi_nocont(i, k) = pclemi_nocont(i, k) * pclc_nocont(i,k)
     739      pcltau_cont(i, k) = pcltau_cont(i, k) * pclc_cont(i,k)
     740      pclemi_cont(i, k) = pclemi_cont(i, k) * pclc_cont(i,k)
     741      pcltau(i, k) = pcltau(i, k) * pclc(i,k)
     742      pclemi(i, k) = pclemi(i, k) * pclc(i,k)
    758743
    759744    ENDDO
    … …  
    883868    DO k = klev, 1, -1
    884869      DO i = 1, klon
    885         zclear(i) = zclear(i)*(1.-max(pclc_lincont(i,k)+pclc_circont(i,k),zcloud(i)))/&
     870        zclear(i) = zclear(i)*(1.-max(pclc_cont(i,k),zcloud(i)))/&
    886871          (1.-min(real(zcloud(i),kind=8),1.-zepsec))
    887872        pct_cont(i) = 1. - zclear(i)
    888         zcloud(i) = pclc_lincont(i,k) + pclc_circont(i,k)
     873        zcloud(i) = pclc_cont(i,k)
    889874        IF (paprs(i,k)<prmhc) THEN
    890875          pch_nocont(i) = pch_nocont(i)*(1.-max(pclc_nocont(i,k),zcloudh(i)))/(1.-min(real( &

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

    r5717 r5790  
    2424  REAL, PROTECTED :: rei_coef, rei_min_temp
    2525  REAL, PROTECTED :: zepsec
    26   REAL, PROTECTED :: rei_linear_contrails=7.5E-6 ! [m] effective radius of ice crystals in linear contrails
    27   REAL, PROTECTED :: rei_cirrus_contrails=10.E-6 ! [m] effective radius of ice crystals in contrails cirrus
     26  REAL, PROTECTED :: eff2vol_radius_contrails=0.7
     27  REAL, PROTECTED :: rho_ice=920. ! Ice crystal density (assuming spherical geometry) [kg/m3]
    2828  REAL, PARAMETER :: thres_tau=0.3, thres_neb=0.001
    2929  REAL, PARAMETER :: prmhc=440.*100. ! Pressure between medium and high level cloud in Pa
    3030  REAL, PARAMETER :: prlmc=680.*100. ! Pressure between low and medium level cloud in Pa
    3131  REAL, PARAMETER :: coef_froi=0.09, coef_chau =0.13
    32   REAL, PROTECTED :: seuil_neb=0.001
     32  REAL, PARAMETER :: seuil_neb=0.001
    3333! if iflag_t_glace=0, old values are used for liquid/ice partitionning:
    3434  REAL, PARAMETER :: t_glace_min_old = 258.
    … …  
    4444!$OMP THREADPRIVATE(rei_coef, rei_min_temp)
    4545!$OMP THREADPRIVATE(zepsec)
    46 !$OMP THREADPRIVATE(rei_linear_contrails, rei_cirrus_contrails, seuil_neb)
     46!$OMP THREADPRIVATE(eff2vol_radius_contrails, rho_ice)
    4747
    4848 
    … …  
    110110    rei_min_temp = 175.
    111111    CALL getin_p('rei_min_temp',rei_min_temp)
    112     CALL getin_p('rei_linear_contrails', rei_linear_contrails)
    113     write(lunout,*)'rei_linear_contrails=',rei_linear_contrails
    114     CALL getin_p('rei_cirrus_contrails', rei_cirrus_contrails)
    115     write(lunout,*)'rei_cirrus_contrails=',rei_cirrus_contrails
    116     CALL getin_p('seuil_neb_rad', seuil_neb)
    117     write(lunout,*)'seuil_neb_rad=',seuil_neb
     112    CALL getin_p('eff2vol_radius_contrails', eff2vol_radius_contrails)
     113    write(lunout,*)'eff2vol_radius_contrails=',eff2vol_radius_contrails
    118114
    119115   

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

    r5779 r5790  
    9797      cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, &
    9898      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, &
    105102      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, dqi_auto, &
    106103      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       dcfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, &
    111       dcfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix, &
    112       dcfl_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, dcfc_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)
    115112
    116113!----------------------------------------------------------------------
    … …  
    140137USE lmdz_lscp_ini,   ONLY: chi_sedim
    141138
    142 USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_lincontrails
    143 USE lmdz_lscp_ini,   ONLY: coef_mixing_lincontrails, coef_shear_lincontrails
    144 USE lmdz_lscp_ini,   ONLY: chi_mixing_lincontrails, linear_contrails_lifetime
    145 USE lmdz_lscp_ini,   ONLY: fallice_linear_contrails, fallice_cirrus_contrails
     139USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_contrails
     140USE lmdz_lscp_ini,   ONLY: coef_mixing_contrails, coef_shear_contrails
     141USE lmdz_lscp_ini,   ONLY: chi_mixing_contrails, eff2vol_radius_contrails
     142USE lmdz_lscp_tools, ONLY: ICECRYSTAL_VELO
    146143USE lmdz_aviation,   ONLY: contrails_formation
    147144
    … …  
    175172REAL,     INTENT(IN)   , DIMENSION(klon) :: flsed_abv     ! sedimentated ice flux FROM THE LAYER ABOVE [kg/s/m2]
    176173REAL,     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 [-]
    183174REAL,     INTENT(INOUT), DIMENSION(klon) :: dzsed         ! sedimentated cloud height [m]
    184175REAL,     INTENT(INOUT), DIMENSION(klon) :: flsed         ! sedimentated ice flux [kg/s/m2]
    185176REAL,     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 [-]
    192177REAL,     INTENT(INOUT), DIMENSION(klon) :: flauto        ! autoconverted ice flux [kg/s/m2]
    193178!
    194179! Input for aviation
    195180!
    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]
     181REAL,     INTENT(IN)   , DIMENSION(klon) :: contfra_in    ! input contrails fraction [-]
     182REAL,     INTENT(IN)   , DIMENSION(klon) :: qtc_in        ! input contrails total specific humidity [kg/kg]
     183REAL,     INTENT(IN)   , DIMENSION(klon) :: nic_in        ! input contrails ice crystals concentration [#/kg]
    200184REAL,     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]
     185REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_fuel   ! aviation fuel consumption concentration [kg/s/m3]
     186REAL,     INTENT(IN)   , DIMENSION(klon) :: dzsed_cont_abv! sedimentated contrails height IN THE LAYER ABOVE [m]
     187REAL,     INTENT(IN)   , DIMENSION(klon) :: flsed_cont_abv! sedimentated ice flux in contrails FROM THE LAYER ABOVE [kg/s/m2]
     188REAL,     INTENT(IN)   , DIMENSION(klon) :: Nflsed_cont_abv! sedimentated ice crystals concentration in contrails FROM THE LAYER ABOVE [#/s/m2]
     189REAL,     INTENT(IN)   , DIMENSION(klon) :: cfsed_cont_abv! contrails fraction IN THE LAYER ABOVE [-]
     190REAL,     INTENT(INOUT), DIMENSION(klon) :: dzsed_cont    ! sedimentated contrails height [m]
     191REAL,     INTENT(INOUT), DIMENSION(klon) :: flsed_cont    ! sedimentated ice flux in contrails [kg/s/m2]
     192REAL,     INTENT(INOUT), DIMENSION(klon) :: Nflsed_cont   ! sedimentated ice crystals concentration flux in contrails [#/s/m2]
     193REAL,     INTENT(INOUT), DIMENSION(klon) :: cfsed_cont    ! sedimentated contrails fraction [-]
    202194!
    203195!  Output
    … …  
    238230!  NB. idem for the INOUT
    239231!
    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]
     232REAL,     INTENT(INOUT), DIMENSION(klon) :: contfra      ! contrail fraction [-]
     233REAL,     INTENT(INOUT), DIMENSION(klon) :: qcont        ! contrail specific humidity [kg/kg]
     234REAL,     INTENT(INOUT), DIMENSION(klon) :: Ncont        ! contrail ice crystals concentration [#/kg]
    244235REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
    245236REAL,     INTENT(INOUT), DIMENSION(klon) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
    246237REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
    247238REAL,     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]
     239REAL,     INTENT(INOUT), DIMENSION(klon) :: AEI_contrails ! Apparent emission index contrails [#/kg]
     240REAL,     INTENT(INOUT), DIMENSION(klon) :: AEI_surv_contrails ! Apparent emission index contrails after vortex loss [#/kg]
     241REAL,     INTENT(INOUT), DIMENSION(klon) :: fsurv_contrails ! Survival fraction after vortex loss [-]
     242REAL,     INTENT(INOUT), DIMENSION(klon) :: section_contrails ! Cross section of newly formed contrails [m2]
     243REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_ini     ! contrails cloud fraction tendency because of initial formation [s-1]
     244REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_ini     ! contrails ice specific humidity tendency because of initial formation [kg/kg/s]
     245REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_ini     ! contrails total specific humidity tendency because of initial formation [kg/kg/s]
     246REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_ini     ! contrails ice crystals concentration tendency because of initial formation [#/kg/s]
     247REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_sub     ! contrails cloud fraction tendency because of sublimation [s-1]
     248REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_sub     ! contrails ice specific humidity tendency because of sublimation [kg/kg/s]
     249REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_sub     ! contrails total specific humidity tendency because of sublimation [kg/kg/s]
     250REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_sub     ! contrails ice crystals concentration tendency because of sublimation [#/kg/s]
     251REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_mix     ! contrails cloud fraction tendency because of mixing [s-1]
     252REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_mix     ! contrails ice specific humidity tendency because of mixing [kg/kg/s]
     253REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_mix     ! contrails total specific humidity tendency because of mixing [kg/kg/s]
     254REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_mix     ! contrails ice crystals concentration tendency because of mixing [#/kg/s]
     255REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_agg     ! contrails ice crystals concentration tendency because of aggregation [#/kg/s]
     256REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_sed     ! contrails cloud fraction tendency because of ice sedimentation [s-1]
     257REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_sed     ! contrails ice specific humidity tendency because of ice sedimentation [kg/kg/s]
     258REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_sed     ! contrails total specific humidity tendency because of ice sedimentation [kg/kg/s]
     259REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_sed     ! contrails ice crystals concentration tendency because of ice sedimentation [#/kg/s]
     260REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_auto    ! contrails cloud fraction tendency because of ice autoconversion [s-1]
     261REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_auto    ! contrails ice specific humidity tendency because of ice autoconversion [kg/kg/s]
     262REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_auto    ! contrails total specific humidity tendency because of ice autoconversion [kg/kg/s]
     263REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_auto    ! contrails ice crystals concentration tendency because of ice autoconversion [#/kg/s]
    277264!
    278265! Local
    … …  
    308295! for sedimentation and autoconversion
    309296REAL, DIMENSION(klon) :: qised_abv
    310 REAL :: iwc, icefall_velo, coef_sed, qice_sedim
     297REAL :: iwc, icefall_velo, coef_sed, qice_sedim, Nice_sedim
    311298REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3, sedfra_tmp
    312299REAL :: dcf_sed1, dcf_sed2, dqvc_sed1, dqvc_sed2, dqi_sed1, dqi_sed2
    313300REAL :: dz_auto, coef_auto, qice_auto
     301!
     302! for aggregation
     303REAL :: dei, sticking_coef
     304REAL :: gamma_agg, dispersion_coef
    314305!
    315306! for mixing
    … …  
    324315!
    325316! for cell properties
    326 REAL :: rho, rhodz, dz
     317REAL, DIMENSION(klon) :: rho
     318REAL :: rhodz, dz
    327319!
    328320! 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
     321REAL :: mice, Niceincld
     322REAL, DIMENSION(klon) :: qised_cont_abv, Nised_cont_abv
     323REAL :: dcfc_sed1, dcfc_sed2, dqtc_sed1, dnic_sed1, dqtc_sed2, dqic_sed1, dqic_sed2, dnic_sed2
     324REAL :: dz_auto_cont
    334325
    335326qzero(:) = 0.
    … …  
    365356    ENDIF
    366357    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.
    371360      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.
    378365      ENDIF
    379366    ENDIF
    … …  
    461448      !--Initialisation of the cell properties
    462449      !--Dry density [kg/m3]
    463       rho = pplay(i) / temp(i) / RD
     450      rho(i) = pplay(i) / temp(i) / RD
    464451      !--Dry air mass [kg/m2]
    465452      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
    466453      !--Cell thickness [m]
    467       dz = rhodz / rho
     454      dz = rhodz / rho(i)
    468455
    469456      !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
    … …  
    537524      !--NB. the greater kappa_depsub, the more efficient is the
    538525      !--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.) &
    541528                   / ( RV * temp(i) / water_vapor_diff / pres_sat &
    542529                     + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
    … …  
    544531      !--If contrails are activated
    545532      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))
    550533        !--The following barriers are needed since the advection scheme does not
    551534        !--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
    567542          !--If ice sedimentation is activated, the quantity of sedimentated ice was added
    568543          !--to the total water vapor in the precipitation routine. Here we remove it
    … …  
    571546          !--inacurracies in the advection scheme might lead to considering this water
    572547          !--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) &
    574549              / ( 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)
    576553        ENDIF
    577554
    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.
    605559        dcfc_sub(i) = 0.
    606560        dqic_sub(i) = 0.
    607561        dqtc_sub(i) = 0.
     562        dnic_sub(i) = 0.
    608563        dcfc_mix(i) = 0.
    609564        dqic_mix(i) = 0.
    610565        dqtc_mix(i) = 0.
     566        dnic_mix(i) = 0.
     567        dnic_agg(i) = 0.
    611568        dcfc_sed(i) = 0.
    612569        dqic_sed(i) = 0.
    613570        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.
    614576      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.
    619579      ENDIF
    620580
    … …  
    624584      !----------------------------------------------------------------------
    625585
    626       !--If there is a linear contrail
    627       IF ( lincontfra(i) .GT. eps ) THEN
     586      !--If there is a contrail
     587      IF ( contfra(i) .GT. eps ) THEN
    628588        !--We remove contrails from the main class
    629         cldfra(i) = MAX(0., cldfra(i) - lincontfra(i))
    630         qcld(i) = MAX(0., qcld(i) - qlincont(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)))
    632592
    633593        !--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)
    635597
    636598        !--If the ice water content is too low, the cloud is purely sublimated
    637599        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)
    641604       
    642605        !--Only a part of the contrail is sublimated
    … …  
    649612
    650613          !--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
    654618        ENDIF ! qiceincld .LT. eps
    655619
    656620        !--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
    699627
    700628
    … …  
    11181046      !--------------------------------------
    11191047
    1120       IF ( ( lincontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps ) ) THEN
     1048      IF ( ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps ) ) THEN
    11211049
    11221050        !-- PART 1 - TURBULENT DIFFUSION
    … …  
    11381066        !-- clouds_perim = P * N_cld_mix
    11391067        !--
    1140         bovera = aspect_ratio_lincontrails
     1068        bovera = aspect_ratio_contrails
    11411069        !--P / a is a function of b / a only, that we can calculate
    11421070        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
    … …  
    11491077        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
    11501078        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
    1151         a_mix = Povera / coef_mixing_lincontrails / RPI / ( 1. - lincontfra(i) ) / bovera
     1079        a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - contfra(i) ) / bovera
    11521080        !--and finally,
    11531081        !-- 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) ) &
    11551083                  * cell_area(i) / Povera / a_mix
    11561084       
    … …  
    11821110        !--With this, the clouds increase in size along b only, by a factor
    11831111        !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind)
    1184         L_shear = coef_shear_lincontrails * shear(i) * dz * dtime
     1112        L_shear = coef_shear_contrails * shear(i) * dz * dtime
    11851113        !--therefore, the fraction of clear sky mixed is
    11861114        !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area
    … …  
    11891117        !--(note that they are equal)
    11901118        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)
    11911120        !--and the environment and cloud mixed fractions are the same,
    11921121        !--which we add to the previous calculated mixed fractions.
    … …  
    12001129       
    12011130        clrfra_mix = MIN(clrfra(i), clrfra_mix)
    1202         cldfra_mix = MIN(lincontfra(i), cldfra_mix)
     1131        cldfra_mix = MIN(contfra(i), cldfra_mix)
    12031132       
    12041133        !--We compute the limit vapor in clear sky where the mixed cloud could not
    … …  
    12091138        !--cloud's vapor without condensing or sublimating ice crystals
    12101139        qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
    1211                       - qlincont(i) / lincontfra(i) * cldfra_mix / clrfra_mix
     1140                      - qcont(i) / contfra(i) * cldfra_mix / clrfra_mix
    12121141
    12131142        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) )
    12201147        ELSE
    12211148          !--We then calculate the clear sky part where the humidity is lower than
    … …  
    12531180          !--decreases the size of the clouds
    12541181          !--For aviation, we increase the chance that the air surrounding contrails
    1255           !--is moist. This is quantified with chi_mixing_lincontrails
     1182          !--is moist. This is quantified with chi_mixing_contrails
    12561183          sigma_mix = MAX(1. - pdf_fra_below_lim / clrfra_mix, &
    12571184                      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 )))
    14271187
    14281188          IF ( pdf_fra_above_lim .GT. eps ) THEN
    … …  
    14371197            dcfc_mix(i) = dcfc_mix(i) - cldfra_mix * ( 1. - sigma_mix )
    14381198            dqtc_mix(i) = dqtc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
    1439                         * qcircont(i) / circontfra(i)
     1199                        * qcont(i) / contfra(i)
    14401200            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) )
    14421204          ENDIF
    14431205
    14441206        ENDIF
    1445       ENDIF ! ( circontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )
    1446 
    1447       IF ( ( lincontfra(i) + circontfra(i) ) .GT. eps ) THEN
     1207      ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )
     1208
     1209      IF ( contfra(i) .GT. eps ) THEN
    14481210        !--We balance the mixing tendencies between the different cloud classes
    1449         mixed_fraction = dcf_mix(i) + dcfl_mix(i) + dcfc_mix(i)
     1211        mixed_fraction = dcf_mix(i) + dcfc_mix(i)
    14501212        IF ( mixed_fraction .GT. clrfra(i) ) THEN
    14511213          mixed_fraction = clrfra(i) / mixed_fraction
    … …  
    14531215          dqvc_mix(i) = dqvc_mix(i) * mixed_fraction
    14541216          dqi_mix(i)  = dqi_mix(i)  * mixed_fraction
    1455           dcfl_mix(i) = dcfl_mix(i) * mixed_fraction
    1456           dqtl_mix(i) = dqtl_mix(i) * mixed_fraction
    1457           dqil_mix(i) = dqil_mix(i) * mixed_fraction
    14581217          dcfc_mix(i) = dcfc_mix(i) * mixed_fraction
    14591218          dqtc_mix(i) = dqtc_mix(i) * mixed_fraction
    14601219          dqic_mix(i) = dqic_mix(i) * mixed_fraction
     1220          dnic_mix(i) = dnic_mix(i) * mixed_fraction
    14611221        ENDIF
    14621222
    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)
    14781229      ENDIF
    14791230
    … …  
    14861237
    14871238
     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
    14881274      IF ( ok_ice_sedim ) THEN
    14891275        !--Initialisation
    14901276        dz_auto = 0.
    1491         dz_auto_lincont = 0.
    1492         dz_auto_circont = 0.
     1277        dz_auto_cont = 0.
    14931278        dcf_sed1  = 0.
    14941279        dqvc_sed1 = 0.
    … …  
    14971282        dqvc_sed2 = 0.
    14981283        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.
    15051284        dcfc_sed1 = 0.
    15061285        dqtc_sed1 = 0.
    15071286        dqic_sed1 = 0.
     1287        dnic_sed1 = 0.
    15081288        dcfc_sed2 = 0.
    15091289        dqtc_sed2 = 0.
    15101290        dqic_sed2 = 0.
    1511         !---------------------------------------
    1512         !--         ICE SEDIMENTATION         --
    1513         !---------------------------------------
     1291        dnic_sed2 = 0.
    15141292        !
    15151293        !--First, the current ice is sedimentated into the layer below. The ice fallspeed
    … …  
    15191297        !
    15201298        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)
    15221300          icefall_velo = fallice_sedim * cice_velo * MAX(0., iwc)**dice_velo
    15231301
    … …  
    15551333        ENDIF
    15561334        !
    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))
    15591338
    15601339          !--Sedimentation
    15611340          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
    15651345          !--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
    15691350          !--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
    15711353
    15721354          !--Autoconversion
    15731355          coef_auto = MAX(0., MIN(1., icefall_velo * dtime / dz - 1.))
    1574           dz_auto_lincont = MIN(dz, icefall_velo * dtime) * coef_auto
    1575           qice_auto = ( qlincont(i) - qsat(i) * lincontfra(i) ) * dz_auto_lincont / dz
     1356          dz_auto_cont = MIN(dz, icefall_velo * dtime) * coef_auto
     1357          qice_auto = ( qcont(i) - qsat(i) * contfra(i) ) * dz_auto_cont / dz
    15761358          !--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
    15801363          !--Convert to flux
    15811364          flauto(i) = flauto(i) + qice_auto * ( paprsdn(i) - paprsup(i) ) / RG / dtime
    15821365
    15831366          !--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)
    16211370          clrfra(i)  = clrfra(i) - dcfc_sed1 - dcfc_auto(i)
    16221371          qclr(i)    = qclr(i) - dqtc_sed1 - dqtc_auto(i)
    … …  
    17071456            ENDIF ! clrfra(i) .GT. eps
    17081457
    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))
    17181461                dqic_sed2  = dqic_sed2 + sedfra_tmp * qiceincld
    17191462            ENDIF
    … …  
    17321475
    17331476        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)
    17351478          IF ( sedfra_abv .GT. eps ) THEN
    17361479
    … …  
    17441487            !--(3) we increase the cloud fraction but use in-cloud water vapor as
    17451488            !--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) / dz
    1749             sedfra3 = MIN(cfsed_lincont(i), cfsed_lincont_abv(i)) &
    1750                     * MIN(MIN(dz, dzsed_lincont_abv(i)), &
    1751                     dzsed_lincont(i) + dz_auto_lincont) / dz
     1489            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
    17521495            sedfra1 = sedfra_abv - sedfra3 - sedfra2
    17531496
    1754             qiceincld = qised_lincont_abv(i) / sedfra_abv
     1497            qiceincld = qised_cont_abv(i) / sedfra_abv
     1498            Niceincld = Nised_cont_abv(i) / sedfra_abv
    17551499
    17561500            !--For region (1)
    … …  
    17871531                        sedfra1 * chi_sedim * pdf_fra_above_lim &
    17881532                        / ( 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))
    18811534
    18821535                IF ( pdf_fra_above_lim .GT. eps ) THEN
    … …  
    18861539                  dqic_sed2 = dqic_sed2 + sedfra_tmp &
    18871540                                   * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
     1541                  dnic_sed2 = dnic_sed2 + sedfra_tmp * Niceincld
    18881542                ENDIF
    18891543              ENDIF ! clrfra(i) .GT. eps
    … …  
    18911545              !--If the sedimented ice falls into a natural cirrus, it increases its content
    18921546              IF ( cldfra(i) .GT. eps ) THEN
    1893                   sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - circontfra(i))
     1547                  sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - contfra(i))
    18941548                  dqi_sed2  = dqi_sed2 + sedfra_tmp * qiceincld
    18951549              ENDIF
    1896 
    1897               !--If the sedimented ice falls into a contrail cirrus, it increases its content
    1898               IF ( lincontfra(i) .GT. eps ) THEN
    1899                   sedfra_tmp = sedfra1 * lincontfra(i) / (totfra_in(i) - circontfra(i))
    1900                   dqil_sed2  = dqil_sed2 + sedfra_tmp * qiceincld
    1901               ENDIF
    19021550            ENDIF ! sedfra1 .GT. eps
    19031551
    19041552            !--For region (2)
    19051553            dqic_sed2 = dqic_sed2 + sedfra2 * qiceincld
     1554            dnic_sed2 = dnic_sed2 + sedfra2 * Niceincld
    19061555
    19071556            !--For region (3)
    1908             IF ( circontfra(i) .GT. eps ) THEN
     1557            IF ( contfra(i) .GT. eps ) THEN
    19091558              dcfc_sed2 = dcfc_sed2 + sedfra3
    19101559              dqtc_sed2 = dqtc_sed2 + sedfra3 * (qsat(i) + qiceincld)
    19111560              dqic_sed2 = dqic_sed2 + sedfra3 * qiceincld
     1561              dnic_sed2 = dnic_sed2 + sedfra3 * Niceincld
    19121562            ENDIF
    19131563          ENDIF ! sedfra_abv .GT. eps
    … …  
    19291579        if ( ok_plane_contrail ) THEN
    19301580          !--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
    19351586          !--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)
    19491588          !--Diagnostics
    19501589          dcfc_sed(i) = dcfc_sed1 + dcfc_sed2
    19511590          dqtc_sed(i) = dqtc_sed1 + dqtc_sed2
    19521591          dqic_sed(i) = dqic_sed1 + dqic_sed2
     1592          dnic_sed(i) = dnic_sed1 + dnic_sed2
    19531593        ENDIF
    19541594
    … …  
    19571597
    19581598      !--We put back contrails in the clouds class
    1959       IF ( ( lincontfra(i) + circontfra(i) ) .GT. 0. ) THEN
    1960         cldfra(i) = cldfra(i) + lincontfra(i) + circontfra(i)
    1961         qcld(i) = qcld(i) + qlincont(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)
    19631603      ENDIF
    19641604
    … …  
    20421682IF ( ok_plane_contrail ) THEN
    20431683
     1684  rho = pplay / temp / RD
     1685
    20441686  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, &
    20471689      keepgoing, pt_pron_clds, &
    20481690      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)
    20501693
    20511694  DO i = 1, klon
    20521695    IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN
    20531696
    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))
    20711699
    20721700      !--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)
    20851712
    20861713
    … …  
    20891716      !-------------------------------------------
    20901717
    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.
    21051725      ENDIF
    21061726     
    … …  
    21101730        qcld(i)   = 0.
    21111731        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.
    21161735        qincld(i) = qsat(i)
    21171736      ELSE
    … …  
    21191738      ENDIF
    21201739
    2121       cldfra(i) = MAX(cldfra(i), lincontfra(i) + circontfra(i))
    2122       qcld(i) = MAX(qcld(i), qlincont(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))
    21241743
    21251744      !--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
    21431749      dcfc_sub(i) = dcfc_sub(i) / dtime
    21441750      dqic_sub(i) = dqic_sub(i) / dtime
    21451751      dqtc_sub(i) = dqtc_sub(i) / dtime
     1752      dnic_sub(i) = dnic_sub(i) / dtime
    21461753      dcfc_mix(i) = dcfc_mix(i) / dtime
    21471754      dqic_mix(i) = dqic_mix(i) / dtime
    21481755      dqtc_mix(i) = dqtc_mix(i) / dtime
     1756      dnic_mix(i) = dnic_mix(i) / dtime
     1757      dnic_agg(i) = dnic_agg(i) / dtime
    21491758      dcfc_sed(i) = dcfc_sed(i) / dtime
    21501759      dqic_sed(i) = dqic_sed(i) / dtime
    21511760      dqtc_sed(i) = dqtc_sed(i) / dtime
     1761      dnic_sed(i) = dnic_sed(i) / dtime
    21521762      dcfc_auto(i) = dcfc_auto(i) / dtime
    21531763      dqic_auto(i) = dqic_auto(i) / dtime
    21541764      dqtc_auto(i) = dqtc_auto(i) / dtime
     1765      dnic_auto(i) = dnic_auto(i) / dtime
    21551766
    21561767    ENDIF ! keepgoing

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

    r5779 r5790  
    232232  !$OMP THREADPRIVATE(ok_plane_contrail)
    233233
    234   LOGICAL, SAVE, PROTECTED :: ok_precip_lincontrails=.TRUE.  ! if True, linear contrails can be autoconverted to snow
    235   !$OMP THREADPRIVATE(ok_precip_lincontrails)
    236 
    237   LOGICAL, SAVE, PROTECTED :: ok_precip_circontrails=.FALSE. ! if True, cirrus contrails can be autoconverted to snow
    238   !$OMP THREADPRIVATE(ok_precip_circontrails)
    239 
    240   REAL, SAVE, PROTECTED :: aspect_ratio_lincontrails=.1      ! [-] aspect ratio of linear contrails
    241   !$OMP THREADPRIVATE(aspect_ratio_lincontrails)
    242 
    243   REAL, SAVE, PROTECTED :: coef_mixing_lincontrails          ! [-] tuning coefficient for the linear contrails mixing process
    244   !$OMP THREADPRIVATE(coef_mixing_lincontrails)
    245  
    246   REAL, SAVE, PROTECTED :: coef_shear_lincontrails           ! [-] additional coefficient for the linear contrails shearing process (subprocess of the contrails mixing process)
    247   !$OMP THREADPRIVATE(coef_shear_lincontrails)
    248  
    249   REAL, SAVE, PROTECTED :: chi_mixing_lincontrails=3.        ! [-] factor for increasing the chance that moist air is surrounding linear contrails
    250   !$OMP THREADPRIVATE(chi_mixing_lincontrails)
     234  REAL, SAVE, PROTECTED :: aspect_ratio_contrails            ! [-] aspect ratio of contrails
     235  !$OMP THREADPRIVATE(aspect_ratio_contrails)
     236
     237  REAL, SAVE, PROTECTED :: coef_mixing_contrails             ! [-] tuning coefficient for the contrails mixing process
     238  !$OMP THREADPRIVATE(coef_mixing_contrails)
     239 
     240  REAL, SAVE, PROTECTED :: coef_shear_contrails              ! [-] additional coefficient for the contrails shearing process (subprocess of the contrails mixing process)
     241  !$OMP THREADPRIVATE(coef_shear_contrails)
     242 
     243  REAL, SAVE, PROTECTED :: chi_mixing_contrails              ! [-] factor for increasing the chance that moist air is surrounding contrails
     244  !$OMP THREADPRIVATE(chi_mixing_contrails)
     245
     246  REAL, SAVE, PROTECTED :: Nice_init_contrails=100.          ! [#/cm3] initial ice crystals concentration in contrails
     247  !$OMP THREADPRIVATE(Nice_init_contrails)
     248
     249  REAL, SAVE, PROTECTED :: N_Brunt_Vaisala_aviation=0.01     ! [s-1] average Brunt Vaisala frequency, for contrail formation
     250  !$OMP THREADPRIVATE(N_Brunt_Vaisala_aviation)
     251
     252  REAL, SAVE, PROTECTED :: circ_0_loss                       ! [m/s] circulation for contrail formation, will be set automatically
     253  !$OMP THREADPRIVATE(circ_0_loss)
     254
     255  REAL, SAVE, PROTECTED :: plume_area_loss                   ! [m2] plume area for contrail formation, will be set automatically
     256  !$OMP THREADPRIVATE(plume_area_loss)
     257
     258  REAL, SAVE, PROTECTED :: nice_init_ref_loss                ! [#/m3] reference ice crystals concentration for contrail formation, will be set automatically
     259  !$OMP THREADPRIVATE(nice_init_ref_loss)
     260
     261  REAL, SAVE, PROTECTED :: Naer_amb=600.                     ! [#/cm3] ambiant background aerosol number concentration, for contrail formation
     262  !$OMP THREADPRIVATE(Naer_amb)
     263
     264  REAL, SAVE, PROTECTED :: raer_amb_mean=15.                 ! [nm] average radius of the background aerosols, for contrail formation
     265  !$OMP THREADPRIVATE(raer_amb_mean)
     266
     267  REAL, SAVE, PROTECTED :: raer_amb_std=2.2                  ! [?] standard deviation of the radius of the background aerosols, for contrail formation
     268  !$OMP THREADPRIVATE(raer_amb_std)
     269
     270  REAL, SAVE, PROTECTED :: r_soot_mean=15.                   ! [nm] average radius of the background aerosols, for contrail formation
     271  !$OMP THREADPRIVATE(r_soot_mean)
     272
     273  REAL, SAVE, PROTECTED :: r_soot_std=1.6                    ! [?] standard deviation of the radius of the background aerosols, for contrail formation
     274  !$OMP THREADPRIVATE(r_soot_std)
     275
     276  REAL, SAVE, PROTECTED :: air_to_fuel_ratio_engine=70.      ! [-] air to fuel ratio engine, for contrail formation
     277  !$OMP THREADPRIVATE(air_to_fuel_ratio_engine)
     278
     279  REAL, SAVE, PROTECTED :: wingspan=60.                      ! [m] average aircraft wingspan, for contrail formation
     280  !$OMP THREADPRIVATE(wingspan)
     281
     282  REAL, SAVE, PROTECTED :: EI_soot_aviation=1.5e15           ! [#/kg] emission index of soot number for a given fuel type
     283  !$OMP THREADPRIVATE(EI_soot_aviation)
    251284
    252285  REAL, SAVE, PROTECTED :: EI_H2O_aviation=1.25              ! [kgH2O/kg] emission index of water vapor for a given fuel type
    … …  
    273306  REAL, SAVE, PROTECTED :: fallice_cirrus_contrails=1.       ! [m/s] Ice fallspeed velocity in cirrus contrails
    274307  !$OMP THREADPRIVATE(fallice_cirrus_contrails)
     308
     309  REAL, SAVE, PROTECTED :: eff2vol_radius_contrails=0.7      ! [-]
     310  !$OMP THREADPRIVATE(eff2vol_radius_contrails)
    275311
    276312  REAL, SAVE, PROTECTED :: aviation_coef=1.                  ! [-] scaling factor for aviation emissions and flown distance
    … …  
    568604    CALL getin_p('chi_mixing',chi_mixing)
    569605    ! for aviation
    570     CALL getin_p('ok_precip_lincontrails',ok_precip_lincontrails)
    571     CALL getin_p('ok_precip_circontrails',ok_precip_circontrails)
    572     CALL getin_p('aspect_ratio_lincontrails',aspect_ratio_lincontrails)
    573     coef_mixing_lincontrails=coef_mixing_lscp
    574     CALL getin_p('coef_mixing_lincontrails',coef_mixing_lincontrails)
    575     coef_shear_lincontrails=coef_shear_lscp
    576     CALL getin_p('coef_shear_lincontrails',coef_shear_lincontrails)
    577     CALL getin_p('chi_mixing_lincontrails',chi_mixing_lincontrails)
     606    aspect_ratio_contrails=aspect_ratio_cirrus
     607    CALL getin_p('aspect_ratio_contrails',aspect_ratio_contrails)
     608    coef_mixing_contrails=coef_mixing_lscp
     609    CALL getin_p('coef_mixing_contrails',coef_mixing_contrails)
     610    coef_shear_contrails=coef_shear_lscp
     611    CALL getin_p('coef_shear_contrails',coef_shear_contrails)
     612    chi_mixing_contrails=chi_mixing
     613    CALL getin_p('chi_mixing_contrails',chi_mixing_contrails)
     614    CALL getin_p('Nice_init_contrails',Nice_init_contrails)
    578615    CALL getin_p('EI_H2O_aviation',EI_H2O_aviation)
     616    CALL getin_p('EI_soot_aviation',EI_soot_aviation)
     617    CALL getin_p('air_to_fuel_ratio_engine',air_to_fuel_ratio_engine)
     618    CALL getin_p('wingspan',wingspan)
     619    CALL getin_p('Naer_amb',Naer_amb)
     620    CALL getin_p('raer_amb_mean',raer_amb_mean)
     621    CALL getin_p('raer_amb_std',raer_amb_std)
     622    CALL getin_p('r_soot_mean',r_soot_mean)
     623    CALL getin_p('r_soot_std',r_soot_std)
     624    CALL getin_p('N_Brunt_Vaisala_aviation',N_Brunt_Vaisala_aviation)
    579625    CALL getin_p('qheat_fuel_aviation',qheat_fuel_aviation)
    580626    CALL getin_p('prop_efficiency_aviation',prop_efficiency_aviation)
    … …  
    584630    CALL getin_p('fallice_linear_contrails',fallice_linear_contrails)
    585631    CALL getin_p('fallice_cirrus_contrails',fallice_cirrus_contrails)
     632    CALL getin_p('eff2vol_radius_contrails',eff2vol_radius_contrails)
    586633    CALL getin_p('aviation_coef',aviation_coef)
    587634    ! for cloudth routines
    … …  
    689736    WRITE(lunout,*) 'lscp_ini, chi_mixing:', chi_mixing
    690737    ! for aviation
    691     WRITE(lunout,*) 'lscp_ini, ok_precip_lincontrails:', ok_precip_lincontrails
    692     WRITE(lunout,*) 'lscp_ini, ok_precip_circontrails:', ok_precip_circontrails
    693     WRITE(lunout,*) 'lscp_ini, aspect_ratio_lincontrails:', aspect_ratio_lincontrails
    694     WRITE(lunout,*) 'lscp_ini, coef_mixing_lincontrails:', coef_mixing_lincontrails
    695     WRITE(lunout,*) 'lscp_ini, coef_shear_lincontrails:', coef_shear_lincontrails
    696     WRITE(lunout,*) 'lscp_ini, chi_mixing_lincontrails:', chi_mixing_lincontrails
     738    WRITE(lunout,*) 'lscp_ini, aspect_ratio_contrails:', aspect_ratio_contrails
     739    WRITE(lunout,*) 'lscp_ini, coef_mixing_contrails:', coef_mixing_contrails
     740    WRITE(lunout,*) 'lscp_ini, coef_shear_contrails:', coef_shear_contrails
     741    WRITE(lunout,*) 'lscp_ini, chi_mixing_contrails:', chi_mixing_contrails
     742    WRITE(lunout,*) 'lscp_ini, Nice_init_contrails:', Nice_init_contrails
    697743    WRITE(lunout,*) 'lscp_ini, EI_H2O_aviation:', EI_H2O_aviation
     744    WRITE(lunout,*) 'lscp_ini, EI_soot_aviation:', EI_soot_aviation
     745    WRITE(lunout,*) 'lscp_ini, air_to_fuel_ratio_engine:', air_to_fuel_ratio_engine
     746    WRITE(lunout,*) 'lscp_ini, wingspan:', wingspan
     747    WRITE(lunout,*) 'lscp_ini, Naer_amb:', Naer_amb
     748    WRITE(lunout,*) 'lscp_ini, raer_amb_mean:', raer_amb_mean
     749    WRITE(lunout,*) 'lscp_ini, raer_amb_std:', raer_amb_std
     750    WRITE(lunout,*) 'lscp_ini, r_soot_mean:', r_soot_mean
     751    WRITE(lunout,*) 'lscp_ini, r_soot_std:', r_soot_std
     752    WRITE(lunout,*) 'lscp_ini, N_Brunt_Vaisala_aviation:', N_Brunt_Vaisala_aviation
    698753    WRITE(lunout,*) 'lscp_ini, qheat_fuel_aviation:', qheat_fuel_aviation
    699754    WRITE(lunout,*) 'lscp_ini, prop_efficiency_aviation:', prop_efficiency_aviation
    … …  
    703758    WRITE(lunout,*) 'lscp_ini, fallice_linear_contrails:', fallice_linear_contrails
    704759    WRITE(lunout,*) 'lscp_ini, fallice_cirrus_contrails:', fallice_cirrus_contrails
     760    WRITE(lunout,*) 'lscp_ini, eff2vol_radius_contrails:', eff2vol_radius_contrails
    705761    WRITE(lunout,*) 'lscp_ini, aviation_coef:', aviation_coef
    706762    ! for cloudth routines
    … …  
    765821                      / nu_iwc_pdf_lscp**(1./3.)
    766822
     823    IF ( ok_plane_contrail ) THEN
     824      !--Calculated here to lighten calculations
     825      !--See Unterstrasser (2016), Lottermoser and Unterstrasser (2025), and lmdz_aviation routine
     826      ! U16, Eq. A16 and LU25, section A2
     827      plume_area_loss = 2. * RPI * (1.5 + 0.314 * wingspan)**2
     828      ! initial circulation [m/s], U16, Eq. A5
     829      circ_0_loss = 10. * wingspan - 70.
     830      ! reference ice crystals number concentration [#/m3], LU25, section A1
     831      nice_init_ref_loss = 3.38e12 / (2. * RPI * (1.5 + 0.314 * 60.3)**2)
     832    ENDIF
     833
    767834    !AA Temporary initialisation
    768835    a_tr_sca(1) = -0.5

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

    r5779 r5790  
    2929     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
    3030     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
    31      cfl_seri, cfc_seri, qtl_seri, qtc_seri,            &
    32      qice_lincont, qice_circont, flight_dist,           &
    33      flight_h2o, qradice_lincont, qradice_circont,      &
     31     cfc_seri, qtc_seri, nic_seri,                      &
     32     qice_cont, flight_dist, flight_fuel, qradice_cont, &
    3433     Tcritcont, qcritcont, potcontfraP, potcontfraNP,   &
    3534     cloudth_sth,                                       &
    … …  
    128127USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
    129128USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato, ok_ice_sedim
    130 USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_precip_lincontrails, ok_precip_circontrails
     129USE lmdz_lscp_ini, ONLY : ok_plane_contrail
    131130USE lmdz_lscp_ini, ONLY : ok_nodeep_lscp, ok_nodeep_lscp_rad
    132131USE lmdz_lscp_ini, ONLY : ok_lscp_mergecond, gamma_mixth
    … …  
    135134USE phys_local_var_mod, ONLY : issrfra100to150, issrfra150to200, issrfra200to250
    136135USE phys_local_var_mod, ONLY : issrfra250to300, issrfra300to400, issrfra400to500
    137 USE phys_local_var_mod, ONLY : dcfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub
    138 USE phys_local_var_mod, ONLY : dcfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix
    139 USE phys_local_var_mod, ONLY : dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix
    140 USE phys_local_var_mod, ONLY : dcfl_sed, dqil_sed, dqtl_sed, dcfc_sed, dqic_sed, dqtc_sed
    141 USE phys_local_var_mod, ONLY : dcfl_auto, dqil_auto, dqtl_auto, dcfc_auto, dqic_auto, dqtc_auto
     136USE phys_local_var_mod, ONLY : AEI_contrails, AEI_surv_contrails
     137USE phys_local_var_mod, ONLY : fsurv_contrails, section_contrails
     138USE phys_local_var_mod, ONLY : dcfc_ini, dqic_ini, dqtc_ini, dnic_ini
     139USE phys_local_var_mod, ONLY : dcfc_sub, dqic_sub, dqtc_sub, dnic_sub
     140USE phys_local_var_mod, ONLY : dcfc_mix, dqic_mix, dqtc_mix, dnic_mix, dnic_agg
     141USE phys_local_var_mod, ONLY : dcfc_sed, dqic_sed, dqtc_sed, dnic_sed
     142USE phys_local_var_mod, ONLY : dcfc_auto, dqic_auto, dqtc_auto, dnic_auto
    142143USE phys_local_var_mod, ONLY : dcf_auto, dqi_auto, dqvc_auto
    143144USE geometry_mod, ONLY: longitude_deg, latitude_deg
    … …  
    212213  ! INPUT/OUTPUT aviation
    213214  !--------------------------------------------------
    214   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfl_seri         ! linear contrails fraction [-]
    215   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfc_seri         ! contrail cirrus fraction [-]
    216   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qtl_seri         ! linear contrails total specific humidity [kg/kg]
    217   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qtc_seri         ! contrail cirrus total specific humidity [kg/kg]
    218   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
    219   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
     215  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfc_seri         ! contrail fraction [-]
     216  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qtc_seri         ! contrail total specific humidity [kg/kg]
     217  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: nic_seri         ! contrail ice crystals concentration [#/kg]
     218  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown concentration [m/s/m3]
     219  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_fuel      ! aviation fuel consumption concentration [kg/s/m3]
    220220
    221221  ! OUTPUT variables
    … …  
    279279  ! for contrails and aviation
    280280
    281   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_lincont   !--condensed water in linear contrails [kg/kg]
    282   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_circont   !--condensed water in contrail cirrus [kg/kg]
    283   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qradice_lincont!--condensed water in linear contrails used in the radiation scheme [kg/kg]
    284   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qradice_circont!--condensed water in contrail cirrus used in the radiation scheme [kg/kg]
     281  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_cont      !--condensed water in contrails [kg/kg]
     282  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qradice_cont   !--condensed water in contrails used in the radiation scheme [kg/kg]
    285283  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
    286284  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
    … …  
    362360  REAL :: delta_z, deepconv_coef
    363361  ! for contrails
    364   REAL, DIMENSION(klon) :: lincontfra, circontfra, qlincont, qcircont
     362  REAL, DIMENSION(klon) :: contfra, qcont, Ncont
    365363  REAL, DIMENSION(klon) :: totfra_in, qtot_in
    366364  LOGICAL, DIMENSION(klon) :: pt_pron_clds
    367   REAL, DIMENSION(klon) :: dzsed_lincont, flsed_lincont, cfsed_lincont
    368   REAL, DIMENSION(klon) :: dzsed_circont, flsed_circont, cfsed_circont
    369   REAL, DIMENSION(klon) :: dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv
    370   REAL, DIMENSION(klon) :: dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv
     365  REAL, DIMENSION(klon) :: dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont
     366  REAL, DIMENSION(klon) :: dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv
    371367  !--for Lamquin et al 2012 diagnostics
    372368  REAL, DIMENSION(klon) :: issrfra100to150UP, issrfra150to200UP, issrfra200to250UP
    … …  
    472468flsed(:)        = 0.
    473469flauto(:)       = 0.
    474 flsed_lincont(:)= 0.
    475 flsed_circont(:)= 0.
     470flsed_cont(:)   = 0.
     471Nflsed_cont(:)  = 0.
    476472pt_pron_clds(:) = .FALSE.
     473
     474!--Contrails
     475dcfc_ini(:,:)  = 0.
     476dqic_ini(:,:)  = 0.
     477dqtc_ini(:,:)  = 0.
     478dnic_ini(:,:)  = 0.
     479dcfc_sub(:,:)  = 0.
     480dqic_sub(:,:)  = 0.
     481dqtc_sub(:,:)  = 0.
     482dnic_sub(:,:)  = 0.
     483dcfc_mix(:,:)  = 0.
     484dqic_mix(:,:)  = 0.
     485dqtc_mix(:,:)  = 0.
     486dnic_mix(:,:)  = 0.
     487dnic_agg(:,:)  = 0.
     488dcfc_sed(:,:)  = 0.
     489dqic_sed(:,:)  = 0.
     490dqtc_sed(:,:)  = 0.
     491dnic_sed(:,:)  = 0.
     492dcfc_auto(:,:) = 0.
     493dqic_auto(:,:) = 0.
     494dqtc_auto(:,:) = 0.
     495dnic_auto(:,:) = 0.
     496Tcritcont(:,:) = missing_val
     497qcritcont(:,:) = missing_val
     498potcontfraP(:,:)  = missing_val
     499potcontfraNP(:,:) = missing_val
     500AEI_contrails(:,:) = missing_val
     501AEI_surv_contrails(:,:) = missing_val
     502fsurv_contrails(:,:) = missing_val
     503section_contrails(:,:) = missing_val
     504
    477505
    478506!--for Lamquin et al (2012) diagnostics
    … …  
    580608      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
    581609                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
    582                         zqvapclr, zqupnew, flsed, flsed_lincont, flsed_circont, &
     610                        zqvapclr, zqupnew, flsed, flsed_cont, &
    583611                        cldfra_in, qvc_in, qliq_in, qice_in, &
    584612                        zrfl, zrflclr, zrflcld, &
    … …  
    592620      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
    593621                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
    594                         flsed, flsed_lincont, flsed_circont, &
     622                        flsed, flsed_cont, &
    595623                        zrfl, zrflclr, zrflcld, &
    596624                        zifl, ziflclr, ziflcld, &
    … …  
    718746          IF ( ok_plane_contrail ) THEN
    719747            IF ( iftop ) THEN
    720               dzsed_lincont_abv(:) = 0.
    721               flsed_lincont_abv(:) = 0.
    722               cfsed_lincont_abv(:) = 0.
    723               dzsed_circont_abv(:) = 0.
    724               flsed_circont_abv(:) = 0.
    725               cfsed_circont_abv(:) = 0.
     748              dzsed_cont_abv(:) = 0.
     749              flsed_cont_abv(:) = 0.
     750              Nflsed_cont_abv(:)= 0.
     751              cfsed_cont_abv(:) = 0.
    726752            ELSE
    727               dzsed_lincont_abv(:) = dzsed_lincont(:)
    728               flsed_lincont_abv(:) = flsed_lincont(:)
    729               cfsed_lincont_abv(:) = cfsed_lincont(:)
    730               dzsed_circont_abv(:) = dzsed_circont(:)
    731               flsed_circont_abv(:) = flsed_circont(:)
    732               cfsed_circont_abv(:) = cfsed_circont(:)
     753              dzsed_cont_abv(:) = dzsed_cont(:)
     754              flsed_cont_abv(:) = flsed_cont(:)
     755              Nflsed_cont_abv(:)= Nflsed_cont(:)
     756              cfsed_cont_abv(:) = cfsed_cont(:)
    733757            ENDIF
    734             lincontfra(:)    = 0.
    735             circontfra(:)    = 0.
    736             qlincont(:)      = 0.
    737             qcircont(:)      = 0.
     758            contfra(:) = 0.
     759            qcont(:)   = 0.
     760            Ncont(:)   = 0.
     761            dzsed_cont(:) = 0.
     762            flsed_cont(:) = 0.
     763            Nflsed_cont(:)= 0.
     764            cfsed_cont(:) = 0.
    738765          ENDIF
    739766
    … …  
    747774            cfsed_abv(:) = cfsed(:)
    748775          ENDIF
     776          dzsed(:) = 0.
     777          flsed(:) = 0.
     778          cfsed(:) = 0.
     779          flauto(:) = 0.
    749780
    750781          DO i = 1, klon
    … …  
    796827                  !--If contrails are activated, their fraction is also reduced when deep
    797828                  !--convection is active
    798                   cfl_seri(i,k) = cfl_seri(i,k) * deepconv_coef
    799                   qtl_seri(i,k) = qtl_seri(i,k) * deepconv_coef
    800829                  cfc_seri(i,k) = cfc_seri(i,k) * deepconv_coef
    801830                  qtc_seri(i,k) = qtc_seri(i,k) * deepconv_coef
     831                  nic_seri(i,k) = nic_seri(i,k) * deepconv_coef
    802832                ENDIF
    803833              ENDIF
    … …  
    925955                        shear, tke_dissip(:,k), cell_area, Tbef, qtot_in, zqs, &
    926956                        gammasat, ratqs(:,k), keepgoing, pt_pron_clds, &
    927                         dzsed_abv, flsed_abv, cfsed_abv, &
    928                         dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv, &
    929                         dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv, &
    930                         dzsed, flsed, cfsed, dzsed_lincont, flsed_lincont, cfsed_lincont, &
    931                         dzsed_circont, flsed_circont, cfsed_circont, flauto, &
     957                        dzsed_abv, flsed_abv, cfsed_abv, dzsed, flsed, cfsed, flauto, &
    932958                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
    933959                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
    … …  
    937963                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
    938964                        dqvc_sed(:,k), dqvc_auto(:,k), &
    939                         cfl_seri(:,k), cfc_seri(:,k), qtl_seri(:,k), qtc_seri(:,k), &
    940                         flight_dist(:,k), flight_h2o(:,k), &
    941                         lincontfra, circontfra, qlincont, qcircont, &
     965                        cfc_seri(:,k), qtc_seri(:,k), nic_seri(:,k), &
     966                        flight_dist(:,k), flight_fuel(:,k), &
     967                        contfra, qcont, Ncont, &
    942968                        Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
    943                         dcfl_ini(:,k), dqil_ini(:,k), dqtl_ini(:,k), &
    944                         dcfl_sub(:,k), dqil_sub(:,k), dqtl_sub(:,k), &
    945                         dcfl_cir(:,k), dqtl_cir(:,k), &
    946                         dcfl_mix(:,k), dqil_mix(:,k), dqtl_mix(:,k), &
    947                         dcfl_sed(:,k), dqil_sed(:,k), dqtl_sed(:,k), &
    948                         dcfl_auto(:,k), dqil_auto(:,k), dqtl_auto(:,k), &
    949                         dcfc_sub(:,k), dqic_sub(:,k), dqtc_sub(:,k), &
    950                         dcfc_mix(:,k), dqic_mix(:,k), dqtc_mix(:,k), &
    951                         dcfc_sed(:,k), dqic_sed(:,k), dqtc_sed(:,k), &
    952                         dcfc_auto(:,k), dqic_auto(:,k), dqtc_auto(:,k))
     969                        AEI_contrails(:,k), AEI_surv_contrails(:,k), &
     970                        fsurv_contrails(:,k), section_contrails(:,k), &
     971                        dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv, &
     972                        dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont, &
     973                        dcfc_ini(:,k), dqic_ini(:,k), dqtc_ini(:,k), dnic_ini(:,k), &
     974                        dcfc_sub(:,k), dqic_sub(:,k), dqtc_sub(:,k), dnic_sub(:,k), &
     975                        dcfc_mix(:,k), dqic_mix(:,k), dqtc_mix(:,k), dnic_mix(:,k), dnic_agg(:,k), &
     976                        dcfc_sed(:,k), dqic_sed(:,k), dqtc_sed(:,k), dnic_sed(:,k), &
     977                        dcfc_auto(:,k), dqic_auto(:,k), dqtc_auto(:,k), dnic_auto(:,k))
    953978
    954979                    IF ( ok_nodeep_lscp ) THEN
    … …  
    11091134
    11101135                        IF ( ok_ice_sedim ) THEN
    1111                           qice_sedim = (flauto(i) + flsed(i) + flsed_lincont(i) &
    1112                               + flsed_circont(i)) / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
     1136                          qice_sedim = (flauto(i) + flsed(i) + flsed_cont(i)) &
     1137                              / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
    11131138                          ! Add the ice that was sedimented, as it is not included in zqn
    11141139                          qlibef(i) = qlibef(i) + qice_sedim
    … …  
    11901215
    11911216                IF ( ok_ice_sedim ) THEN
    1192                   qice_sedim = (flauto(i) + flsed(i) + flsed_lincont(i) + flsed_circont(i)) &
     1217                  qice_sedim = (flauto(i) + flsed(i) + flsed_cont(i)) &
    11931218                      / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
    11941219                  ! Remove the ice that was sedimented. As it is not included in zqn,
    … …  
    12251250
    12261251      !--Ice water content of contrails
    1227       qice_lincont(:,k) = qlincont(:) - zqs(:) * lincontfra(:)
    1228       qice_circont(:,k) = qcircont(:) - zqs(:) * circontfra(:)
    1229 
    1230       !--Contrails precipitate as natural clouds. We save the partition of ice
    1231       !--between natural clouds and contrails
    1232       !--NB. we use qlincont / qcircont as a temporary variable to save this partition
    1233       IF ( ok_precip_lincontrails ) THEN
    1234         DO i = 1, klon
    1235           IF ( zoliqi(i) .GT. 0. ) THEN
    1236             qlincont(i) = qice_lincont(i,k) / zoliqi(i)
    1237           ELSE
    1238             qlincont(i) = 0.
    1239           ENDIF
    1240         ENDDO
    1241       ELSE
    1242         !--If linear contrails do not precipitate, they are removed temporarily from
    1243         !--the cloud variables
    1244         DO i = 1, klon
    1245           rneb(i,k) = rneb(i,k) - lincontfra(i)
    1246           zoliq(i) = zoliq(i) - qice_lincont(i,k)
    1247           zoliqi(i) = zoliqi(i) - qice_lincont(i,k)
    1248         ENDDO
    1249       ENDIF
    1250       IF ( ok_precip_circontrails ) THEN
    1251         DO i = 1, klon
    1252           IF ( zoliqi(i) .GT. 0. ) THEN
    1253             qcircont(i) = qice_circont(i,k) / zoliqi(i)
    1254           ELSE
    1255             qcircont(i) = 0.
    1256           ENDIF
    1257         ENDDO
    1258       ELSE
    1259         !--If contrails cirrus do not precipitate, they are removed temporarily from
    1260         !--the cloud variables
    1261         DO i = 1, klon
    1262           rneb(i,k) = rneb(i,k) - circontfra(i)
    1263           zoliq(i) = zoliq(i) - qice_circont(i,k)
    1264           zoliqi(i) = zoliqi(i) - qice_circont(i,k)
    1265         ENDDO
    1266       ENDIF
     1252      qice_cont(:,k) = qcont(:) - zqs(:) * contfra(:)
     1253
     1254      !--If contrails do not precipitate
     1255      DO i = 1, klon
     1256        rneb(i,k) = rneb(i,k) - contfra(i)
     1257        zoliq(i) = zoliq(i) - qice_cont(i,k)
     1258        zoliqi(i) = zoliqi(i) - qice_cont(i,k)
     1259      ENDDO
    12671260    ENDIF
    12681261
    … …  
    13061299
    13071300    IF ( ok_plane_contrail ) THEN
    1308       !--Contrails fraction is left unchanged, but contrails water has changed
    1309       !--We alse compute the ice content that will be seen by radiation
    1310       !--(qradice_lincont/circont)
    1311       IF ( ok_precip_circontrails ) THEN
    1312         DO i = 1, klon
    1313           IF ( zoliqi(i) .GT. 0. ) THEN
    1314             qradice_circont(i,k) = zradocond(i) * qcircont(i)
    1315             qcircont(i) = zqs(i) * circontfra(i) + zoliqi(i) * qcircont(i)
    1316           ELSE
    1317             qradice_circont(i,k) = 0.
    1318             circontfra(i) = 0.
    1319             qcircont(i) = 0.
    1320           ENDIF
    1321         ENDDO
    1322       ELSE
    1323         !--If contrails do not precipitate, they are put back into
    1324         !--the cloud variables
    1325         DO i = 1, klon
    1326           rneb(i,k) = rneb(i,k) + circontfra(i)
    1327           zoliq(i) = zoliq(i) + qice_circont(i,k)
    1328           zoliqi(i) = zoliqi(i) + qice_circont(i,k)
    1329           zradocond(i) = zradocond(i) + qice_circont(i,k)
    1330           zradoice(i) = zradoice(i) + qice_circont(i,k)
    1331           qradice_circont(i,k) = qice_circont(i,k)
    1332         ENDDO
    1333       ENDIF
    1334       IF ( ok_precip_lincontrails ) THEN
    1335         DO i = 1, klon
    1336           IF ( zoliqi(i) .GT. 0. ) THEN
    1337             qradice_lincont(i,k) = zradocond(i) * qlincont(i)
    1338             qlincont(i) = zqs(i) * lincontfra(i) + zoliqi(i) * qlincont(i)
    1339           ELSE
    1340             qradice_lincont(i,k) = 0.
    1341             lincontfra(i) = 0.
    1342             qlincont(i) = 0.
    1343           ENDIF
    1344         ENDDO
    1345       ELSE
    1346         !--If contrails do not precipitate, they are put back into
    1347         !--the cloud variables
    1348         DO i = 1, klon
    1349           rneb(i,k) = rneb(i,k) + lincontfra(i)
    1350           zoliq(i) = zoliq(i) + qice_lincont(i,k)
    1351           zoliqi(i) = zoliqi(i) + qice_lincont(i,k)
    1352           zradocond(i) = zradocond(i) + qice_lincont(i,k)
    1353           zradoice(i) = zradoice(i) + qice_lincont(i,k)
    1354           qradice_lincont(i,k) = qice_lincont(i,k)
    1355         ENDDO
    1356       ENDIF
     1301      !--Contrails do not precipitate
     1302      DO i = 1, klon
     1303        rneb(i,k) = rneb(i,k) + contfra(i)
     1304        zoliq(i) = zoliq(i) + qice_cont(i,k)
     1305        zoliqi(i) = zoliqi(i) + qice_cont(i,k)
     1306        zradocond(i) = zradocond(i) + qice_cont(i,k)
     1307        zradoice(i) = zradoice(i) + qice_cont(i,k)
     1308        qradice_cont(i,k) = qice_cont(i,k)
     1309      ENDDO
    13571310    ENDIF
    13581311
    … …  
    14531406      DO i = 1, klon
    14541407        IF ( zoliq(i) .LE. 0. ) THEN
    1455           lincontfra(i) = 0.
    1456           circontfra(i) = 0.
    1457           qlincont(i) = 0.
    1458           qcircont(i) = 0.
     1408          contfra(i) = 0.
     1409          qcont(i) = 0.
     1410          Ncont(i) = 0.
    14591411        ENDIF
    14601412      ENDDO
    1461       cfl_seri(:,k) = lincontfra(:)
    1462       cfc_seri(:,k) = circontfra(:)
    1463       qtl_seri(:,k) = qlincont(:)
    1464       qtc_seri(:,k) = qcircont(:)
     1413      cfc_seri(:,k) = contfra(:)
     1414      qtc_seri(:,k) = qcont(:)
     1415      nic_seri(:,k) = Ncont(:)
    14651416    ENDIF
    14661417
    … …  
    16191570  IF ( ok_ice_sedim ) THEN
    16201571    DO i = 1, klon
    1621       snow(i) = snow(i) + flsed(i) + flsed_lincont(i) + flsed_circont(i)
     1572      snow(i) = snow(i) + flsed(i) + flsed_cont(i)
    16221573    ENDDO
    16231574  ENDIF

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

    r5779 r5790  
    1818SUBROUTINE histprecip_precld( &
    1919           klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, &
    20            zmqc, zneb, znebprecip, znebprecipclr, flsed, flsed_lincont, flsed_circont, &
     20           zmqc, zneb, znebprecip, znebprecipclr, flsed, flsed_cont, &
    2121           zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
    2222           )
    … …  
    4545REAL,    INTENT(INOUT), DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
    4646REAL,    INTENT(IN),    DIMENSION(klon) :: flsed          !--sedimentated ice flux [kg/s/m2]
    47 REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_lincont  !--linear contrails sedimentated ice flux [kg/s/m2]
    48 REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_circont  !--contrail cirrus sedimentated ice flux [kg/s/m2]
     47REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_cont     !--contrails sedimentated ice flux [kg/s/m2]
    4948
    5049REAL,    INTENT(IN),    DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
    … …  
    134133    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
    135134    !--added to the total water content of the gridbox
    136     IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
    137       qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) &
     135    IF ( (flsed(i) + flsed_cont(i)) .GT. 0. ) THEN
     136      qice_sedim = (flsed(i) + flsed_cont(i)) &
    138137          / ( paprsdn(i) - paprsup(i) ) * RG * dtime
    139138
    … …  
    757756           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
    758757           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
    759            flsed, flsed_lincont, flsed_circont, &
     758           flsed, flsed_cont, &
    760759           cldfra, qvc, qliq, qice, &
    761760           rain, rainclr, raincld, snow, snowclr, snowcld, &
    … …  
    793792REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
    794793REAL,    INTENT(IN),    DIMENSION(klon) :: flsed          !--sedimentated ice water flux [kg/s/m2]
    795 REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_lincont  !--sedimentated ice water flux [kg/s/m2]
    796 REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_circont  !--sedimentated ice water flux [kg/s/m2]
     794REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_cont     !--sedimentated ice water flux [kg/s/m2]
    797795
    798796REAL,    INTENT(INOUT), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
    … …  
    896894  cpw = RCPD * RVTMP2
    897895  DO i = 1, klon
    898     IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
    899       qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) / dhum_to_dflux(i)
     896    IF ( (flsed(i) + flsed_cont(i)) .GT. 0. ) THEN
     897      qice_sedim = (flsed(i) + flsed_cont(i)) / dhum_to_dflux(i)
    900898      !--No condensed water so cp=cp(vapor+dry air)
    901899      !-- RVTMP2=rcpv/rcpd-1
    … …  
    10051003    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
    10061004    !--added to the total water content of the gridbox
    1007     IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
    1008       qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) / dhum_to_dflux(i)
     1005    IF ( (flsed(i) + flsed_cont(i)) .GT. 0. ) THEN
     1006      qice_sedim = (flsed(i) + flsed_cont(i)) / dhum_to_dflux(i)
    10091007      !--Vapor is updated after evaporation/sublimation (it is increased)
    10101008      qvap(i) = qvap(i) + qice_sedim

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

    r5717 r5790  
    108108
    109109END SUBROUTINE FALLICE_VELOCITY
     110!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     111
     112!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     113FUNCTION ICECRYSTAL_VELO(mice, temp, pres)
     114
     115    ! Ref:
     116    ! Spichtinger and Gierens (2009)
     117
     118    USE lmdz_lscp_ini, ONLY: RPI, rho_ice
     119
     120    IMPLICIT NONE
     121
     122    REAL :: mice      ! ice crystal mass [kg]
     123    REAL :: temp      ! temperature [K]
     124    REAL :: pres      ! air pressure [Pa]
     125    REAL :: icecrystal_velo    ! fallspeed velocity of crystals [m/s]
     126
     127    !--Local
     128    REAL :: c, g, d
     129
     130    c = (pres / 30000)**(-0.178) * (temp / 233.)**(-0.394)
     131    IF ( mice .LT. 2.146e-13 ) THEN
     132        g = 735.4
     133        d = 0.42
     134    ELSEIF ( mice .LT. 2.166e-9) THEN
     135        g = 63292.4
     136        d = 0.57
     137    ELSEIF ( mice .LT. 4.264e-8 ) THEN
     138        g = 329.8
     139        d = 0.31
     140    ELSE
     141        g = 8.8
     142        d = 0.096
     143    ENDIF
     144    icecrystal_velo = g * mice**d * c
     145
     146END FUNCTION ICECRYSTAL_VELO
    110147!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    111148

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

    r5717 r5790  
    2323       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
    2424       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, &
    25        cf_ancien, qvc_ancien, qvcon, qccon, cfl_ancien, cfc_ancien, qtl_ancien, qtc_ancien, &
     25       cf_ancien, qvc_ancien, qvcon, qccon, cfc_ancien, qtc_ancien, nic_ancien, &
    2626       radpas, radsol, rain_fall, &
    2727       ratqs, rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
    … …  
    482482  ! cas specifique des variables de l'aviation
    483483  IF ( ok_plane_contrail ) THEN
    484     ancien_ok=ancien_ok.AND.phyetat0_get(cfl_ancien,"CFLANCIEN","CFLANCIEN",0.)
    485484    ancien_ok=ancien_ok.AND.phyetat0_get(cfc_ancien,"CFCANCIEN","CFCANCIEN",0.)
    486     ancien_ok=ancien_ok.AND.phyetat0_get(qtl_ancien,"QTLANCIEN","QTLANCIEN",0.)
    487485    ancien_ok=ancien_ok.AND.phyetat0_get(qtc_ancien,"QTCANCIEN","QTCANCIEN",0.)
     486    ancien_ok=ancien_ok.AND.phyetat0_get(nic_ancien,"NICANCIEN","NICANCIEN",0.)
    488487  ELSE
    489     cfl_ancien(:,:)=0.
    490488    cfc_ancien(:,:)=0.
    491     qtl_ancien(:,:)=0.
    492489    qtc_ancien(:,:)=0.
     490    nic_ancien(:,:)=0.
    493491  ENDIF
    494492
    … …  
    521519
    522520  IF ( ok_plane_contrail ) THEN
    523     IF ( ( maxval(cfl_ancien).EQ.minval(cfl_ancien) ) .OR. &
    524          ( maxval(cfc_ancien).EQ.minval(cfc_ancien) ) .OR. &
    525          ( maxval(qtl_ancien).EQ.minval(qtl_ancien) ) .OR. &
    526          ( maxval(qtc_ancien).EQ.minval(qtc_ancien) ) ) THEN
     521    IF ( ( maxval(cfc_ancien).EQ.minval(cfc_ancien) ) .OR. &
     522         ( maxval(qtc_ancien).EQ.minval(qtc_ancien) ) .OR. &
     523         ( maxval(nic_ancien).EQ.minval(nic_ancien) ) ) THEN
    527524       ancien_ok=.false.
    528525     ENDIF

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

    r5717 r5790  
    2424                                ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
    2525                                qvcon, qccon,                                &
    26                                 qvc_ancien, cfl_ancien, cfc_ancien,          &
    27                                 qtl_ancien, qtc_ancien, u_ancien, v_ancien,  &
     26                                qvc_ancien, cfc_ancien, qtc_ancien, nic_ancien, &
     27                                u_ancien, v_ancien,                          &
    2828                                clwcon, rnebcon, ratqs, pbl_tke,             &
    2929                                wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
    … …  
    287287
    288288    IF ( ok_plane_contrail ) THEN
    289       CALL put_field(pass,"CFLANCIEN", "CFLANCIEN", cfl_ancien)
    290289      CALL put_field(pass,"CFCANCIEN", "CFCANCIEN", cfc_ancien)
    291       CALL put_field(pass,"QTLANCIEN", "QTLANCIEN", qtl_ancien)
    292290      CALL put_field(pass,"QTCANCIEN", "QTCANCIEN", qtc_ancien)
     291      CALL put_field(pass,"NICANCIEN", "NICANCIEN", nic_ancien)
    293292    ENDIF
    294293

  • LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90

    r5779 r5790  
    685685      REAL, SAVE, ALLOCATABLE :: d_q_avi(:,:)
    686686      !$OMP THREADPRIVATE(d_q_avi)
    687       REAL, SAVE, ALLOCATABLE :: cfl_seri(:,:), d_cfl_dyn(:,:)
    688       !$OMP THREADPRIVATE(cfl_seri, d_cfl_dyn)
    689687      REAL, SAVE, ALLOCATABLE :: cfc_seri(:,:), d_cfc_dyn(:,:)
    690688      !$OMP THREADPRIVATE(cfc_seri, d_cfc_dyn)
    691       REAL, SAVE, ALLOCATABLE :: qtl_seri(:,:), d_qtl_dyn(:,:)
    692       !$OMP THREADPRIVATE(qtl_seri, d_qtl_dyn)
    693689      REAL, SAVE, ALLOCATABLE :: qtc_seri(:,:), d_qtc_dyn(:,:)
    694690      !$OMP THREADPRIVATE(qtc_seri, d_qtc_dyn)
    695       REAL, SAVE, ALLOCATABLE :: flight_dist(:,:), flight_h2o(:,:)
    696       !$OMP THREADPRIVATE(flight_dist, flight_h2o)
     691      REAL, SAVE, ALLOCATABLE :: nic_seri(:,:), d_nic_dyn(:,:)
     692      !$OMP THREADPRIVATE(nic_seri, d_nic_dyn)
     693      REAL, SAVE, ALLOCATABLE :: flight_dist(:,:), flight_fuel(:,:)
     694      !$OMP THREADPRIVATE(flight_dist, flight_fuel)
    697695      REAL, SAVE, ALLOCATABLE :: Tcritcont(:,:), qcritcont(:,:)
    698696      !$OMP THREADPRIVATE(Tcritcont, qcritcont)
    699697      REAL, SAVE, ALLOCATABLE :: potcontfraP(:,:), potcontfraNP(:,:)
    700698      !$OMP THREADPRIVATE(potcontfraP, potcontfraNP)
    701       REAL, SAVE, ALLOCATABLE :: qice_lincont(:,:), qice_circont(:,:)
    702       !$OMP THREADPRIVATE(qice_lincont, qice_circont)
    703       REAL, SAVE, ALLOCATABLE :: qradice_lincont(:,:), qradice_circont(:,:)
    704       !$OMP THREADPRIVATE(qradice_lincont, qradice_circont)
    705       REAL, SAVE, ALLOCATABLE :: dcfl_cir(:,:), dqtl_cir(:,:)
    706       !$OMP THREADPRIVATE(dcfl_cir, dqtl_cir)
    707       REAL, SAVE, ALLOCATABLE :: dcfl_ini(:,:), dqil_ini(:,:), dqtl_ini(:,:)
    708       !$OMP THREADPRIVATE(dcfl_ini, dqil_ini, dqtl_ini)
    709       REAL, SAVE, ALLOCATABLE :: dcfl_sub(:,:), dqil_sub(:,:), dqtl_sub(:,:)
    710       !$OMP THREADPRIVATE(dcfl_sub, dqil_sub, dqtl_sub)
    711       REAL, SAVE, ALLOCATABLE :: dcfl_mix(:,:), dqil_mix(:,:), dqtl_mix(:,:)
    712       !$OMP THREADPRIVATE(dcfl_mix, dqil_mix, dqtl_mix)
    713       REAL, SAVE, ALLOCATABLE :: dcfl_sed(:,:), dqil_sed(:,:), dqtl_sed(:,:)
    714       !$OMP THREADPRIVATE(dcfl_sed, dqil_sed, dqtl_sed)
    715       REAL, SAVE, ALLOCATABLE :: dcfc_sed(:,:), dqic_sed(:,:), dqtc_sed(:,:)
    716       !$OMP THREADPRIVATE(dcfc_sed, dqic_sed, dqtc_sed)
    717       REAL, SAVE, ALLOCATABLE :: dcfl_auto(:,:), dqil_auto(:,:), dqtl_auto(:,:)
    718       !$OMP THREADPRIVATE(dcfl_auto, dqil_auto, dqtl_auto)
    719       REAL, SAVE, ALLOCATABLE :: dcfc_auto(:,:), dqic_auto(:,:), dqtc_auto(:,:)
    720       !$OMP THREADPRIVATE(dcfc_auto, dqic_auto, dqtc_auto)
    721       REAL, SAVE, ALLOCATABLE :: dcfc_sub(:,:), dqic_sub(:,:), dqtc_sub(:,:)
    722       !$OMP THREADPRIVATE(dcfc_sub, dqic_sub, dqtc_sub)
    723       REAL, SAVE, ALLOCATABLE :: dcfc_mix(:,:), dqic_mix(:,:), dqtc_mix(:,:)
    724       !$OMP THREADPRIVATE(dcfc_mix, dqic_mix, dqtc_mix)
     699      REAL, SAVE, ALLOCATABLE :: AEI_contrails(:,:), AEI_surv_contrails(:,:)
     700      !$OMP THREADPRIVATE(AEI_contrails, AEI_surv_contrails)
     701      REAL, SAVE, ALLOCATABLE :: fsurv_contrails(:,:), section_contrails(:,:)
     702      !$OMP THREADPRIVATE(fsurv_contrails, section_contrails)
     703      REAL, SAVE, ALLOCATABLE :: qice_cont(:,:)
     704      !$OMP THREADPRIVATE(qice_cont)
     705      REAL, SAVE, ALLOCATABLE :: qradice_cont(:,:)
     706      !$OMP THREADPRIVATE(qradice_cont)
     707      REAL, SAVE, ALLOCATABLE :: dcfc_ini(:,:), dqic_ini(:,:), dqtc_ini(:,:), dnic_ini(:,:)
     708      !$OMP THREADPRIVATE(dcfc_ini, dqic_ini, dqtc_ini, dnic_ini)
     709      REAL, SAVE, ALLOCATABLE :: dcfc_sed(:,:), dqic_sed(:,:), dqtc_sed(:,:), dnic_sed(:,:)
     710      !$OMP THREADPRIVATE(dcfc_sed, dqic_sed, dqtc_sed, dnic_sed)
     711      REAL, SAVE, ALLOCATABLE :: dcfc_auto(:,:), dqic_auto(:,:), dqtc_auto(:,:), dnic_auto(:,:)
     712      !$OMP THREADPRIVATE(dcfc_auto, dqic_auto, dqtc_auto, dnic_auto)
     713      REAL, SAVE, ALLOCATABLE :: dcfc_sub(:,:), dqic_sub(:,:), dqtc_sub(:,:), dnic_sub(:,:)
     714      !$OMP THREADPRIVATE(dcfc_sub, dqic_sub, dqtc_sub, dnic_sub)
     715      REAL, SAVE, ALLOCATABLE :: dcfc_mix(:,:), dqic_mix(:,:), dqtc_mix(:,:), dnic_mix(:,:)
     716      !$OMP THREADPRIVATE(dcfc_mix, dqic_mix, dqtc_mix, dnic_mix)
     717      REAL, SAVE, ALLOCATABLE :: dnic_agg(:,:)
     718      !$OMP THREADPRIVATE(dnic_agg)
    725719      REAL, SAVE, ALLOCATABLE :: cldfra_nocont(:,:), cldtau_nocont(:,:), cldemi_nocont(:,:)
    726720      !$OMP THREADPRIVATE(cldfra_nocont, cldtau_nocont, cldemi_nocont)
    … …  
    12991293!-- LSCP - aviation and contrails variables
    13001294      ALLOCATE(d_q_avi(klon,klev))
    1301       ALLOCATE(cfl_seri(klon,klev), d_cfl_dyn(klon,klev))
    13021295      ALLOCATE(cfc_seri(klon,klev), d_cfc_dyn(klon,klev))
    1303       ALLOCATE(qtl_seri(klon,klev), d_qtl_dyn(klon,klev))
    13041296      ALLOCATE(qtc_seri(klon,klev), d_qtc_dyn(klon,klev))
    1305       ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
     1297      ALLOCATE(nic_seri(klon,klev), d_nic_dyn(klon,klev))
     1298      ALLOCATE(flight_dist(klon,klev), flight_fuel(klon,klev))
    13061299      ALLOCATE(Tcritcont(klon,klev), qcritcont(klon,klev))
    13071300      ALLOCATE(potcontfraP(klon,klev), potcontfraNP(klon,klev))
    1308       ALLOCATE(qice_lincont(klon,klev), qice_circont(klon,klev))
    1309       ALLOCATE(qradice_lincont(klon,klev), qradice_circont(klon,klev))
    1310       ALLOCATE(dcfl_cir(klon,klev), dqtl_cir(klon,klev))
    1311       ALLOCATE(dcfl_ini(klon,klev), dqil_ini(klon,klev), dqtl_ini(klon,klev))
    1312       ALLOCATE(dcfl_sub(klon,klev), dqil_sub(klon,klev), dqtl_sub(klon,klev))
    1313       ALLOCATE(dcfl_mix(klon,klev), dqil_mix(klon,klev), dqtl_mix(klon,klev))
    1314       ALLOCATE(dcfl_sed(klon,klev), dqil_sed(klon,klev), dqtl_sed(klon,klev))
    1315       ALLOCATE(dcfl_auto(klon,klev), dqil_auto(klon,klev), dqtl_auto(klon,klev))
    1316       ALLOCATE(dcfc_sub(klon,klev), dqic_sub(klon,klev), dqtc_sub(klon,klev))
    1317       ALLOCATE(dcfc_mix(klon,klev), dqic_mix(klon,klev), dqtc_mix(klon,klev))
    1318       ALLOCATE(dcfc_sed(klon,klev), dqic_sed(klon,klev), dqtc_sed(klon,klev))
    1319       ALLOCATE(dcfc_auto(klon,klev), dqic_auto(klon,klev), dqtc_auto(klon,klev))
     1301      ALLOCATE(AEI_contrails(klon,klev), AEI_surv_contrails(klon,klev))
     1302      ALLOCATE(fsurv_contrails(klon,klev), section_contrails(klon,klev))
     1303      ALLOCATE(qice_cont(klon,klev))
     1304      ALLOCATE(qradice_cont(klon,klev))
     1305      ALLOCATE(dcfc_ini(klon,klev), dqic_ini(klon,klev), dqtc_ini(klon,klev), dnic_ini(klon,klev))
     1306      ALLOCATE(dcfc_sub(klon,klev), dqic_sub(klon,klev), dqtc_sub(klon,klev), dnic_sub(klon,klev))
     1307      ALLOCATE(dcfc_mix(klon,klev), dqic_mix(klon,klev), dqtc_mix(klon,klev), dnic_mix(klon,klev))
     1308      ALLOCATE(dnic_agg(klon,klev))
     1309      ALLOCATE(dcfc_sed(klon,klev), dqic_sed(klon,klev), dqtc_sed(klon,klev), dnic_sed(klon,klev))
     1310      ALLOCATE(dcfc_auto(klon,klev), dqic_auto(klon,klev), dqtc_auto(klon,klev), dnic_auto(klon,klev))
    13201311      ALLOCATE(cldfra_nocont(klon,klev), cldtau_nocont(klon,klev), cldemi_nocont(klon,klev))
    13211312      ALLOCATE(cldh_nocont(klon), contcov(klon), conttau(klon,klev), contemi(klon,klev))
    … …  
    17411732
    17421733!-- LSCP - aviation and contrails variables
    1743       DEALLOCATE(cfl_seri, d_cfl_dyn, cfc_seri, d_cfc_dyn)
    1744       DEALLOCATE(qtl_seri, d_qtl_dyn, qtc_seri, d_qtc_dyn)
    1745       DEALLOCATE(d_q_avi, flight_dist, flight_h2o)
     1734      DEALLOCATE(cfc_seri, d_cfc_dyn, qtc_seri, d_qtc_dyn, nic_seri, d_nic_dyn)
     1735      DEALLOCATE(d_q_avi, flight_dist, flight_fuel)
    17461736      DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
    1747       DEALLOCATE(qice_lincont, qice_circont, qradice_lincont, qradice_circont)
    1748       DEALLOCATE(dcfl_cir, dqtl_cir, dcfl_ini, dqil_ini)
    1749       DEALLOCATE(dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, dcfl_mix, dqil_mix, dqtl_mix)
    1750       DEALLOCATE(dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix)
    1751       DEALLOCATE(dcfl_sed, dqil_sed, dqtl_sed, dcfc_sed, dqic_sed, dqtc_sed)
    1752       DEALLOCATE(dcfl_auto, dqil_auto, dqtl_auto, dcfc_auto, dqic_auto, dqtc_auto)
     1737      DEALLOCATE(AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails)
     1738      DEALLOCATE(qice_cont, qradice_cont)
     1739      DEALLOCATE(dcfc_ini, dqic_ini, dqtc_ini, dnic_ini)
     1740      DEALLOCATE(dcfc_sub, dqic_sub, dqtc_sub, dnic_sub)
     1741      DEALLOCATE(dcfc_mix, dqic_mix, dqtc_mix, dnic_mix)
     1742      DEALLOCATE(dnic_agg)
     1743      DEALLOCATE(dcfc_sed, dqic_sed, dqtc_sed, dnic_sed)
     1744      DEALLOCATE(dcfc_auto, dqic_auto, dqtc_auto, dnic_auto)
    17531745      DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont, conttau, contemi)
    17541746      DEALLOCATE(cldh_nocont, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)

  • LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90

    r5779 r5790  
    25152515  TYPE(ctrl_out), SAVE :: o_dqavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    25162516    'dqavi', 'Water vapor emissions from aviation tendency', 'kg/kg/s', (/ ('',i=1,10) /))
    2517   TYPE(ctrl_out), SAVE :: o_cflseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2518     'cflseri', 'Linear contrail fraction', '-', (/ ('',i=1,10) /))
    2519   TYPE(ctrl_out), SAVE :: o_dcfldyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2520     'dcfldyn', 'Dynamics linear contrail fraction tendency', 's-1', (/ ('',i=1,10) /))
    25212517  TYPE(ctrl_out), SAVE :: o_cfcseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2522     'cfcseri', 'Contrail cirrus fraction', '-', (/ ('',i=1,10) /))
     2518    'cfcseri', 'Contrail fraction', '-', (/ ('',i=1,10) /))
    25232519  TYPE(ctrl_out), SAVE :: o_dcfcdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2524     'dcfcdyn', 'Dynamics contrail cirrus fraction tendency', 's-1', (/ ('',i=1,10) /))
    2525   TYPE(ctrl_out), SAVE :: o_qtlseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2526     'qtlseri', 'Linear contrail total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
    2527   TYPE(ctrl_out), SAVE :: o_dqtldyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2528     'dqtldyn', 'Dynamics linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
     2520    'dcfcdyn', 'Dynamics contrail fraction tendency', 's-1', (/ ('',i=1,10) /))
    25292521  TYPE(ctrl_out), SAVE :: o_qtcseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2530     'qtcseri', 'Contrail cirrus total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
     2522    'qtcseri', 'Contrail total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
    25312523  TYPE(ctrl_out), SAVE :: o_dqtcdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2532     'dqtcdyn', 'Dynamics contrail cirrus total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
     2524    'dqtcdyn', 'Dynamics contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
     2525  TYPE(ctrl_out), SAVE :: o_nicseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2526    'nicseri', 'Contrail ice crystals number concentration', '#/kg', (/ ('',i=1,10) /))
     2527  TYPE(ctrl_out), SAVE :: o_dnicdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2528    'dnicdyn', 'Dynamics contrail ice crystals number concentration tendency', '#/kg/s', (/ ('',i=1,10) /))
    25332529  TYPE(ctrl_out), SAVE :: o_Tcritcont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    25342530    'Tcritcont', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
    … …  
    25392535  TYPE(ctrl_out), SAVE :: o_potcontfraNP = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    25402536    'potcontfraNP', 'Potential non-persistent contrail fraction', '-', (/ ('', i=1,10)/))
    2541   TYPE(ctrl_out), SAVE :: o_qice_lincont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2542     'qice_lincont', 'Linear contrails ice specific humidity', 'kg/kg', (/ ('', i=1, 10) /))
    2543   TYPE(ctrl_out), SAVE :: o_qice_circont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2544     'qice_circont', 'Contrail cirrus ice specific humidity', 'kg/kg', (/ ('', i=1, 10) /))
    2545   TYPE(ctrl_out), SAVE :: o_dcflini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2546     'dcflini', 'Initial formation linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    2547   TYPE(ctrl_out), SAVE :: o_dqilini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2548     'dqilini', 'Initial formation linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2549   TYPE(ctrl_out), SAVE :: o_dqtlini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2550     'dqtlini', 'Initial formation linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2551   TYPE(ctrl_out), SAVE :: o_dcflcir = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2552     'dcflcir', 'Conversion to cirrus linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    2553   TYPE(ctrl_out), SAVE :: o_dqtlcir = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2554     'dqtlcir', 'Conversion to cirrus linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2555   TYPE(ctrl_out), SAVE :: o_dcflsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2556     'dcflsub', 'Sublimation linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    2557   TYPE(ctrl_out), SAVE :: o_dqilsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2558     'dqilsub', 'Sublimation linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2559   TYPE(ctrl_out), SAVE :: o_dqtlsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2560     'dqtlsub', 'Sublimation linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2561   TYPE(ctrl_out), SAVE :: o_dcflmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2562     'dcflmix', 'Mixing linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    2563   TYPE(ctrl_out), SAVE :: o_dqilmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2564     'dqilmix', 'Mixing linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2565   TYPE(ctrl_out), SAVE :: o_dqtlmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2566     'dqtlmix', 'Mixing linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2567   TYPE(ctrl_out), SAVE :: o_dcflsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2568     'dcflsed', 'Ice sedimentation linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    2569   TYPE(ctrl_out), SAVE :: o_dqilsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2570     'dqilsed', 'Ice sedimentation linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2571   TYPE(ctrl_out), SAVE :: o_dqtlsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2572     'dqtlsed', 'Ice sedimentation linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2537  TYPE(ctrl_out), SAVE :: o_qice_cont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2538    'qice_cont', 'Contrails ice specific humidity', 'kg/kg', (/ ('', i=1, 10) /))
     2539  TYPE(ctrl_out), SAVE :: o_dcfcini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2540    'dcfcini', 'Initial formation contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2541  TYPE(ctrl_out), SAVE :: o_dqicini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2542    'dqicini', 'Initial formation contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2543  TYPE(ctrl_out), SAVE :: o_dqtcini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2544    'dqtcini', 'Initial formation contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2545  TYPE(ctrl_out), SAVE :: o_dnicini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2546    'dnicini', 'Initial formation contrail ice crystals concentration tendency', '#/kg/s', (/ ('', i=1, 10) /))
    25732547  TYPE(ctrl_out), SAVE :: o_dcfcsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2574     'dcfcsed', 'Ice sedimentation contrail cirrus fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2548    'dcfcsed', 'Ice sedimentation contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    25752549  TYPE(ctrl_out), SAVE :: o_dqicsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2576     'dqicsed', 'Ice sedimentation contrail cirrus ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2550    'dqicsed', 'Ice sedimentation contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    25772551  TYPE(ctrl_out), SAVE :: o_dqtcsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2578     'dqtcsed', 'Ice sedimentation contrail cirrus total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2579   TYPE(ctrl_out), SAVE :: o_dcflauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2580     'dcflauto', 'Ice autoconversion linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    2581   TYPE(ctrl_out), SAVE :: o_dqilauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2582     'dqilauto', 'Ice autoconversion linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2583   TYPE(ctrl_out), SAVE :: o_dqtlauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2584     'dqtlauto', 'Ice autoconversion linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2552    'dqtcsed', 'Ice sedimentation contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2553  TYPE(ctrl_out), SAVE :: o_dnicsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2554    'dnicsed', 'Ice sedimentation contrail ice crystals concentration tendency', '#/kg/s', (/ ('', i=1, 10) /))
    25852555  TYPE(ctrl_out), SAVE :: o_dcfcauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2586     'dcfcauto', 'Ice autoconversion contrail cirrus fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2556    'dcfcauto', 'Ice autoconversion contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    25872557  TYPE(ctrl_out), SAVE :: o_dqicauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2588     'dqicauto', 'Ice autoconversion contrail cirrus ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2558    'dqicauto', 'Ice autoconversion contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    25892559  TYPE(ctrl_out), SAVE :: o_dqtcauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2590     'dqtcauto', 'Ice autoconversion contrail cirrus total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2560    'dqtcauto', 'Ice autoconversion contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2561  TYPE(ctrl_out), SAVE :: o_dnicauto = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2562    'dnicauto', 'Ice autoconversion contrail ice crystals concentration tendency', '#/kg/s', (/ ('', i=1, 10) /))
    25912563  TYPE(ctrl_out), SAVE :: o_dcfcsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2592     'dcfcsub', 'Sublimation contrail cirrus fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2564    'dcfcsub', 'Sublimation contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    25932565  TYPE(ctrl_out), SAVE :: o_dqicsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2594     'dqicsub', 'Sublimation contrail cirrus ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2566    'dqicsub', 'Sublimation contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    25952567  TYPE(ctrl_out), SAVE :: o_dqtcsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2596     'dqtcsub', 'Sublimation contrail cirrus total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2568    'dqtcsub', 'Sublimation contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2569  TYPE(ctrl_out), SAVE :: o_dnicsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2570    'dnicsub', 'Sublimation contrail ice crystals concentration tendency', '#/kg/s', (/ ('', i=1, 10) /))
    25972571  TYPE(ctrl_out), SAVE :: o_dcfcmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2598     'dcfcmix', 'Mixing contrail cirrus fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2572    'dcfcmix', 'Mixing contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
    25992573  TYPE(ctrl_out), SAVE :: o_dqicmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2600     'dqicmix', 'Mixing contrail cirrus ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2574    'dqicmix', 'Mixing contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    26012575  TYPE(ctrl_out), SAVE :: o_dqtcmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2602     'dqtcmix', 'Mixing contrail cirrus total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2576    'dqtcmix', 'Mixing contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2577  TYPE(ctrl_out), SAVE :: o_dnicmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2578    'dnicmix', 'Mixing contrail ice crystals concentration tendency', '#/kg/s', (/ ('', i=1, 10) /))
     2579  TYPE(ctrl_out), SAVE :: o_dnicagg = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2580    'dnicagg', 'Aggregation contrail ice crystals concentration tendency', '#/kg/s', (/ ('', i=1, 10) /))
    26032581  TYPE(ctrl_out), SAVE :: o_flight_dist = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    26042582    'flightdist', 'Aviation flown distance', 'm/s/m^3', (/ ('', i=1, 10) /))
    2605   TYPE(ctrl_out), SAVE :: o_flight_h2o = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2606     'flighth2o', 'Aviation H2O flight emission', 'kg H2O/s/m^3', (/ ('', i=1, 10) /))
     2583  TYPE(ctrl_out), SAVE :: o_flight_fuel = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2584    'flightfuel', 'Aviation fuel consumption', 'kg/s/m^3', (/ ('', i=1, 10) /))
     2585  TYPE(ctrl_out), SAVE :: o_AEI_contrails = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2586    'AEI_cont', 'Apparent emission index contrails', '#/kg', (/ ('', i=1, 10) /))
     2587  TYPE(ctrl_out), SAVE :: o_AEI_surv_contrails = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2588    'AEI_surv_cont', 'Apparent emission index contrails after vortex loss', '#/kg', (/ ('', i=1, 10) /))
     2589  TYPE(ctrl_out), SAVE :: o_fsurv_contrails = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2590    'fsurv_cont', 'Survival fraction after vortex loss', '-', (/ ('', i=1, 10) /))
     2591  TYPE(ctrl_out), SAVE :: o_section_contrails = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2592    'section_cont', 'Cross section of newly formed contrails', 'm2', (/ ('', i=1, 10) /))
    26072593  TYPE(ctrl_out), SAVE :: o_cldfra_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    26082594    'cldfra_nocont', 'Cloud fraction w/o contrails', '-', (/ ('', i=1, 10) /))

  • LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90

    r5779 r5790  
    239239         o_issrfra250to300, o_issrfra300to400, o_issrfra400to500, &
    240240!-- LSCP - aviation variables
    241          o_cflseri, o_dcfldyn, o_cfcseri, o_dcfcdyn, &
    242          o_qtlseri, o_dqtldyn, o_qtcseri, o_dqtcdyn, o_dqavi, &
     241         o_cfcseri, o_dcfcdyn, o_qtcseri, o_dqtcdyn, o_nicseri, o_dnicdyn, o_dqavi, &
    243242         o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
    244          o_flight_dist, o_flight_h2o, o_qice_lincont, o_qice_circont, o_dcflcir, o_dqtlcir, &
    245          o_dcflini, o_dqilini, o_dqtlini, o_dcflsub, o_dqilsub, o_dqtlsub, &
    246          o_dcflmix, o_dqilmix, o_dqtlmix, &
    247          o_dcflsed, o_dqilsed, o_dqtlsed, o_dcfcsed, o_dqicsed, o_dqtcsed, &
    248          o_dcflauto, o_dqilauto, o_dqtlauto, o_dcfcauto, o_dqicauto, o_dqtcauto, &
    249          o_dcfcsub, o_dqicsub, o_dqtcsub, o_dcfcmix, o_dqicmix, o_dqtcmix, &
     243         o_AEI_contrails, o_AEI_surv_contrails, o_fsurv_contrails, o_section_contrails, &
     244         o_flight_dist, o_flight_fuel, o_qice_cont, &
     245         o_dcfcini, o_dqicini, o_dqtcini, o_dnicini, &
     246         o_dcfcmix, o_dqicmix, o_dqtcmix, o_dnicmix, o_dnicagg, &
     247         o_dcfcsub, o_dqicsub, o_dqtcsub, o_dnicsub, &
     248         o_dcfcsed, o_dqicsed, o_dqtcsed, o_dnicsed, &
     249         o_dcfcauto, o_dqicauto, o_dqtcauto, o_dnicauto, &
    250250         o_cldfra_nocont, o_cldtau_nocont, o_cldemi_nocont, o_cldh_nocont, &
    251251         o_contcov, o_conttau, o_contemi, o_iwp_nocont, o_iwc_nocont, o_ref_ice_nocont, &
    … …  
    379379         issrfra100to150, issrfra150to200, issrfra200to250, &
    380380         issrfra250to300, issrfra300to400, issrfra400to500, &
    381          cfl_seri, d_cfl_dyn, cfc_seri, d_cfc_dyn, &
    382          qtl_seri, d_qtl_dyn, qtc_seri, d_qtc_dyn, d_q_avi, &
     381         cfc_seri, d_cfc_dyn, qtc_seri, d_qtc_dyn, nic_seri, d_nic_dyn, d_q_avi, &
    383382         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    384          flight_dist, flight_h2o, qice_lincont, qice_circont, dcfl_cir, dqtl_cir, &
    385          dcfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, &
    386          dcfl_mix, dqil_mix, dqtl_mix, &
    387          dcfl_sed, dqil_sed, dqtl_sed, dcfc_sed, dqic_sed, dqtc_sed, &
    388          dcfl_auto, dqil_auto, dqtl_auto, dcfc_auto, dqic_auto, dqtc_auto, &
    389          dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix, &
     383         AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, &
     384         flight_dist, flight_fuel, qice_cont, &
     385         dcfc_ini, dqic_ini, dqtc_ini, dnic_ini, &
     386         dcfc_sub, dqic_sub, dqtc_sub, dnic_sub, &
     387         dcfc_mix, dqic_mix, dqtc_mix, dnic_mix, dnic_agg, &
     388         dcfc_sed, dqic_sed, dqtc_sed, dnic_sed, &
     389         dcfc_auto, dqic_auto, dqtc_auto, dnic_auto, &
    390390         cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
    391391         contcov, conttau, contemi, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
    … …  
    22372237!-- LSCP - aviation variables
    22382238       IF (ok_plane_contrail) THEN
    2239          CALL histwrite_phy(o_cflseri, cfl_seri)
    2240          CALL histwrite_phy(o_dcfldyn, d_cfl_dyn)
    22412239         CALL histwrite_phy(o_cfcseri, cfc_seri)
    22422240         CALL histwrite_phy(o_dcfcdyn, d_cfc_dyn)
    2243          CALL histwrite_phy(o_qtlseri, qtl_seri)
    2244          CALL histwrite_phy(o_dqtldyn, d_qtl_dyn)
    22452241         CALL histwrite_phy(o_qtcseri, qtc_seri)
    22462242         CALL histwrite_phy(o_dqtcdyn, d_qtc_dyn)
     2243         CALL histwrite_phy(o_nicseri, nic_seri)
     2244         CALL histwrite_phy(o_dnicdyn, d_nic_dyn)
    22472245         CALL histwrite_phy(o_flight_dist, flight_dist)
     2246         CALL histwrite_phy(o_flight_fuel, flight_fuel)
    22482247         CALL histwrite_phy(o_Tcritcont, Tcritcont)
    22492248         CALL histwrite_phy(o_qcritcont, qcritcont)
    22502249         CALL histwrite_phy(o_potcontfraP, potcontfraP)
    22512250         CALL histwrite_phy(o_potcontfraNP, potcontfraNP)
    2252          CALL histwrite_phy(o_qice_lincont, qice_lincont)
    2253          CALL histwrite_phy(o_qice_circont, qice_circont)
    2254          CALL histwrite_phy(o_dcflcir, dcfl_cir)
    2255          CALL histwrite_phy(o_dqtlcir, dqtl_cir)
    2256          CALL histwrite_phy(o_dcflini, dcfl_ini)
    2257          CALL histwrite_phy(o_dqilini, dqil_ini)
    2258          CALL histwrite_phy(o_dqtlini, dqtl_ini)
    2259          CALL histwrite_phy(o_dcflsub, dcfl_sub)
    2260          CALL histwrite_phy(o_dqilsub, dqil_sub)
    2261          CALL histwrite_phy(o_dqtlsub, dqtl_sub)
    2262          CALL histwrite_phy(o_dcflmix, dcfl_mix)
    2263          CALL histwrite_phy(o_dqilmix, dqil_mix)
    2264          CALL histwrite_phy(o_dqtlmix, dqtl_mix)
    2265          CALL histwrite_phy(o_dcflsed, dcfl_sed)
    2266          CALL histwrite_phy(o_dqilsed, dqil_sed)
    2267          CALL histwrite_phy(o_dqtlsed, dqtl_sed)
    2268          CALL histwrite_phy(o_dcflauto, dcfl_auto)
    2269          CALL histwrite_phy(o_dqilauto, dqil_auto)
    2270          CALL histwrite_phy(o_dqtlauto, dqtl_auto)
     2251         CALL histwrite_phy(o_AEI_contrails, AEI_contrails)
     2252         CALL histwrite_phy(o_AEI_surv_contrails, AEI_surv_contrails)
     2253         CALL histwrite_phy(o_fsurv_contrails, fsurv_contrails)
     2254         CALL histwrite_phy(o_section_contrails, section_contrails)
     2255         CALL histwrite_phy(o_qice_cont, qice_cont)
     2256         CALL histwrite_phy(o_dcfcini, dcfc_ini)
     2257         CALL histwrite_phy(o_dqicini, dqic_ini)
     2258         CALL histwrite_phy(o_dqtcini, dqtc_ini)
     2259         CALL histwrite_phy(o_dnicini, dnic_ini)
    22712260         CALL histwrite_phy(o_dcfcsub, dcfc_sub)
    22722261         CALL histwrite_phy(o_dqicsub, dqic_sub)
    22732262         CALL histwrite_phy(o_dqtcsub, dqtc_sub)
     2263         CALL histwrite_phy(o_dnicsub, dnic_sub)
    22742264         CALL histwrite_phy(o_dcfcmix, dcfc_mix)
    22752265         CALL histwrite_phy(o_dqicmix, dqic_mix)
    22762266         CALL histwrite_phy(o_dqtcmix, dqtc_mix)
     2267         CALL histwrite_phy(o_dnicmix, dnic_mix)
     2268         CALL histwrite_phy(o_dnicagg, dnic_agg)
    22772269         CALL histwrite_phy(o_dcfcsed, dcfc_sed)
    22782270         CALL histwrite_phy(o_dqicsed, dqic_sed)
    22792271         CALL histwrite_phy(o_dqtcsed, dqtc_sed)
     2272         CALL histwrite_phy(o_dnicsed, dnic_sed)
    22802273         CALL histwrite_phy(o_dcfcauto, dcfc_auto)
    22812274         CALL histwrite_phy(o_dqicauto, dqic_auto)
    22822275         CALL histwrite_phy(o_dqtcauto, dqtc_auto)
     2276         CALL histwrite_phy(o_dnicauto, dnic_auto)
    22832277         CALL histwrite_phy(o_cldfra_nocont, cldfra_nocont)
    22842278         CALL histwrite_phy(o_cldtau_nocont, cldtau_nocont)
    … …  
    23012295           CALL histwrite_phy(o_soll_nocont, sollw_nocont)
    23022296         ENDIF
    2303        ENDIF
    2304        IF (ok_plane_h2o) THEN
    2305          CALL histwrite_phy(o_flight_h2o, flight_h2o)
    2306          CALL histwrite_phy(o_dqavi, d_q_avi)
     2297         IF (ok_plane_h2o) THEN
     2298           CALL histwrite_phy(o_dqavi, d_q_avi)
     2299         ENDIF
    23072300       ENDIF
    23082301       

  • LMDZ6/branches/contrails/libf/phylmd/phys_state_var_mod.F90

    r5717 r5790  
    142142      REAL, ALLOCATABLE, SAVE :: qvcon(:,:), qccon(:,:)
    143143!$OMP THREADPRIVATE(qvcon, qccon)
    144       REAL, ALLOCATABLE, SAVE :: cfl_ancien(:,:), cfc_ancien(:,:)
    145 !$OMP THREADPRIVATE(cfl_ancien, cfc_ancien)
    146       REAL, ALLOCATABLE, SAVE :: qtl_ancien(:,:), qtc_ancien(:,:)
    147 !$OMP THREADPRIVATE(qtl_ancien, qtc_ancien)
     144      REAL, ALLOCATABLE, SAVE :: cfc_ancien(:,:), qtc_ancien(:,:), nic_ancien(:,:)
     145!$OMP THREADPRIVATE(cfc_ancien, qtc_ancien, nic_ancien)
    148146!!! RomP >>>
    149147      REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
    … …  
    662660      ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
    663661      ALLOCATE(cf_ancien(klon,klev), qvc_ancien(klon,klev))
    664       ALLOCATE(cfl_ancien(klon,klev), cfc_ancien(klon,klev))
    665       ALLOCATE(qtl_ancien(klon,klev), qtc_ancien(klon,klev))
     662      ALLOCATE(cfc_ancien(klon,klev), qtc_ancien(klon,klev), nic_ancien(klon,klev))
    666663      ALLOCATE(qvcon(klon,klev), qccon(klon,klev))
    667664!!! Rom P >>>
    … …  
    914911      DEALLOCATE(u_ancien, v_ancien)
    915912      DEALLOCATE(cf_ancien, qvc_ancien)
    916       DEALLOCATE(cfl_ancien, cfc_ancien, qtl_ancien, qtc_ancien)
     913      DEALLOCATE(cfc_ancien, qtc_ancien, nic_ancien)
    917914      DEALLOCATE(qvcon, qccon)
    918915      DEALLOCATE(tr_ancien)                           !RomP

  • LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90

    r5779 r5790  
    3939    USE ioipsl_getin_p_mod, ONLY : getin_p
    4040    USE indice_sol_mod
    41     USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, addPhase, ivap, iliq, isol, ibs, icf, iqvc, icfl, icfc, iqtl, iqtc
     41    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, tracers, type_trac, addPhase, ivap, iliq, isol, ibs, icf, iqvc, icfc, iqtc, inic
    4242    USE strings_mod,  ONLY: strIdx
    4343    USE iophy
    … …  
    337337       dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    338338       !-- LSCP - aviation and contrails variables
    339        cfl_seri, cfc_seri, qtl_seri, qtc_seri, d_cfl_dyn, d_cfc_dyn, d_qtl_dyn, d_qtc_dyn, &
    340        d_q_avi, flight_dist, flight_h2o, &
    341        qice_lincont, qice_circont, qradice_lincont, qradice_circont, &
     339       cfc_seri, qtc_seri, nic_seri, d_cfc_dyn, d_qtc_dyn, d_nic_dyn, &
     340       d_q_avi, flight_dist, flight_fuel, &
     341       qice_cont, qradice_cont, &
    342342       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    343343       cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
    … …  
    14271427       ENDIF
    14281428
    1429        IF (ok_plane_contrail.AND.((icfl.EQ.0).OR.(icfc.EQ.0).OR.(iqtl.EQ.0).OR.(iqtc.EQ.0))) THEN
     1429       IF (ok_plane_contrail.AND.((icfc.EQ.0).OR.(iqtc.EQ.0).OR.(inic.EQ.0))) THEN
    14301430          WRITE (lunout, *) ' ok_plane_contrail=y requires 9 tracers ', &
    1431               '(H2O_g, H2O_l, H2O_s, CLDFRA, CLDVAP, LINCONTFRA, ', &
    1432               'CIRCONTFRA, LINCONTWATER, CIRCONTWATER) ', &
     1431              '(H2O_g, H2O_l, H2O_s, CLDFRA, CLDVAP, CONTFRA, ', &
     1432              'CONTWATER, CONTNICE) ', &
    14331433              'but not all are provided. Might as well stop here.'
    14341434          abort_message='see above'
    … …  
    24572457          cf_seri(i,k) = 0.
    24582458          qvc_seri(i,k)= 0.
    2459           cfl_seri(i,k)= 0.
    24602459          cfc_seri(i,k)= 0.
    2461           qtl_seri(i,k)= 0.
    24622460          qtc_seri(i,k)= 0.
     2461          nic_seri(i,k)= 0.
    24632462          !CR: ATTENTION, on rajoute la variable glace
    24642463          IF (nqo .GE. 3) THEN        !--vapour, liquid and ice
    … …  
    24752474          ENDIF
    24762475          IF (ok_plane_contrail) THEN
    2477              cfl_seri(i,k) = qx(i,k,icfl)
    24782476             cfc_seri(i,k) = qx(i,k,icfc)
    2479              qtl_seri(i,k) = qx(i,k,iqtl)
    24802477             qtc_seri(i,k) = qx(i,k,iqtc)
     2478             nic_seri(i,k) = qx(i,k,inic)
    24812479             !--DYNAMICO can return NaNs for children tracers
    2482              IF (ISNAN(cfl_seri(i,k))) cfl_seri(i,k) = 0.
    24832480             IF (ISNAN(cfc_seri(i,k))) cfc_seri(i,k) = 0.
    2484              IF (ISNAN(qtl_seri(i,k))) qtl_seri(i,k) = 0.
    24852481             IF (ISNAN(qtc_seri(i,k))) qtc_seri(i,k) = 0.
     2482             IF (ISNAN(nic_seri(i,k))) nic_seri(i,k) = 0.
    24862483          ENDIF
    24872484       ENDDO
    … …  
    25582555       d_cf_dyn(:,:) = (cf_seri(:,:)-cf_ancien(:,:))/phys_tstep
    25592556       d_qvc_dyn(:,:)= (qvc_seri(:,:)-qvc_ancien(:,:))/phys_tstep
    2560        d_cfl_dyn(:,:)= (cfl_seri(:,:)-cfl_ancien(:,:))/phys_tstep
    25612557       d_cfc_dyn(:,:)= (cfc_seri(:,:)-cfc_ancien(:,:))/phys_tstep
    2562        d_qtl_dyn(:,:)= (qtl_seri(:,:)-qtl_ancien(:,:))/phys_tstep
    25632558       d_qtc_dyn(:,:)= (qtc_seri(:,:)-qtc_ancien(:,:))/phys_tstep
     2559       d_nic_dyn(:,:)= (nic_seri(:,:)-nic_ancien(:,:))/phys_tstep
    25642560       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
    25652561       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
    … …  
    25832579       d_cf_dyn(:,:) = 0.0
    25842580       d_qvc_dyn(:,:)= 0.0
    2585        d_cfl_dyn(:,:)= 0.0
    25862581       d_cfc_dyn(:,:)= 0.0
    2587        d_qtl_dyn(:,:)= 0.0
    25882582       d_qtc_dyn(:,:)= 0.0
     2583       d_nic_dyn(:,:)= 0.0
    25892584       d_q_dyn2d(:)  = 0.0
    25902585       d_ql_dyn2d(:) = 0.0
    … …  
    27112706            qi_seri_lscp(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27122707            qvc_seri(i,k) = qvc_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    2713             qtl_seri(i,k) = qtl_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27142708            qtc_seri(i,k) = qtc_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27152709          ELSE
    … …  
    27172711            qi_seri_lscp(i,k) = 0.
    27182712            qvc_seri(i,k) = 0.
    2719             qtl_seri(i,k) = 0.
    27202713            qtc_seri(i,k) = 0.
    27212714          ENDIF
    … …  
    39463939         ratqs,ratqsc,ratqs_inter_,sigma_qtherm)
    39473940
     3941    ! Vertical interpolation is done at each physical timestep
     3942    IF (ok_plane_contrail) CALL vertical_interpolation_aviation( &
     3943            klon, klev, paprs, pplay, t_seri, flight_dist, flight_fuel)
     3944
    39483945    !--Add the water emissions from aviation
    39493946    IF ( ok_plane_h2o ) THEN
    39503947       CALL aviation_water_emissions(klon, klev, phys_tstep, pplay, &
    3951             t_seri, q_seri, flight_h2o, d_q_avi)
     3948            t_seri, q_seri, flight_fuel, d_q_avi)
    39523949       CALL add_phys_tend(du0, dv0, dt0, d_q_avi, dql0, dqi0, dqbs0, paprs, &
    39533950            'avi', abortphy, flag_inhib_tend, itap, 0)
    … …  
    39693966
    39703967    IF (ok_no_issr_strato) CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
    3971     ! Vertical interpolation is done at each physical timestep
    3972     IF (ok_plane_contrail) CALL vertical_interpolation_aviation(klon, klev, paprs, pplay, t_seri, flight_dist, flight_h2o)
    39733968
    39743969    IF (ok_ice_supersat) THEN
    … …  
    39823977      DO k = 1, klev
    39833978        DO i = 1, klon
    3984           qtl_seri(i,k) = qtl_seri(i,k) * q_seri(i,k)
    39853979          qtc_seri(i,k) = qtc_seri(i,k) * q_seri(i,k)
    39863980        ENDDO
    … …  
    40144008         dqi_adj, dqi_sub, dqi_con, dqi_mix, &
    40154009         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    4016          cfl_seri, cfc_seri, qtl_seri, qtc_seri, qice_lincont, qice_circont, &
    4017          flight_dist, flight_h2o, qradice_lincont, qradice_circont, &
     4010         cfc_seri, qtc_seri, nic_seri, qice_cont, &
     4011         flight_dist, flight_fuel, qradice_cont, &
    40184012         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    40194013         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
    … …  
    45654559               zfice, dNovrN, ptconv, rnebcon, clwcon, &
    45664560               !--AB contrails
    4567                cfl_seri, cfc_seri, qradice_lincont, qradice_circont, cldfra_nocont, &
     4561               cfc_seri, qradice_cont, nic_seri, cldfra_nocont, &
    45684562               cldtau_nocont, cldemi_nocont, conttau, contemi, cldh_nocont, contcov, &
    45694563               fiwp_nocont, fiwc_nocont, ref_ice_nocont)
    … …  
    57745768          ENDIF
    57755769          IF (ok_plane_contrail) THEN
    5776              d_qx(i,k,icfl) = ( cfl_seri(i,k) - qx(i,k,icfl) ) / phys_tstep
    57775770             d_qx(i,k,icfc) = ( cfc_seri(i,k) - qx(i,k,icfc) ) / phys_tstep
    5778              d_qx(i,k,iqtl) = ( qtl_seri(i,k) - qx(i,k,iqtl) ) / phys_tstep
    57795771             d_qx(i,k,iqtc) = ( qtc_seri(i,k) - qx(i,k,iqtc) ) / phys_tstep
     5772             d_qx(i,k,inic) = ( nic_seri(i,k) - qx(i,k,inic) ) / phys_tstep
    57805773          ENDIF
    57815774       ENDDO
    … …  
    58095802    cf_ancien(:,:) = cf_seri(:,:)
    58105803    qvc_ancien(:,:)= qvc_seri(:,:)
    5811     cfl_ancien(:,:)= cfl_seri(:,:)
    58125804    cfc_ancien(:,:)= cfc_seri(:,:)
    5813     qtl_ancien(:,:)= qtl_seri(:,:)
    58145805    qtc_ancien(:,:)= qtc_seri(:,:)
     5806    nic_ancien(:,:)= nic_seri(:,:)
    58155807    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
    58165808    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
Note: See TracChangeset for help on using the changeset viewer.