Changeset 5797 for LMDZ6


Ignore:
Timestamp:
Aug 5, 2025, 2:22:14 PM (8 days ago)
Author:
aborella
Message:

Bugfix for saturation adjustment in cirrus mixing + bugfix for contrails sedimentation + new diagnostics + support for unadjusted contrails

Location:
LMDZ6/branches/contrails
Files:
15 edited

Legend:

Unmodified
Added
Removed
  • LMDZ6/branches/contrails/DefLists/field_def_lmdz.xml

    r5796 r5797  
    993993        <field id="qtcseri"    long_name="Contrail total specific humidity" unit="kg/kg" />
    994994        <field id="dqtcdyn"    long_name="Dynamics contrail total specific humidity tendency" unit="kg/kg/s" />
     995        <field id="qicseri"    long_name="Contrail ice specific humidity" unit="kg/kg" />
     996        <field id="dqicdyn"    long_name="Dynamics contrail ice specific humidity tendency" unit="kg/kg/s" />
    995997        <field id="nicseri"    long_name="Contrail ice crystals concentration" unit="#/kg" />
    996998        <field id="dnicdyn"    long_name="Dynamics contrail ice crystals concentration tendency" unit="#/kg/s" />
     
    10131015        <field id="dnicmix"    long_name="Mixing contrail ice crystals concentration tendency"    unit="#/kg/s" />
    10141016        <field id="dnicagg"    long_name="Aggregation contrail ice crystals concentration tendency"    unit="#/kg/s" />
     1017        <field id="dqicadj"    long_name="Phase adjustment contrail ice specific humidity tendency"    unit="kg/kg/s" />
     1018        <field id="dqtcadj"    long_name="Phase adjustment contrail total specific humidity tendency"    unit="kg/kg/s" />
    10151019        <field id="dcfcsed"    long_name="Ice sedimentation contrail fraction tendency"    unit="s-1" />
    10161020        <field id="dqicsed"    long_name="Ice sedimentation contrail ice specific humidity tendency"    unit="kg/kg/s" />
     
    10311035        <field id="rvol_ygcont"   long_name="Ice crystals volumic radius in young contrails"    unit="microns" detect_missing_value=".true." />
    10321036        <field id="tau_ygcont"    long_name="Young contrails optical depth"    unit="-" detect_missing_value=".true." />
     1037        <field id="nice_vscont"   long_name="Ice particle number concentration in visible contrails"    unit="#/cm3" detect_missing_value=".true." />
     1038        <field id="iwc_vscont"    long_name="Ice water content in visible contrails"    unit="g/m3" detect_missing_value=".true." />
     1039        <field id="rvol_vscont"   long_name="Ice crystals volumic radius in visible contrails"    unit="microns" detect_missing_value=".true." />
     1040        <field id="tau_vscont"    long_name="Visible contrails optical depth"    unit="-" detect_missing_value=".true." />
    10331041        <field id="nice_cont"     long_name="Ice particle number concentration in contrails"    unit="#/cm3" detect_missing_value=".true." />
    10341042        <field id="iwc_cont"      long_name="Ice water content in contrails"    unit="g/m3" detect_missing_value=".true." />
  • LMDZ6/branches/contrails/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90

    r5791 r5797  
    5050    ale_bl, ale_bl_trig, alp_bl, &
    5151    ale_wake, ale_bl_stat, AWAKE_S, &
    52     cf_ancien, qvc_ancien, qvcon, qccon, cfc_ancien, qtc_ancien, nic_ancien
     52    cf_ancien, qvc_ancien, qvcon, qccon, cfc_ancien, qtc_ancien, qic_ancien, nic_ancien
    5353
    5454  USE comconst_mod, ONLY: pi, dtvr
     
    251251  cfc_ancien = 0.
    252252  qtc_ancien = 0.
     253  qic_ancien = 0.
    253254  nic_ancien = 0.
    254255 
  • LMDZ6/branches/contrails/libf/phylmd/infotrac_phy.F90

    r5791 r5797  
    2727   !=== FOR ISOTOPES: Specific to water
    2828   PUBLIC :: iH2O                                          !--- Value of "ixIso" for "H2O" isotopes class
    29    PUBLIC :: ivap, iliq, isol, ibs, icf, iqvc, itke, icfc, iqtc, inic
     29   PUBLIC :: ivap, iliq, isol, ibs, icf, iqvc, itke, icfc, iqtc, iqic, 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, itke, icfc, iqtc, inic
    106 !$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, itke, icfc, iqtc, inic)
     105   INTEGER, SAVE :: ivap, iliq, isol, ibs, icf, iqvc, itke, icfc, iqtc, iqic, inic
     106!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, itke, icfc, iqtc, iqic, inic)
    107107
    108108   !=== DIMENSIONS OF THE TRACERS TABLES, TRACERS TYPE(S)
     
    384384   icfc = strIdx(tracers(:)%name, 'CONTFRA')
    385385   iqtc = strIdx(tracers(:)%name, 'CONTWATER')
     386   iqic = strIdx(tracers(:)%name, 'CONTICE')
    386387   inic = strIdx(tracers(:)%name, 'CONTNICE')
    387388   !--Two ways of declaring tracers - the way below should be deleted in the future
     
    390391   IF (icfc.EQ.0) icfc = strIdx(tracers(:)%name, addPhase('H2O', 'x'))
    391392   IF (iqtc.EQ.0) iqtc = strIdx(tracers(:)%name, addPhase('H2O', 'y'))
     393   IF (iqic.EQ.0) iqic = strIdx(tracers(:)%name, addPhase('H2O', 'w'))
    392394   IF (inic.EQ.0) inic = strIdx(tracers(:)%name, addPhase('H2O', 'z'))
    393395   !--- Compute indices for TKE when it is advected
  • LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90

    r5796 r5797  
    238238        dist_contrails = flight_dist(i) * dtime * potcontfraP(i)
    239239        dcfc_ini(i) = dist_contrails * section_contrails(i)
    240         dqtc_ini(i) = icesat_ratio * qsat(i) * dcfc_ini(i)
    241         dqic_ini(i) = dqtc_ini(i) - qsat(i) * dcfc_ini(i)
     240        dqtc_ini(i) = qsat(i) * dcfc_ini(i)
     241        dqic_ini(i) = ( icesat_ratio - 1. ) * qsat(i) * dcfc_ini(i)
    242242        dnic_ini(i) = Nice_per_m_init_contrails * dist_contrails / rho(i)
    243243      ENDIF
  • LMDZ6/branches/contrails/libf/phylmd/lmdz_cloud_optics_prop.f90

    r5796 r5797  
    730730          rei_cont = (mice_cont / rho_ice * 3. / 4. / RPI)**(1./3.)
    731731          !--Effective radius (in microns)
    732           rei_cont = MIN(100., MAX(rei_min_contrails, rei_cont / eff2vol_radius_contrails * 1e6))
     732          rei_cont = MIN(200., MAX(rei_min_contrails, rei_cont / eff2vol_radius_contrails * 1e6))
    733733          zfiwp_var = 1000.*xfiwc_cont(i, k)/pclc_cont(i, k)*rhodz(i, k)
    734734
  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90

    r5796 r5797  
    102102      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, dqi_auto, &
    103103      dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, dqvc_sed, dqvc_auto, &
    104       contfra_in, qtc_in, nic_in, flight_dist, flight_fuel, &
    105       contfra, qcont, Ncont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
     104      contfra_in, qvcont_in, qicont_in, nic_in, flight_dist, flight_fuel, &
     105      contfra, qcont, qvcont, Ncont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    106106      AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, &
    107107      dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv, &
    108108      dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont, &
    109109      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, &
     110      dcfc_mix, dqic_mix, dqtc_mix, dnic_mix, dnic_agg, dqic_adj, dqtc_adj, &
    111111      dcfc_sed, dqic_sed, dqtc_sed, dnic_sed, dcfc_auto, dqic_auto, dqtc_auto, dnic_auto)
    112112
     
    137137USE lmdz_lscp_ini,   ONLY: chi_sedim
    138138
    139 USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_contrails
    140 USE lmdz_lscp_ini,   ONLY: coef_mixing_contrails, coef_shear_contrails
     139USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, ok_unadjusted_contrails
     140USE lmdz_lscp_ini,   ONLY: aspect_ratio_contrails, coef_mixing_contrails, coef_shear_contrails
    141141USE lmdz_lscp_ini,   ONLY: chi_mixing_contrails, eff2vol_radius_contrails
    142142USE lmdz_lscp_tools, ONLY: ICECRYSTAL_VELO
     
    180180!
    181181REAL,     INTENT(IN)   , DIMENSION(klon) :: contfra_in    ! input contrails fraction [-]
    182 REAL,     INTENT(IN)   , DIMENSION(klon) :: qtc_in        ! input contrails total specific humidity [kg/kg]
     182REAL,     INTENT(IN)   , DIMENSION(klon) :: qvcont_in     ! input contrails specific humidity [kg/kg]
     183REAL,     INTENT(IN)   , DIMENSION(klon) :: qicont_in     ! input contrails ice specific humidity [kg/kg]
    183184REAL,     INTENT(IN)   , DIMENSION(klon) :: nic_in        ! input contrails ice crystals concentration [#/kg]
    184185REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_dist   ! aviation distance flown concentration [m/s/m3]
     
    231232!
    232233REAL,     INTENT(INOUT), DIMENSION(klon) :: contfra      ! contrail fraction [-]
    233 REAL,     INTENT(INOUT), DIMENSION(klon) :: qcont        ! contrail specific humidity [kg/kg]
     234REAL,     INTENT(INOUT), DIMENSION(klon) :: qcont        ! contrail total specific humidity [kg/kg]
     235REAL,     INTENT(INOUT), DIMENSION(klon) :: qvcont       ! contrail specific humidity [kg/kg]
    234236REAL,     INTENT(INOUT), DIMENSION(klon) :: Ncont        ! contrail ice crystals concentration [#/kg]
    235237REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
     
    249251REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_sub     ! contrails total specific humidity tendency because of sublimation [kg/kg/s]
    250252REAL,     INTENT(INOUT), DIMENSION(klon) :: dnic_sub     ! contrails ice crystals concentration tendency because of sublimation [#/kg/s]
     253REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_adj     ! contrails ice specific humidity tendency because of phase adjustment [kg/kg/s]
     254REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_adj     ! contrails total specific humidity tendency because of phase adjustment [kg/kg/s]
    251255REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_mix     ! contrails cloud fraction tendency because of mixing [s-1]
    252256REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_mix     ! contrails ice specific humidity tendency because of mixing [kg/kg/s]
     
    319323!
    320324! for contrails
    321 REAL :: mice, Niceincld
     325REAL :: mice, Niceincld, qvapincont
    322326REAL, DIMENSION(klon) :: qised_cont_abv, Nised_cont_abv
    323327REAL :: dcfc_sed1, dcfc_sed2, dqtc_sed1, dnic_sed1, dqtc_sed2, dqic_sed1, dqic_sed2, dnic_sed2
     
    358362      contfra(i) = 0.
    359363      qcont(i) = 0.
     364      qvcont(i) = 0.
     365      Ncont(i) = 0.
    360366      IF ( ok_ice_sedim ) THEN
    361367        flsed_cont(i) = 0.
     
    534540        !--conserve order relations
    535541        contfra(i) = MAX(0., MIN(cldfra(i), contfra_in(i)))
    536         qcont(i) = MAX(0., MIN(qcld(i), qtc_in(i)))
    537542        Ncont(i) = MAX(0., nic_in(i))
     543        IF ( ok_unadjusted_contrails ) THEN
     544          qcont(i) = MAX(0., MIN(qcld(i), qvcont_in(i) + qicont_in(i)))
     545          qvcont(i) = MAX(0., MIN(qcont(i), qvcont_in(i)))
     546        ELSE
     547          qcont(i) = MAX(0., MIN(qcld(i), qvcont_in(i)))
     548          qvcont(i) = qsat(i) * contfra(i)
     549        ENDIF
    538550
    539551        qised_cont_abv(i) = 0.
     
    548560          qised_cont_abv(i) = MIN(qclr(i), flsed_cont_abv(i) &
    549561              / ( paprsdn(i) - paprsup(i) ) * RG * dtime)
    550           Nised_cont_abv(i) = MIN(qclr(i), Nflsed_cont_abv(i) &
    551               / ( paprsdn(i) - paprsup(i) ) * RG * dtime)
     562          Nised_cont_abv(i) = Nflsed_cont_abv(i) &
     563              / ( paprsdn(i) - paprsup(i) ) * RG * dtime
    552564          qclr(i) = qclr(i) - qised_cont_abv(i)
    553565        ENDIF
     
    561573        dqtc_sub(i) = 0.
    562574        dnic_sub(i) = 0.
     575        dqic_adj(i) = 0.
     576        dqtc_adj(i) = 0.
    563577        dcfc_mix(i) = 0.
    564578        dqic_mix(i) = 0.
     
    577591        contfra(i) = 0.
    578592        qcont(i) = 0.
     593        qvcont(i) = 0.
     594        Ncont(i) = 0.
    579595      ENDIF
    580596
     
    589605        cldfra(i) = MAX(0., cldfra(i) - contfra(i))
    590606        qcld(i) = MAX(0., qcld(i) - qcont(i))
    591         qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qsat(i) * contfra(i)))
    592 
    593         !--The contrail is always adjusted to saturation
    594         qiceincld = ( qcont(i) / contfra(i) - qsat(i) )
     607        qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qvcont(i)))
     608
     609        qvapincld = qvcont(i) / contfra(i)
     610        qiceincld = ( qcont(i) / contfra(i) - qvapincld )
    595611        !--The in-cloud ice crystals concentration is conserved
    596612        Niceincld = Ncont(i) / contfra(i)
    597613
    598614        !--If the ice water content is too low, the cloud is purely sublimated
    599         IF ( qiceincld .LT. eps ) THEN
     615        !--Most probably, we advected a cloud with no ice water content (possible
     616        !--if the entire cloud precipited for example)
     617        !--Same if the number concentration is too low (less than 1 #/m3)
     618        IF ( ( qiceincld .LT. eps ) .OR. ( ( Niceincld * rho(i) ) .LT. 1. ) ) THEN
    600619          dcfc_sub(i) = - contfra(i)
    601           dqic_sub(i) = - qiceincld * contfra(i)
    602           dqtc_sub(i) = - qcont(i)
     620          dqtc_sub(i) = - qvcont(i)
     621          dqic_sub(i) = - ( qcont(i) - qvcont(i) )
    603622          dnic_sub(i) = - Ncont(i)
    604        
    605         !--Only a part of the contrail is sublimated
     623
     624          contfra(i) = 0.
     625          qcont(i) = 0.
     626          qvcont(i) = 0.
     627          Ncont(i) = 0.
     628          clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcfc_sub(i))
     629          qclr(i) = qclr(i) - dqtc_sub(i) - dqic_sub(i)
     630
     631        !--Else, the cloud is adjusted and sublimated
    606632        ELSE
    607           !--Gamma distribution starting at qvapincld (everything that is below qiceincld_min)
    608           pdf_shape = nu_iwc_pdf_lscp / qiceincld
    609           pdf_y = pdf_shape * qiceincld_min(i)
    610           pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp      , pdf_y )
    611           pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y )
    612 
    613           !--Tendencies and diagnostics
    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
    618         ENDIF ! qiceincld .LT. eps
    619 
    620         !--Add tendencies
    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)
     633
     634          IF ( ok_unadjusted_contrails ) THEN
     635            IF ( qvapincld .GE. qsat(i) ) THEN
     636              !--If the cloud is initially supersaturated
     637              !--Exact explicit integration (qvcont exact, qice explicit)
     638              tauinv_depsub = depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub &
     639                  * (Niceincld / N_ice_volume)**(2./3.)
     640              qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub )
     641            ELSE
     642              !--If the cloud is initially subsaturated
     643              !--Exact explicit integration (qice exact, qvcont explicit)
     644              !--The barrier is set so that the resulting vapor in cloud
     645              !--cannot be greater than qsat
     646              !--qice_ratio is the ratio between the new ice content and
     647              !--the old one, it is comprised between 0 and 1
     648              tauinv_depsub = qiceincld**(1./3.) * kappa_depsub &
     649                  * (Niceincld / N_ice_volume)**(2./3.)
     650              qice_ratio = tauinv_depsub * dtime / 1.5 / qiceincld * ( qsat(i) - qvapincld )
     651              !--The new vapor in the cloud is increased with the
     652              !--sublimated ice
     653              qvapincld_new = qvapincld + qiceincld * ( 1. - MAX(0., 1. - qice_ratio)**1.5 )
     654              !--The new vapor in the cloud cannot be greater than qsat
     655              qvapincld_new = MIN(qvapincld_new, qsat(i))
     656              !--If all the ice is sublimated
     657              IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0.
     658            ENDIF ! qvapincld .GT. qsat
     659          ELSE
     660            !--We keep the saturation adjustment hypothesis, and the vapor in the
     661            !--cloud is set equal to the saturation vapor
     662            IF ( ( qvapincld + qiceincld ) .GT. qsat(i) ) THEN
     663              qvapincld_new = qsat(i)
     664            ELSE
     665              qvapincld_new = 0.
     666            ENDIF
     667          ENDIF ! ok_unadjusted_contrails
     668         
     669
     670          !---------------------------------------
     671          !--    DISSIPATION OF THE CONTRAIL    --
     672          !---------------------------------------
     673
     674          !--If the dissipation process must be activated
     675          IF ( ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) ) .GT. qvapincld ) THEN
     676            !--Gamma distribution starting at qvapincld
     677            pdf_shape = nu_iwc_pdf_lscp / qiceincld
     678            pdf_y = pdf_shape * ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) - qvapincld )
     679            pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp      , pdf_y )
     680            pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y )
     681
     682            !--Tendencies and diagnostics
     683            dcfc_sub(i) = - contfra(i) * pdf_e1
     684            dqic_sub(i) = - contfra(i) * pdf_e2 / pdf_shape
     685            dqtc_sub(i) = dcfc_sub(i) * qvapincld
     686            dnic_sub(i) = dcfc_sub(i) * Niceincld
     687
     688            !--Add tendencies
     689            contfra(i) = MAX(0., contfra(i) + dcfc_sub(i))
     690            qcont(i)   = qcont(i) + dqtc_sub(i) + dqic_sub(i)
     691            qvcont(i)  = qvcont(i) + dqtc_sub(i)
     692            Ncont(i)   = Ncont(i) + dnic_sub(i)
     693            clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcfc_sub(i))
     694            qclr(i)   = qclr(i) - dqtc_sub(i) - dqic_sub(i)
     695          ELSEIF ( qvapincld_new .EQ. 0. ) THEN
     696            !--If all the ice has been sublimated, we sublimate
     697            !--completely the cloud and do not activate the dissipation
     698            !--process
     699            !--Tendencies and diagnostics
     700            dcfc_sub(i) = - contfra(i)
     701            dqtc_sub(i) = - qvcont(i)
     702            dqic_sub(i) = - ( qcont(i) - qvcont(i) )
     703            dnic_sub(i) = - Ncont(i)
     704
     705            !--Add tendencies
     706            contfra(i) = 0.
     707            qcont(i) = 0.
     708            qvcont(i) = 0.
     709            Ncont(i) = 0.
     710            clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcfc_sub(i))
     711            qclr(i) = qclr(i) - dqtc_sub(i) - dqic_sub(i)
     712          ENDIF ! ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) ) .GT. qvapincld
     713
     714
     715          !------------------------------------
     716          !--       PHASE ADJUSTMENT        --
     717          !------------------------------------
     718
     719          IF ( qvapincld_new .GT. 0. ) THEN
     720            !--Adjustment of the IWC to the new vapor in cloud
     721            !--(this can be either positive or negative)
     722            dqtc_adj(i) = ( qvapincld_new * contfra(i) - qvcont(i) )
     723            dqic_adj(i) = - dqtc_adj(i)
     724
     725            !--Add tendencies
     726            !--The vapor in the cloud is updated, but not qcont as it is constant
     727            !--through this process, as well as contfra which is unmodified
     728            qvcont(i) = MAX(0., MIN(qcont(i), qvcont(i) + dqtc_adj(i)))
     729          ENDIF
     730
     731        ENDIF   ! qiceincld .LT. eps
    626732      ENDIF ! contfra(i) .GT. eps
    627733
     
    9751081          !--Whatever we do, the cloud will increase in size
    9761082          dcf_mix(i)  = clrfra_mix
    977           dqvc_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
     1083          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     1084            dqvc_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
     1085          ELSE
     1086            dqvc_mix(i) = clrfra_mix * qsat(i)
     1087            dqi_mix(i)  = clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) )
     1088          ENDIF
    9781089        ELSE
    9791090          !--We then calculate the clear sky part where the humidity is lower than
     
    10271138              dcf_mix(i)  = clrfra_mix * sigma_mix
    10281139              dqvc_mix(i) = clrfra_mix * sigma_mix * qsat(i)
    1029               dqi_mix(i)  = clrfra_mix * sigma_mix * (pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
     1140              dqi_mix(i)  = clrfra_mix * sigma_mix &
     1141                          * (pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
    10301142            ENDIF
    10311143          ENDIF
     
    11371249        !--diagnostic, and if the cloud size is increased, we add the new vapor to the
    11381250        !--cloud's vapor without condensing or sublimating ice crystals
    1139         qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
    1140                       - qcont(i) / contfra(i) * cldfra_mix / clrfra_mix
     1251        IF ( ok_unadjusted_contrails ) THEN
     1252          qiceinmix = ( qcont(i) - qvcont(i) ) / contfra(i) / ( 1. + clrfra_mix / cldfra_mix )
     1253          tauinv_depsub = qiceinmix**(1./3.) * kappa_depsub &
     1254              * (Niceincld / N_ice_volume)**(2./3.)
     1255          !qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime * tauinv_depsub ) )
     1256          qvapinmix_lim = qsat(i) - qiceinmix * MAX(1., 1.5 / ( dtime * tauinv_depsub ))
     1257          qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) &
     1258                        - qvcont(i) / contfra(i) * cldfra_mix / clrfra_mix
     1259        ELSE
     1260          !--NB. if tau_depsub = 0 (ie tauinv_depsub = inf), we get the same result as above
     1261          qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
     1262                        - qcont(i) / contfra(i) * cldfra_mix / clrfra_mix
     1263        ENDIF
    11411264
    11421265        IF ( qvapinclr_lim .LT. 0. ) THEN
    11431266          !--Whatever we do, contrails will increase in size
    11441267          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) )
     1268          IF ( ok_unadjusted_contrails ) THEN
     1269            dqtc_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
     1270          ELSE
     1271            dqtc_mix(i) = clrfra_mix * qsat(i)
     1272            dqic_mix(i) = clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) )
     1273          ENDIF
    11471274        ELSE
    11481275          !--We then calculate the clear sky part where the humidity is lower than
     
    11871314
    11881315          IF ( pdf_fra_above_lim .GT. eps ) THEN
    1189             dcfc_mix(i) = dcfc_mix(i) + clrfra_mix * sigma_mix
    1190             dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * sigma_mix &
    1191                         * pdf_q_above_lim / pdf_fra_above_lim
    1192             dqic_mix(i) = dqic_mix(i) + clrfra_mix * sigma_mix &
    1193                         * ( pdf_q_above_lim / pdf_fra_above_lim - qsat(i) )
     1316            IF ( ok_unadjusted_contrails ) THEN
     1317              dcfc_mix(i) = clrfra_mix * sigma_mix
     1318              dqtc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
     1319            ELSE
     1320              dcfc_mix(i) = clrfra_mix * sigma_mix
     1321              dqtc_mix(i) = clrfra_mix * sigma_mix * qsat(i)
     1322              dqic_mix(i) = clrfra_mix * sigma_mix &
     1323                          * ( pdf_q_above_lim / pdf_fra_above_lim - qsat(i) )
     1324            ENDIF
    11941325          ENDIF
    11951326
     
    11971328            dcfc_mix(i) = dcfc_mix(i) - cldfra_mix * ( 1. - sigma_mix )
    11981329            dqtc_mix(i) = dqtc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
    1199                         * qcont(i) / contfra(i)
     1330                        * qvcont(i) / contfra(i)
    12001331            dqic_mix(i) = dqic_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
    1201                         * ( qcont(i) / contfra(i) - qsat(i) )
     1332                        * ( qcont(i) - qvcont(i) ) / contfra(i)
    12021333            dnic_mix(i) = dnic_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
    1203                         * ( Ncont(i) / contfra(i) )
     1334                        * Ncont(i) / contfra(i)
    12041335          ENDIF
    12051336
     
    12231354        !--Add tendencies
    12241355        contfra(i) = contfra(i) + dcfc_mix(i)
    1225         qcont(i)   = qcont(i) + dqtc_mix(i)
     1356        qcont(i)   = qcont(i) + dqtc_mix(i) + dqic_mix(i)
     1357        qvcont(i)  = qvcont(i) + dqtc_mix(i)
    12261358        Ncont(i)   = Ncont(i) + dnic_mix(i)
    12271359        clrfra(i)  = clrfra(i) - dcfc_mix(i)
     
    12451377      IF ( ok_plane_contrail ) THEN
    12461378        IF ( contfra(i) .GT. eps ) THEN
    1247           mice = qcont(i) / MAX(eps, Ncont(i))
     1379          mice = ( qcont(i) - qvcont(i) ) / MAX(eps, Ncont(i))
    12481380          icefall_velo = ICECRYSTAL_VELO(mice, temp(i), pplay(i))
    12491381          !--Effective radius (in meters)
    12501382          dei = (mice / rho_ice * 3. / 4. / RPI)**(1./3.) / eff2vol_radius_contrails
    1251           dei = MIN(1e-4, MAX(1e-6, dei))
    12521383          !--Effective radius to effective diameter
    12531384          dei = 8. / 3. / SQRT(3.) * dei
     
    13331464        !
    13341465        IF ( contfra(i) .GT. eps ) THEN
    1335           mice = qcont(i) / MAX(eps, Ncont(i))
     1466          mice = ( qcont(i) - qvcont(i) ) / MAX(eps, Ncont(i))
    13361467          icefall_velo = ICECRYSTAL_VELO(mice, temp(i), pplay(i))
    13371468
     
    13401471          dzsed_cont(i) = MIN(dz, icefall_velo * dtime) * coef_sed
    13411472          cfsed_cont(i) = contfra(i)
    1342           qice_sedim = ( qcont(i) - qsat(i) * contfra(i) ) * dzsed_cont(i) / dz
     1473          qice_sedim = ( qcont(i) - qvcont(i) ) * dzsed_cont(i) / dz
    13431474          Nice_sedim = Ncont(i) * dzsed_cont(i) / dz
    13441475          !--Tendencies
    13451476          dcfc_sed1 = - contfra(i) * dzsed_cont(i) / dz
    13461477          dqic_sed1 = - qice_sedim
    1347           dqtc_sed1 = - qcont(i) * dzsed_cont(i) / dz
     1478          dqtc_sed1 = - qvcont(i) * dzsed_cont(i) / dz
    13481479          dnic_sed1 = - Nice_sedim
    13491480          !--Convert to flux
     
    13541485          coef_auto = MAX(0., MIN(1., icefall_velo * dtime / dz - 1.))
    13551486          dz_auto_cont = MIN(dz, icefall_velo * dtime) * coef_auto
    1356           qice_auto = ( qcont(i) - qsat(i) * contfra(i) ) * dz_auto_cont / dz
     1487          qice_auto = ( qcont(i) - qvcont(i) ) * dz_auto_cont / dz
    13571488          !--Tendencies
    13581489          dcfc_auto(i) = - contfra(i) * dz_auto_cont / dz
    13591490          dqic_auto(i) = - qice_auto
    1360           dqtc_auto(i) = - qcont(i) * dz_auto_cont / dz
     1491          dqtc_auto(i) = - qvcont(i) * dz_auto_cont / dz
    13611492          dnic_auto(i) = - Ncont(i) * dz_auto_cont / dz
    13621493          !--Convert to flux
    13631494          flauto(i) = flauto(i) + qice_auto * ( paprsdn(i) - paprsup(i) ) / RG / dtime
    13641495
     1496          !--Save the vapor in the cloud (will be needed later)
     1497          qvapincont = qvcont(i) / contfra(i)
    13651498          !--Add tendencies
    13661499          contfra(i)  = contfra(i) + dcfc_sed1 + dcfc_auto(i)
    1367           qcont(i)    = qcont(i) + dqtc_sed1 + dqtc_auto(i)
     1500          qcont(i)    = qcont(i) + dqtc_sed1 + dqtc_auto(i) + dqic_sed1 + dqic_auto(i)
     1501          qvcont(i)   = qvcont(i) + dqtc_sed1 + dqtc_auto(i)
    13681502          Ncont(i)    = Ncont(i) + dnic_sed1 + dnic_auto(i)
    13691503          clrfra(i)  = clrfra(i) - dcfc_sed1 - dcfc_auto(i)
    1370           qclr(i)    = qclr(i) - dqtc_sed1 - dqtc_auto(i)
     1504          qclr(i)    = qclr(i) - dqtc_sed1 - dqtc_auto(i) - dqic_sed1 - dqic_auto(i)
    13711505        ENDIF
    13721506        !
     
    15001634            IF ( sedfra1 .GT. eps ) THEN
    15011635              IF ( clrfra(i) .GT. eps ) THEN
    1502                 qvapinclr_lim = qsat(i) - qiceincld
     1636                IF ( ok_unadjusted_contrails ) THEN
     1637                  tauinv_depsub = qiceincld**(1./3.) * kappa_depsub &
     1638                      * (Niceincld / N_ice_volume)**(2./3.)
     1639                  qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) )
     1640                ELSE
     1641                  qvapinclr_lim = qsat(i) - qiceincld
     1642                ENDIF
    15031643
    15041644                rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     
    15331673
    15341674                IF ( pdf_fra_above_lim .GT. eps ) THEN
    1535                   dcfc_sed2 = dcfc_sed2 + sedfra_tmp
    1536                   dqtc_sed2 = dqtc_sed2 + sedfra_tmp &
    1537                                    * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim)
    1538                   dqic_sed2 = dqic_sed2 + sedfra_tmp &
     1675                  IF ( ok_unadjusted_contrails ) THEN
     1676                    dcfc_sed2 = dcfc_sed2 + sedfra_tmp
     1677                    dqtc_sed2 = dqtc_sed2 + sedfra_tmp * pdf_q_above_lim / pdf_fra_above_lim
     1678                    dqic_sed2 = dqic_sed2 + sedfra_tmp * qiceincld
     1679                    dnic_sed2 = dnic_sed2 + sedfra_tmp * Niceincld
     1680                  ELSE
     1681                    dcfc_sed2 = dcfc_sed2 + sedfra_tmp
     1682                    dqtc_sed2 = dqtc_sed2 + sedfra_tmp * qsat(i)
     1683                    dqic_sed2 = dqic_sed2 + sedfra_tmp &
    15391684                                   * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
    1540                   dnic_sed2 = dnic_sed2 + sedfra_tmp * Niceincld
     1685                    dnic_sed2 = dnic_sed2 + sedfra_tmp * Niceincld
     1686                  ENDIF
    15411687                ENDIF
    15421688              ENDIF ! clrfra(i) .GT. eps
     
    15561702            IF ( contfra(i) .GT. eps ) THEN
    15571703              dcfc_sed2 = dcfc_sed2 + sedfra3
    1558               dqtc_sed2 = dqtc_sed2 + sedfra3 * (qsat(i) + qiceincld)
     1704              dqtc_sed2 = dqtc_sed2 + sedfra3 * qvapincont
    15591705              dqic_sed2 = dqic_sed2 + sedfra3 * qiceincld
    15601706              dnic_sed2 = dnic_sed2 + sedfra3 * Niceincld
     
    15791725          !--Add tendencies
    15801726          contfra(i) = contfra(i) + dcfc_sed2
    1581           qcont(i)   = qcont(i)   + dqtc_sed2
     1727          qcont(i)   = qcont(i)   + dqtc_sed2 + dqic_sed2
     1728          qvcont(i)  = qvcont(i)  + dqtc_sed2
    15821729          Ncont(i)   = Ncont(i)   + dnic_sed2
    1583           clrfra(i)  = clrfra(i) - dcfc_sed2
    1584           qclr(i)    = qclr(i)   - dqtc_sed2
     1730          clrfra(i)  = clrfra(i)  - dcfc_sed2
     1731          qclr(i)    = qclr(i)    - dqtc_sed2 - dqic_sed2
    15851732          !--We re-include sublimated sedimentated ice into clear sky water vapor
    15861733          qclr(i)    = qclr(i) + qised_cont_abv(i)
     
    15991746        cldfra(i) = cldfra(i) + contfra(i)
    16001747        qcld(i) = qcld(i) + qcont(i)
    1601         qvc(i) = qvc(i) + qsat(i) * contfra(i)
     1748        qvc(i) = qvc(i) + qvcont(i)
    16021749      ENDIF
    16031750
     
    16961843      dcfc_ini(i) = MIN(dcfc_ini(i), issrfra(i))
    16971844      dqtc_ini(i) = MIN(dqtc_ini(i), qissr(i))
     1845      dqic_ini(i) = MIN(dqic_ini(i) + dqtc_ini(i), qissr(i)) - dqtc_ini(i)
    16981846
    16991847      !--Add tendencies
    17001848      cldfra(i)  = cldfra(i) + dcfc_ini(i)
    1701       qcld(i)    = qcld(i) + dqtc_ini(i)
    1702       qvc(i)     = qvc(i) + dcfc_ini(i) * qsat(i)
     1849      qcld(i)    = qcld(i) + dqtc_ini(i) + dqic_ini(i)
     1850      qvc(i)     = qvc(i) + dqtc_ini(i)
    17031851      issrfra(i) = issrfra(i) - dcfc_ini(i)
    1704       qissr(i)   = qissr(i) - dqtc_ini(i)
     1852      qissr(i)   = qissr(i) - dqtc_ini(i) - dqic_ini(i)
    17051853      clrfra(i)  = clrfra(i) - dcfc_ini(i)
    1706       qclr(i)    = qclr(i) - dqtc_ini(i)
     1854      qclr(i)    = qclr(i) - dqtc_ini(i) - dqic_ini(i)
    17071855
    17081856      contfra(i) = contfra(i) + dcfc_ini(i)
    1709       qcont(i)   = qcont(i) + dqtc_ini(i)
     1857      qcont(i)   = qcont(i) + dqtc_ini(i) + dqic_ini(i)
     1858      qvcont(i)  = qvcont(i) + dqtc_ini(i)
    17101859      Ncont(i)   = Ncont(i) + dnic_ini(i)
    17111860
     
    17151864      !-------------------------------------------
    17161865
    1717       IF ( (contfra(i) .LT. eps) .OR. (qcont(i) .LT. (qsat(i) * contfra(i))) &
    1718               & .OR. (Ncont(i) .LT. eps) ) THEN
     1866      IF ( (contfra(i) .LT. eps) .OR. (qcont(i) .LT. qvcont(i)) .OR. (Ncont(i) .LT. eps) ) THEN
    17191867        cldfra(i) = cldfra(i) - contfra(i)
    17201868        qcld(i) = qcld(i) - qcont(i)
    1721         qvc(i) = qvc(i) - qsat(i) * contfra(i)
     1869        qvc(i) = qvc(i) - qvcont(i)
    17221870        contfra(i) = 0.
    17231871        qcont(i) = 0.
     1872        qvcont(i) = 0.
    17241873        Ncont(i) = 0.
    17251874      ENDIF
     
    17321881        contfra(i)= 0.
    17331882        qcont(i)  = 0.
     1883        qvcont(i) = 0.
    17341884        Ncont(i)  = 0.
    17351885        qincld(i) = qsat(i)
     
    17401890      cldfra(i) = MAX(cldfra(i), contfra(i))
    17411891      qcld(i) = MAX(qcld(i), qcont(i))
    1742       qvc(i) = MAX(qvc(i), qsat(i) * contfra(i))
     1892      qvc(i) = MAX(qvc(i), qvcont(i))
    17431893
    17441894      !--Diagnostics
  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90

    r5796 r5797  
    232232  !$OMP THREADPRIVATE(ok_plane_contrail)
    233233
    234   INTEGER, SAVE, PROTECTED :: iflag_cross_sec_contrail=0      ! choice of the initial cross section parameterization
     234  LOGICAL, SAVE, PROTECTED :: ok_unadjusted_contrails=.FALSE.! if .TRUE., contrails are not adjusted to saturation. Requires an additional tracer
     235  !$OMP THREADPRIVATE(ok_unadjusted_contrails)
     236
     237  INTEGER, SAVE, PROTECTED :: iflag_cross_sec_contrail=0     ! choice of the initial cross section parameterization
    235238  !$OMP THREADPRIVATE(iflag_cross_sec_contrail)
    236239
     
    616619    CALL getin_p('chi_mixing',chi_mixing)
    617620    ! for aviation
     621    CALL getin_p('ok_unadjusted_contrails',ok_unadjusted_contrails)
    618622    CALL getin_p('iflag_cross_sec_contrail',iflag_cross_sec_contrail)
    619623    CALL getin_p('iflag_AEI_contrail',iflag_AEI_contrail)
     
    752756    WRITE(lunout,*) 'lscp_ini, chi_mixing:', chi_mixing
    753757    ! for aviation
     758    WRITE(lunout,*) 'lscp_ini, ok_unadjusted_contrails:', ok_unadjusted_contrails
    754759    WRITE(lunout,*) 'lscp_ini, iflag_cross_sec_contrail:', iflag_cross_sec_contrail
    755760    WRITE(lunout,*) 'lscp_ini, iflag_AEI_contrail:', iflag_AEI_contrail
  • LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_main.f90

    r5796 r5797  
    3030     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
    3131     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
    32      cfc_seri, qtc_seri, nic_seri,                      &
     32     cfc_seri, qtc_seri, qic_seri, nic_seri,            &
    3333     qice_cont, flight_dist, flight_fuel,               &
    3434     contfrarad, qradice_cont,                          &
     
    129129USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
    130130USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato, ok_ice_sedim
    131 USE lmdz_lscp_ini, ONLY : ok_plane_contrail
     131USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_unadjusted_contrails
    132132USE lmdz_lscp_ini, ONLY : ok_nodeep_lscp, ok_higher_cirrus_cover
    133133USE lmdz_lscp_ini, ONLY : ok_lscp_mergecond, gamma_mixth
     
    141141USE phys_local_var_mod, ONLY : dcfc_ini, dqic_ini, dqtc_ini, dnic_ini
    142142USE phys_local_var_mod, ONLY : dcfc_sub, dqic_sub, dqtc_sub, dnic_sub
    143 USE phys_local_var_mod, ONLY : dcfc_mix, dqic_mix, dqtc_mix, dnic_mix, dnic_agg
     143USE phys_local_var_mod, ONLY : dcfc_mix, dqic_mix, dqtc_mix, dnic_mix
     144USE phys_local_var_mod, ONLY : dnic_agg, dqic_adj, dqtc_adj
    144145USE phys_local_var_mod, ONLY : dcfc_sed, dqic_sed, dqtc_sed, dnic_sed
    145146USE phys_local_var_mod, ONLY : dcfc_auto, dqic_auto, dqtc_auto, dnic_auto
    146147USE phys_local_var_mod, ONLY : dcf_auto, dqi_auto, dqvc_auto
    147148USE phys_local_var_mod, ONLY : nice_ygcont, iwc_ygcont, rvol_ygcont, tau_ygcont
     149USE phys_local_var_mod, ONLY : nice_vscont, iwc_vscont, rvol_vscont, tau_vscont
    148150USE phys_local_var_mod, ONLY : nice_cont, iwc_cont, rvol_cont, tau_cont
    149151
     
    219221  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfc_seri         ! contrail fraction [-]
    220222  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qtc_seri         ! contrail total specific humidity [kg/kg]
     223  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qic_seri         ! contrail ice specific humidity [kg/kg]
    221224  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: nic_seri         ! contrail ice crystals concentration [#/kg]
    222225  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown concentration [m/s/m3]
     
    367370  REAL :: delta_z, deepconv_coef
    368371  ! for contrails
    369   REAL, DIMENSION(klon) :: contfra, qcont, Ncont
     372  REAL, DIMENSION(klon) :: contfra, qcont, qvcont, Ncont
    370373  REAL, DIMENSION(klon) :: totfra_in, qtot_in
    371374  LOGICAL, DIMENSION(klon) :: pt_pron_clds
     
    489492dqtc_sub(:,:)  = 0.
    490493dnic_sub(:,:)  = 0.
     494dqic_adj(:,:)  = 0.
     495dqtc_adj(:,:)  = 0.
    491496dcfc_mix(:,:)  = 0.
    492497dqic_mix(:,:)  = 0.
     
    787792            contfra(:) = 0.
    788793            qcont(:)   = 0.
     794            qvcont(:)  = 0.
    789795            Ncont(:)   = 0.
    790796            dzsed_cont(:) = 0.
     
    858864                  cfc_seri(i,k) = cfc_seri(i,k) * deepconv_coef
    859865                  qtc_seri(i,k) = qtc_seri(i,k) * deepconv_coef
     866                  qic_seri(i,k) = qic_seri(i,k) * deepconv_coef
    860867                  nic_seri(i,k) = nic_seri(i,k) * deepconv_coef
    861868                ENDIF
     
    973980                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
    974981                        dqvc_sed(:,k), dqvc_auto(:,k), &
    975                         cfc_seri(:,k), qtc_seri(:,k), nic_seri(:,k), &
     982                        cfc_seri(:,k), qtc_seri(:,k), qic_seri(:,k), nic_seri(:,k), &
    976983                        flight_dist(:,k), flight_fuel(:,k), &
    977                         contfra, qcont, Ncont, &
     984                        contfra, qcont, qvcont, Ncont, &
    978985                        Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
    979986                        AEI_contrails(:,k), AEI_surv_contrails(:,k), &
     
    983990                        dcfc_ini(:,k), dqic_ini(:,k), dqtc_ini(:,k), dnic_ini(:,k), &
    984991                        dcfc_sub(:,k), dqic_sub(:,k), dqtc_sub(:,k), dnic_sub(:,k), &
    985                         dcfc_mix(:,k), dqic_mix(:,k), dqtc_mix(:,k), dnic_mix(:,k), dnic_agg(:,k), &
     992                        dcfc_mix(:,k), dqic_mix(:,k), dqtc_mix(:,k), dnic_mix(:,k), &
     993                        dnic_agg(:,k), dqic_adj(:,k), dqtc_adj(:,k), &
    986994                        dcfc_sed(:,k), dqic_sed(:,k), dqtc_sed(:,k), dnic_sed(:,k), &
    987995                        dcfc_auto(:,k), dqic_auto(:,k), dqtc_auto(:,k), dnic_auto(:,k))
     
    11131121                        ENDIF
    11141122                       
    1115                         IF ( ok_unadjusted_clouds ) THEN
     1123                        IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN
    11161124                          !--AB We relax the saturation adjustment assumption
    11171125                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
     
    11711179                    zqn(i) = zq(i)
    11721180                    rneb(i,k) = 1.0
    1173                     IF ( ok_unadjusted_clouds ) THEN
     1181                    IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN
    11741182                      !--AB We relax the saturation adjustment assumption
    11751183                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
     
    11801188                    rhcl(i,k)=1.0
    11811189                ELSE
    1182                     IF ( ok_unadjusted_clouds ) THEN
     1190                    IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN
    11831191                      !--AB We relax the saturation adjustment assumption
    11841192     
     
    13751383          contfra(i) = 0.
    13761384          qcont(i) = 0.
     1385          qvcont(i) = 0.
    13771386          Ncont(i) = 0.
    13781387        ENDIF
    13791388      ENDDO
    13801389      cfc_seri(:,k) = contfra(:)
    1381       qtc_seri(:,k) = qcont(:)
    13821390      nic_seri(:,k) = Ncont(:)
    13831391      !--Ice water content of contrails
    1384       qice_cont(:,k) = qcont(:) - zqs(:) * contfra(:)
     1392      qice_cont(:,k) = qcont(:) - qvcont(:)
     1393      IF ( ok_unadjusted_contrails ) THEN
     1394        qtc_seri(:,k) = qvcont(:)
     1395        qic_seri(:,k) = qice_cont(:,k)
     1396      ELSE
     1397        qtc_seri(:,k) = qcont(:)
     1398      ENDIF
    13851399      !--Radiative properties
    13861400      contfrarad(:,k) = contfra(:)
     
    15081522          rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
    15091523          iwp_cont = 1e3 * dqic_ini(i,k) / dcfc_ini(i,k) * rhodz
    1510           rei_cont = MIN(100., MAX(10., rvol_ygcont(i,k) / eff2vol_radius_contrails))
     1524          rei_cont = MIN(200., MAX(10., rvol_ygcont(i,k) / eff2vol_radius_contrails))
    15111525          tau_ygcont(i,k) = iwp_cont*(3.448e-3+2.431/rei_cont)
    15121526        ELSE
     
    15171531        ENDIF
    15181532        !--All contrails
    1519         IF ( cfc_seri(i,k) .GT. 1e-3 ) THEN
     1533        IF ( cfc_seri(i,k) .GT. 1e-8 ) THEN
    15201534          rho = pplay(i,k) / zt(i) / RD
    15211535          nice_cont(i,k) = nic_seri(i,k) / cfc_seri(i,k) / 1e6 * rho
     
    15251539          rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
    15261540          iwp_cont = 1e3 * qice_cont(i,k) / contfrarad(i,k) * rhodz
    1527           rei_cont = MIN(100., MAX(10., rvol_cont(i,k) / eff2vol_radius_contrails))
     1541          rei_cont = MIN(200., MAX(10., rvol_cont(i,k) / eff2vol_radius_contrails))
    15281542          tau_cont(i,k) = iwp_cont*(3.448e-3+2.431/rei_cont)
     1543          IF ( tau_cont(i,k) .GT. 0.05 ) THEN
     1544            nice_vscont(i,k) = nice_cont(i,k)
     1545            iwc_vscont(i,k) = iwc_cont(i,k)
     1546            rvol_vscont(i,k) = rvol_cont(i,k)
     1547            tau_vscont(i,k) = tau_cont(i,k)
     1548          ELSE
     1549            nice_vscont(i,k) = missing_val
     1550            iwc_vscont(i,k) = missing_val
     1551            rvol_vscont(i,k) = missing_val
     1552            tau_vscont(i,k) = missing_val
     1553          ENDIF
    15291554        ELSE
    15301555          nice_cont(i,k) = missing_val
     
    15321557          rvol_cont(i,k) = missing_val
    15331558          tau_cont(i,k) = missing_val
     1559          nice_vscont(i,k) = missing_val
     1560          iwc_vscont(i,k) = missing_val
     1561          rvol_vscont(i,k) = missing_val
     1562          tau_vscont(i,k) = missing_val
    15341563        ENDIF
    15351564      ENDDO
  • LMDZ6/branches/contrails/libf/phylmd/phyetat0_mod.f90

    r5791 r5797  
    2626       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
    2727       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, &
    28        cf_ancien, qvc_ancien, qvcon, qccon, cfc_ancien, qtc_ancien, nic_ancien, &
     28       cf_ancien, qvc_ancien, qvcon, qccon, cfc_ancien, qtc_ancien, qic_ancien, nic_ancien, &
    2929       tke_ancien, radpas, radsol, rain_fall, ratqs, &
    3030       rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
     
    498498    ancien_ok=ancien_ok.AND.phyetat0_get(cfc_ancien,"CFCANCIEN","CFCANCIEN",0.)
    499499    ancien_ok=ancien_ok.AND.phyetat0_get(qtc_ancien,"QTCANCIEN","QTCANCIEN",0.)
     500    ancien_ok=ancien_ok.AND.phyetat0_get(qic_ancien,"QICANCIEN","QICANCIEN",0.)
    500501    ancien_ok=ancien_ok.AND.phyetat0_get(nic_ancien,"NICANCIEN","NICANCIEN",0.)
    501502  ELSE
    502503    cfc_ancien(:,:)=0.
    503504    qtc_ancien(:,:)=0.
     505    qic_ancien(:,:)=0.
    504506    nic_ancien(:,:)=0.
    505507  ENDIF
     
    535537    IF ( ( maxval(cfc_ancien).EQ.minval(cfc_ancien) ) .OR. &
    536538         ( maxval(qtc_ancien).EQ.minval(qtc_ancien) ) .OR. &
     539         ( maxval(qic_ancien).EQ.minval(qic_ancien) ) .OR. &
    537540         ( maxval(nic_ancien).EQ.minval(nic_ancien) ) ) THEN
    538541       ancien_ok=.false.
  • LMDZ6/branches/contrails/libf/phylmd/phyredem.f90

    r5791 r5797  
    2424                                ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
    2525                                qvcon, qccon,                                &
    26                                 qvc_ancien, cfc_ancien, qtc_ancien, nic_ancien, &
     26                                qvc_ancien, cfc_ancien, qtc_ancien,          &
     27                                qic_ancien, nic_ancien,                      &
    2728                                u_ancien, v_ancien, tke_ancien,              &
    2829                                clwcon, rnebcon, ratqs, pbl_tke,             &
     
    312313      CALL put_field(pass,"CFCANCIEN", "CFCANCIEN", cfc_ancien)
    313314      CALL put_field(pass,"QTCANCIEN", "QTCANCIEN", qtc_ancien)
     315      CALL put_field(pass,"QICANCIEN", "QICANCIEN", qic_ancien)
    314316      CALL put_field(pass,"NICANCIEN", "NICANCIEN", nic_ancien)
    315317    ENDIF
  • LMDZ6/branches/contrails/libf/phylmd/phys_local_var_mod.F90

    r5796 r5797  
    710710      REAL, SAVE, ALLOCATABLE :: qtc_seri(:,:), d_qtc_dyn(:,:)
    711711      !$OMP THREADPRIVATE(qtc_seri, d_qtc_dyn)
     712      REAL, SAVE, ALLOCATABLE :: qic_seri(:,:), d_qic_dyn(:,:)
     713      !$OMP THREADPRIVATE(qic_seri, d_qic_dyn)
    712714      REAL, SAVE, ALLOCATABLE :: nic_seri(:,:), d_nic_dyn(:,:)
    713715      !$OMP THREADPRIVATE(nic_seri, d_nic_dyn)
     
    726728      REAL, SAVE, ALLOCATABLE :: rvol_ygcont(:,:), tau_ygcont(:,:)
    727729      !$OMP THREADPRIVATE(rvol_ygcont, tau_ygcont)
     730      REAL, SAVE, ALLOCATABLE :: nice_vscont(:,:), iwc_vscont(:,:)
     731      !$OMP THREADPRIVATE(nice_vscont, iwc_vscont)
     732      REAL, SAVE, ALLOCATABLE :: rvol_vscont(:,:), tau_vscont(:,:)
     733      !$OMP THREADPRIVATE(rvol_vscont, tau_vscont)
    728734      REAL, SAVE, ALLOCATABLE :: nice_cont(:,:), iwc_cont(:,:)
    729735      !$OMP THREADPRIVATE(nice_cont, iwc_cont)
     
    744750      REAL, SAVE, ALLOCATABLE :: dcfc_mix(:,:), dqic_mix(:,:), dqtc_mix(:,:), dnic_mix(:,:)
    745751      !$OMP THREADPRIVATE(dcfc_mix, dqic_mix, dqtc_mix, dnic_mix)
    746       REAL, SAVE, ALLOCATABLE :: dnic_agg(:,:)
    747       !$OMP THREADPRIVATE(dnic_agg)
     752      REAL, SAVE, ALLOCATABLE :: dnic_agg(:,:), dqic_adj(:,:), dqtc_adj(:,:)
     753      !$OMP THREADPRIVATE(dnic_agg, dqic_adj, dqtc_adj)
    748754      REAL, SAVE, ALLOCATABLE :: cldfra_cont(:,:), cldtau_nocont(:,:), cldemi_nocont(:,:)
    749755      !$OMP THREADPRIVATE(cldfra_cont, cldtau_nocont, cldemi_nocont)
     
    13351341      ALLOCATE(cfc_seri(klon,klev), d_cfc_dyn(klon,klev))
    13361342      ALLOCATE(qtc_seri(klon,klev), d_qtc_dyn(klon,klev))
     1343      ALLOCATE(qic_seri(klon,klev), d_qic_dyn(klon,klev))
    13371344      ALLOCATE(nic_seri(klon,klev), d_nic_dyn(klon,klev))
    13381345      ALLOCATE(flight_dist(klon,klev), flight_fuel(klon,klev))
     
    13431350      ALLOCATE(nice_ygcont(klon,klev), iwc_ygcont(klon,klev))
    13441351      ALLOCATE(rvol_ygcont(klon,klev), tau_ygcont(klon,klev))
     1352      ALLOCATE(nice_vscont(klon,klev), iwc_vscont(klon,klev))
     1353      ALLOCATE(rvol_vscont(klon,klev), tau_vscont(klon,klev))
    13451354      ALLOCATE(nice_cont(klon,klev), iwc_cont(klon,klev))
    13461355      ALLOCATE(rvol_cont(klon,klev), tau_cont(klon,klev))
     
    13501359      ALLOCATE(dcfc_sub(klon,klev), dqic_sub(klon,klev), dqtc_sub(klon,klev), dnic_sub(klon,klev))
    13511360      ALLOCATE(dcfc_mix(klon,klev), dqic_mix(klon,klev), dqtc_mix(klon,klev), dnic_mix(klon,klev))
    1352       ALLOCATE(dnic_agg(klon,klev))
     1361      ALLOCATE(dnic_agg(klon,klev), dqic_adj(klon,klev), dqtc_adj(klon,klev))
    13531362      ALLOCATE(dcfc_sed(klon,klev), dqic_sed(klon,klev), dqtc_sed(klon,klev), dnic_sed(klon,klev))
    13541363      ALLOCATE(dcfc_auto(klon,klev), dqic_auto(klon,klev), dqtc_auto(klon,klev), dnic_auto(klon,klev))
     
    17871796
    17881797!-- LSCP - aviation and contrails variables
    1789       DEALLOCATE(cfc_seri, d_cfc_dyn, qtc_seri, d_qtc_dyn, nic_seri, d_nic_dyn)
     1798      DEALLOCATE(cfc_seri, d_cfc_dyn, qtc_seri, d_qtc_dyn)
     1799      DEALLOCATE(qic_seri, d_qic_dyn, nic_seri, d_nic_dyn)
    17901800      DEALLOCATE(d_q_avi, flight_dist, flight_fuel)
    17911801      DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
    17921802      DEALLOCATE(AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails)
    17931803      DEALLOCATE(nice_ygcont, iwc_ygcont, rvol_ygcont, tau_ygcont)
     1804      DEALLOCATE(nice_vscont, iwc_vscont, rvol_vscont, tau_vscont)
    17941805      DEALLOCATE(nice_cont, iwc_cont, rvol_cont, tau_cont)
    17951806      DEALLOCATE(qice_cont, contfra, qradice_cont)
     
    17971808      DEALLOCATE(dcfc_sub, dqic_sub, dqtc_sub, dnic_sub)
    17981809      DEALLOCATE(dcfc_mix, dqic_mix, dqtc_mix, dnic_mix)
    1799       DEALLOCATE(dnic_agg)
     1810      DEALLOCATE(dnic_agg, dqic_adj, dqtc_adj)
    18001811      DEALLOCATE(dcfc_sed, dqic_sed, dqtc_sed, dnic_sed)
    18011812      DEALLOCATE(dcfc_auto, dqic_auto, dqtc_auto, dnic_auto)
  • LMDZ6/branches/contrails/libf/phylmd/phys_output_ctrlout_mod.F90

    r5796 r5797  
    26462646  TYPE(ctrl_out), SAVE :: o_dqtcdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    26472647    'dqtcdyn', 'Dynamics contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
     2648  TYPE(ctrl_out), SAVE :: o_qicseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2649    'qicseri', 'Contrail ice specific humidity', 'kg/kg', (/ ('',i=1,10) /))
     2650  TYPE(ctrl_out), SAVE :: o_dqicdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2651    'dqicdyn', 'Dynamics contrail ice specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
    26482652  TYPE(ctrl_out), SAVE :: o_nicseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    26492653    'nicseri', 'Contrail ice crystals number concentration', '#/kg', (/ ('',i=1,10) /))
     
    27022706  TYPE(ctrl_out), SAVE :: o_dnicagg = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    27032707    'dnicagg', 'Aggregation contrail ice crystals concentration tendency', '#/kg/s', (/ ('', i=1, 10) /))
     2708  TYPE(ctrl_out), SAVE :: o_dqicadj = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2709    'dqicadj', 'Phase adjustment contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2710  TYPE(ctrl_out), SAVE :: o_dqtcadj = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2711    'dqtcadj', 'Phase adjustment contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    27042712  TYPE(ctrl_out), SAVE :: o_flight_dist = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    27052713    'flightdist', 'Aviation flown distance', 'm/s/m^3', (/ ('', i=1, 10) /))
     
    27222730  TYPE(ctrl_out), SAVE :: o_tau_ygcont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    27232731    'tau_ygcont', 'Young contrails optical depth', '-', (/ ('', i=1, 10) /))
     2732  TYPE(ctrl_out), SAVE :: o_nice_vscont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2733    'nice_vscont', 'Ice particle number concentration in visible contrails', '#/cm3', (/ ('', i=1, 10) /))
     2734  TYPE(ctrl_out), SAVE :: o_iwc_vscont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2735    'iwc_vscont', 'Ice water content in visible contrails', 'g/m3', (/ ('', i=1, 10) /))
     2736  TYPE(ctrl_out), SAVE :: o_rvol_vscont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2737    'rvol_vscont', 'Ice crystals volumic radius in visible contrails', 'microns', (/ ('', i=1, 10) /))
     2738  TYPE(ctrl_out), SAVE :: o_tau_vscont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2739    'tau_vscont', 'Visible contrails optical depth', '-', (/ ('', i=1, 10) /))
    27242740  TYPE(ctrl_out), SAVE :: o_nice_cont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    27252741    'nice_cont', 'Ice particle number concentration in contrails', '#/cm3', (/ ('', i=1, 10) /))
  • LMDZ6/branches/contrails/libf/phylmd/phys_output_write_mod.F90

    r5796 r5797  
    242242         o_issrfra250to300, o_issrfra300to400, o_issrfra400to500, &
    243243!-- LSCP - aviation variables
    244          o_cfcseri, o_dcfcdyn, o_qtcseri, o_dqtcdyn, o_nicseri, o_dnicdyn, o_dqavi, &
     244         o_cfcseri, o_dcfcdyn, o_qtcseri, o_dqtcdyn, &
     245         o_qicseri, o_dqicdyn, o_nicseri, o_dnicdyn, o_dqavi, &
    245246         o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
    246247         o_AEI_contrails, o_AEI_surv_contrails, o_fsurv_contrails, o_section_contrails, &
    247248         o_nice_ygcont, o_iwc_ygcont, o_rvol_ygcont, o_tau_ygcont, &
     249         o_nice_vscont, o_iwc_vscont, o_rvol_vscont, o_tau_vscont, &
    248250         o_nice_cont, o_iwc_cont, o_rvol_cont, o_tau_cont, &
    249251         o_flight_dist, o_flight_fuel, o_qice_cont, &
    250252         o_dcfcini, o_dqicini, o_dqtcini, o_dnicini, &
    251          o_dcfcmix, o_dqicmix, o_dqtcmix, o_dnicmix, o_dnicagg, &
     253         o_dcfcmix, o_dqicmix, o_dqtcmix, o_dnicmix, &
     254         o_dnicagg, o_dqicadj, o_dqtcadj, &
    252255         o_dcfcsub, o_dqicsub, o_dqtcsub, o_dnicsub, &
    253256         o_dcfcsed, o_dqicsed, o_dqtcsed, o_dnicsed, &
     
    410413         issrfra100to150, issrfra150to200, issrfra200to250, &
    411414         issrfra250to300, issrfra300to400, issrfra400to500, &
    412          cfc_seri, d_cfc_dyn, qtc_seri, d_qtc_dyn, nic_seri, d_nic_dyn, d_q_avi, &
     415         cfc_seri, d_cfc_dyn, qtc_seri, d_qtc_dyn, &
     416         qic_seri, d_qic_dyn, nic_seri, d_nic_dyn, d_q_avi, &
    413417         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    414418         AEI_contrails, AEI_surv_contrails, fsurv_contrails, section_contrails, &
    415419         nice_ygcont, iwc_ygcont, rvol_ygcont, tau_ygcont, &
     420         nice_vscont, iwc_vscont, rvol_vscont, tau_vscont, &
    416421         nice_cont, iwc_cont, rvol_cont, tau_cont, &
    417422         flight_dist, flight_fuel, qice_cont, &
    418423         dcfc_ini, dqic_ini, dqtc_ini, dnic_ini, &
    419424         dcfc_sub, dqic_sub, dqtc_sub, dnic_sub, &
    420          dcfc_mix, dqic_mix, dqtc_mix, dnic_mix, dnic_agg, &
     425         dcfc_mix, dqic_mix, dqtc_mix, dnic_mix, &
     426         dnic_agg, dqic_adj, dqtc_adj, &
    421427         dcfc_sed, dqic_sed, dqtc_sed, dnic_sed, &
    422428         dcfc_auto, dqic_auto, dqtc_auto, dnic_auto, &
     
    24552461         CALL histwrite_phy(o_qtcseri, qtc_seri)
    24562462         CALL histwrite_phy(o_dqtcdyn, d_qtc_dyn)
     2463         CALL histwrite_phy(o_qicseri, qic_seri)
     2464         CALL histwrite_phy(o_dqicdyn, d_qic_dyn)
    24572465         CALL histwrite_phy(o_nicseri, nic_seri)
    24582466         CALL histwrite_phy(o_dnicdyn, d_nic_dyn)
     
    24712479         CALL histwrite_phy(o_rvol_ygcont, rvol_ygcont)
    24722480         CALL histwrite_phy(o_tau_ygcont, tau_ygcont)
     2481         CALL histwrite_phy(o_nice_vscont, nice_vscont)
     2482         CALL histwrite_phy(o_iwc_vscont, iwc_vscont)
     2483         CALL histwrite_phy(o_rvol_vscont, rvol_vscont)
     2484         CALL histwrite_phy(o_tau_vscont, tau_vscont)
    24732485         CALL histwrite_phy(o_nice_cont, nice_cont)
    24742486         CALL histwrite_phy(o_iwc_cont, iwc_cont)
     
    24892501         CALL histwrite_phy(o_dnicmix, dnic_mix)
    24902502         CALL histwrite_phy(o_dnicagg, dnic_agg)
     2503         CALL histwrite_phy(o_dqicadj, dqic_adj)
     2504         CALL histwrite_phy(o_dqtcadj, dqtc_adj)
    24912505         CALL histwrite_phy(o_dcfcsed, dcfc_sed)
    24922506         CALL histwrite_phy(o_dqicsed, dqic_sed)
  • LMDZ6/branches/contrails/libf/phylmd/phys_state_var_mod.F90

    r5791 r5797  
    142142      REAL, ALLOCATABLE, SAVE :: qvcon(:,:), qccon(:,:)
    143143!$OMP THREADPRIVATE(qvcon, qccon)
    144       REAL, ALLOCATABLE, SAVE :: cfc_ancien(:,:), qtc_ancien(:,:), nic_ancien(:,:)
    145 !$OMP THREADPRIVATE(cfc_ancien, qtc_ancien, nic_ancien)
     144      REAL, ALLOCATABLE, SAVE :: cfc_ancien(:,:), qtc_ancien(:,:)
     145!$OMP THREADPRIVATE(cfc_ancien, qtc_ancien)
     146      REAL, ALLOCATABLE, SAVE :: qic_ancien(:,:), nic_ancien(:,:)
     147!$OMP THREADPRIVATE(qic_ancien, nic_ancien)
    146148       REAL, ALLOCATABLE, SAVE :: tke_ancien(:,:)
    147149!$OMP THREADPRIVATE(tke_ancien)
     
    680682      ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
    681683      ALLOCATE(cf_ancien(klon,klev), qvc_ancien(klon,klev))
    682       ALLOCATE(cfc_ancien(klon,klev), qtc_ancien(klon,klev), nic_ancien(klon,klev))
     684      ALLOCATE(cfc_ancien(klon,klev), qtc_ancien(klon,klev))
     685      ALLOCATE(qic_ancien(klon,klev), nic_ancien(klon,klev))
    683686      ALLOCATE(qvcon(klon,klev), qccon(klon,klev))
    684687      ALLOCATE(tke_ancien(klon,klev+1))
     
    951954      DEALLOCATE(u_ancien, v_ancien)
    952955      DEALLOCATE(cf_ancien, qvc_ancien)
    953       DEALLOCATE(cfc_ancien, qtc_ancien, nic_ancien)
     956      DEALLOCATE(cfc_ancien, qtc_ancien, qic_ancien, nic_ancien)
    954957      DEALLOCATE(qvcon, qccon)
    955958      DEALLOCATE(tke_ancien)
  • LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90

    r5796 r5797  
    3939    USE ioipsl_getin_p_mod, ONLY : getin_p
    4040    USE indice_sol_mod
    41     USE infotrac_phy, ONLY: nqtot, nbtr, nqo, nqtke, tracers, type_trac, addPhase, ivap, iliq, isol, ibs, icf, iqvc, itke, icfc, iqtc, inic
     41    USE infotrac_phy, ONLY: nqtot, nbtr, nqo, nqtke, tracers, type_trac, addPhase, &
     42        ivap, iliq, isol, ibs, icf, iqvc, itke, icfc, iqtc, iqic, inic
    4243    USE strings_mod,  ONLY: strIdx
    4344    USE iophy
     
    341342       dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    342343       !-- LSCP - aviation and contrails variables
    343        cfc_seri, qtc_seri, nic_seri, d_cfc_dyn, d_qtc_dyn, d_nic_dyn, &
     344       cfc_seri, qtc_seri, qic_seri, nic_seri, d_cfc_dyn, d_qtc_dyn, d_qic_dyn, d_nic_dyn, &
    344345       d_q_avi, flight_dist, flight_fuel, &
    345346       qice_cont, contfra, qradice_cont, &
     
    24912492          cfc_seri(i,k)= 0.
    24922493          qtc_seri(i,k)= 0.
     2494          qic_seri(i,k)= 0.
    24932495          nic_seri(i,k)= 0.
    24942496          !CR: ATTENTION, on rajoute la variable glace
     
    25082510             cfc_seri(i,k) = qx(i,k,icfc)
    25092511             qtc_seri(i,k) = qx(i,k,iqtc)
     2512             qic_seri(i,k) = qx(i,k,iqic)
    25102513             nic_seri(i,k) = qx(i,k,inic)
    25112514             !--DYNAMICO can return NaNs for children tracers
    25122515             IF (ISNAN(cfc_seri(i,k))) cfc_seri(i,k) = 0.
    25132516             IF (ISNAN(qtc_seri(i,k))) qtc_seri(i,k) = 0.
     2517             IF (ISNAN(qic_seri(i,k))) qic_seri(i,k) = 0.
    25142518             IF (ISNAN(nic_seri(i,k))) nic_seri(i,k) = 0.
    25152519          ENDIF
     
    26032607       d_cfc_dyn(:,:)= (cfc_seri(:,:)-cfc_ancien(:,:))/phys_tstep
    26042608       d_qtc_dyn(:,:)= (qtc_seri(:,:)-qtc_ancien(:,:))/phys_tstep
     2609       d_qic_dyn(:,:)= (qic_seri(:,:)-qic_ancien(:,:))/phys_tstep
    26052610       d_nic_dyn(:,:)= (nic_seri(:,:)-nic_ancien(:,:))/phys_tstep
    26062611       d_tke_dyn(:,:)= (pbl_tke(:,:,is_ave)-tke_ancien(:,:))/phys_tstep
     
    26282633       d_cfc_dyn(:,:)= 0.0
    26292634       d_qtc_dyn(:,:)= 0.0
     2635       d_qic_dyn(:,:)= 0.0
    26302636       d_nic_dyn(:,:)= 0.0
    26312637       d_tke_dyn(:,:)= 0.0
     
    27552761            qvc_seri(i,k) = qvc_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27562762            qtc_seri(i,k) = qtc_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2763            qic_seri(i,k) = qic_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27572764          ELSE
    27582765            ql_seri_lscp(i,k) = 0.
     
    27602767            qvc_seri(i,k) = 0.
    27612768            qtc_seri(i,k) = 0.
     2769            qic_seri(i,k) = 0.
    27622770          ENDIF
    27632771        ENDDO
     
    40334041        DO i = 1, klon
    40344042          qtc_seri(i,k) = qtc_seri(i,k) * q_seri(i,k)
     4043          qic_seri(i,k) = qic_seri(i,k) * q_seri(i,k)
    40354044        ENDDO
    40364045      ENDDO
     
    40634072         dqi_adj, dqi_sub, dqi_con, dqi_mix, &
    40644073         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    4065          cfc_seri, qtc_seri, nic_seri, qice_cont, &
     4074         cfc_seri, qtc_seri, qic_seri, nic_seri, qice_cont, &
    40664075         flight_dist, flight_fuel, contfra, qradice_cont, &
    40674076         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
     
    58235832             d_qx(i,k,icfc) = ( cfc_seri(i,k) - qx(i,k,icfc) ) / phys_tstep
    58245833             d_qx(i,k,iqtc) = ( qtc_seri(i,k) - qx(i,k,iqtc) ) / phys_tstep
     5834             d_qx(i,k,iqic) = ( qic_seri(i,k) - qx(i,k,iqic) ) / phys_tstep
    58255835             d_qx(i,k,inic) = ( nic_seri(i,k) - qx(i,k,inic) ) / phys_tstep
    58265836          ENDIF
     
    58665876    cfc_ancien(:,:)= cfc_seri(:,:)
    58675877    qtc_ancien(:,:)= qtc_seri(:,:)
     5878    qic_ancien(:,:)= qic_seri(:,:)
    58685879    nic_ancien(:,:)= nic_seri(:,:)
    58695880    tke_ancien(:,:)= pbl_tke(:,:,is_ave)
Note: See TracChangeset for help on using the changeset viewer.