Changeset 5601


Ignore:
Timestamp:
Apr 2, 2025, 4:04:40 PM (2 months ago)
Author:
aborella
Message:

Multiple changes:

  • tracers which were ratios are now absolute quantities. This is needed because when the ratio

is not defined, some aberrations may occur

  • added a new tracer for total specific humidity in contrails
  • rework of the mixing process for cirrus clouds (and contrails)
  • changed the numerical integration of ice crystals' sublimation
  • subroutines do not take real inputs anymore (at least klon tables)
  • added more radiative diagnostics for contrails
Location:
LMDZ6/branches/contrails
Files:
20 edited

Legend:

Unmodified
Added
Removed

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

    r5598 r5601  
    898898        <field id="cfseri"     long_name="Cloud fraction"    unit="-" />
    899899        <field id="dcfdyn"     long_name="Dynamics cloud fraction tendency"    unit="s-1" />
    900         <field id="rvcseri"    long_name="Cloudy water vapor to total water vapor ratio"    unit="-" />
    901         <field id="drvcdyn"    long_name="Dynamics cloudy water vapor to total water vapor ratio tendency"    unit="s-1" />
     900        <field id="qvcseri"    long_name="Cloudy water vapor"    unit="kg/kg" />
     901        <field id="dqvcdyn"    long_name="Dynamics cloudy water vapor tendency"    unit="kg/kg/s" />
    902902        <field id="qsub"       long_name="Subsaturated clear sky total water"    unit="kg/kg" />
    903903        <field id="qissr"      long_name="Supersaturated clear sky total water"    unit="kg/kg" />
    … …  
    920920        <field id="qsati"      long_name="Saturation with respect to ice"    unit="kg/kg" />
    921921
    922             <field id="rcontseri"  long_name="Linear contrail fraction to cloud fraction ratio" unit="-" />
    923             <field id="drcontdyn"  long_name="Dynamics linear contrail fraction to cloud fraction ratio tendency" unit="s-1" />
     922            <field id="cfaseri"    long_name="Linear contrail fraction" unit="-" />
     923            <field id="dcfadyn"    long_name="Dynamics linear contrail fraction tendency" unit="s-1" />
     924            <field id="qtaseri"    long_name="Linear contrail total specific humidity" unit="s-1" />
     925        <field id="dqtadyn"    long_name="Dynamics linear contrail total specific humidity tendency" unit="kg/kg/s" />
    924926            <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
    925927            <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
    926928        <field id="potcontfraP"  long_name="Fraction with pontential persistent contrail"    unit="-" />
    927929        <field id="potcontfraNP" long_name="Fraction with potential non-persistent contrail"    unit="-" />
    928         <field id="contfra"    long_name="Linear contrail fraction"    unit="-" />
    929         <field id="qice_cont"  long_name="Contrails ice specific humidity seen by radiation"    unit="kg/kg" />
    930         <field id="dcontfracir"  long_name="Linear contrail fraction to cirrus cloud fraction tendency"    unit="-" />
    931         <field id="dcfavi"     long_name="Aviation cloud fraction tendency"    unit="s-1" />
    932         <field id="dqiavi"     long_name="Aviation ice tendency"    unit="kg/kg/s" />
    933         <field id="dqvcavi"    long_name="Aviation cloudy water vapor tendency"    unit="kg/kg/s" />
     930        <field id="qice_cont"  long_name="Contrails ice specific humidity"    unit="kg/kg" />
     931        <field id="dcfaini"    long_name="Initial formation linear contrail fraction tendency"    unit="s-1" />
     932        <field id="dqiaini"    long_name="Initial formation linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
     933        <field id="dqtaini"    long_name="Initial formation linear contrail total specific humidity tendency"    unit="kg/kg/s" />
     934        <field id="dcfacir"    long_name="Conversion to cirrus linear contrail fraction tendency"    unit="s-1" />
     935        <field id="dqtacir"    long_name="Conversion to cirrus linear contrail total specific humidity tendency"    unit="kg/kg/s" />
     936        <field id="dcfasub"    long_name="Sublimation linear contrail fraction tendency"    unit="s-1" />
     937        <field id="dqiasub"    long_name="Sublimation linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
     938        <field id="dqtasub"    long_name="Sublimation linear contrail total specific humidity tendency"    unit="kg/kg/s" />
     939        <field id="dcfamix"    long_name="Mixing linear contrail fraction tendency"    unit="s-1" />
     940        <field id="dqiamix"    long_name="Mixing linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
     941        <field id="dqtamix"    long_name="Mixing linear contrail total specific humidity tendency"    unit="kg/kg/s" />
    934942        <field id="flightdist" long_name="Aviation flown distance concentration"    unit="m/s/m3" />
    935943        <field id="flighth2o"  long_name="Aviation emitted H2O concentration"    unit="kg H2O/s/m3" />
    … …  
    937945        <field id="cldfra_nocont" long_name="Cloud fraction w/o contrails"    unit="-" />
    938946        <field id="cldtau_nocont" long_name="Cloud optical thickness w/o contrails"    unit="1" />
     947        <field id="conttau" long_name="Contrails optical thickness"    unit="1" />
     948        <field id="contemi" long_name="Contrails optical emissivity"    unit="1" />
    939949        <field id="cldemi_nocont" long_name="Cloud optical emissivity w/o contrails"    unit="1" />
    940950        <field id="iwc_nocont"    long_name="Cloud ice water content seen by radiation w/o contrails"    unit="kg/m3" />

  • LMDZ6/branches/contrails/libf/dyn3d/conf_gcm.f90

    r5285 r5601  
    799799     CALL getin('type_trac',type_trac)
    800800
     801     !Config  Key  = adv_qsat_liq
     802     !Config  Desc = option for qsat calculation in the dynamics
     803     !Config  Def  = n
     804     !Config  Help = controls which phase is considered for qsat calculation
     805     !Config         
     806     adv_qsat_liq = .FALSE.
     807     CALL getin('adv_qsat_liq',adv_qsat_liq)
     808
    801809     !Config  Key  = ok_dynzon
    802810     !Config  Desc = sortie des transports zonaux dans la dynamique
    … …  
    905913     write(lunout,*)' offline = ', offline
    906914     write(lunout,*)' type_trac = ', type_trac
     915     write(lunout,*)' adv_qsat_liq = ', adv_qsat_liq
    907916     write(lunout,*)' ok_dynzon = ', ok_dynzon
    908917     write(lunout,*)' ok_dyn_ins = ', ok_dyn_ins

  • LMDZ6/branches/contrails/libf/dyn3d/qminimum.f90

    r5285 r5601  
    2222  INTEGER, SAVE :: iq_vap, iq_liq        ! indices pour l'eau vapeur/liquide
    2323  REAL, PARAMETER :: seuil_vap = 1.0e-10 ! seuil pour l'eau vapeur
    24   REAL, PARAMETER :: seuil_liq = 1.0e-11 ! seuil pour l'eau liquide
     24  REAL, PARAMETER :: seuil_liq = 1.0e-15 ! seuil pour l'eau liquide
    2525  !
    2626  !  NB. ....( Il est souhaitable mais non obligatoire que les valeurs des

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

    r5296 r5601  
    5050    ale_bl, ale_bl_trig, alp_bl, &
    5151    ale_wake, ale_bl_stat, AWAKE_S, &
    52     cf_ancien, rvc_ancien
     52    cf_ancien, qvc_ancien, cfa_ancien, qta_ancien
    5353
    5454  USE comconst_mod, ONLY: pi, dtvr
    … …  
    242242
    243243  cf_ancien = 0.
    244   rvc_ancien = 0.
     244  qvc_ancien = 0.
     245  cfa_ancien = 0.
     246  qta_ancien = 0.
    245247 
    246248  z0m(:,:)=0 ! ym missing 5th subsurface initialization

  • LMDZ6/branches/contrails/libf/misc/readTracFiles_mod.f90

    r5536 r5601  
    122122  !--- SOME PARAMETERS THAT ARE NOT LIKELY TO CHANGE OFTEN
    123123  CHARACTER(LEN=maxlen), SAVE      :: tran0        = 'air'      !--- Default transporting fluid
    124   CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlibfca'    !--- Old phases for water (no separator)
    125   CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsbfca'    !--- Known phases initials
     124  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlibfcaq'    !--- Old phases for water (no separator)
     125  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsbfcaq'    !--- Known phases initials
    126126  INTEGER, PARAMETER :: nphases = LEN_TRIM(known_phases)        !--- Number of phases
    127127  CHARACTER(LEN=maxlen), SAVE      :: phases_names(nphases) &   !--- Known phases names
    128                                 = ['gaseous  ', 'liquid   ', 'solid    ','blownSnow', 'fracCloud', 'cldVapRat', 'aviContRa']
     128                                = ['gaseous  ', 'liquid   ', 'solid    ','blownSnow', 'fracCloud', 'cldVapor ', 'aviContra', 'qtContra ']
    129129  CHARACTER(LEN=1),      SAVE :: phases_sep  =  '_'             !--- Phase separator
    130130  CHARACTER(LEN=maxlen), SAVE :: isoFile = 'isotopes_params.def'!--- Name of the isotopes parameters file

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

    r5296 r5601  
    259259    ! LSCP condensation and ice supersaturation
    260260    cf_ancien = 0.
    261     rvc_ancien = 0.
     261    qvc_ancien = 0.
    262262
    263263    wake_delta_pbl_TKE(:,:,:)=0

  • LMDZ6/branches/contrails/libf/phylmd/dyn1d/old_lmdz1d.f90

    r5368 r5601  
    2222       zgam, zmax0, zmea, zpic, zsig, &
    2323       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
    24        ql_ancien, qs_ancien, qbs_ancien, cf_ancien, rvc_ancien, &
     24       ql_ancien, qs_ancien, qbs_ancien, &
     25       cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, &
    2526       prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
    2627       u10m,v10m,ale_wake,ale_bl_stat
    … …  
    872873        IF ( ok_ice_supersat ) THEN
    873874          cf_ancien = 0.
    874           rvc_ancien = 0.
     875          qvc_ancien = 0.
     876        ENDIF
     877        IF ( ok_plane_contrail ) THEN
     878          cfa_ancien = 0.
     879          qta_ancien = 0.
    875880        ENDIF
    876881!jyg<

  • LMDZ6/branches/contrails/libf/phylmd/dyn1d/scm.f90

    r5450 r5601  
    1818       zgam, zmax0, zmea, zpic, zsig, &
    1919       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
    20        ql_ancien, qs_ancien, qbs_ancien, cf_ancien, rvc_ancien, &
     20       ql_ancien, qs_ancien, qbs_ancien, &
     21       cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, &
    2122       prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
    2223       u10m,v10m,ale_wake,ale_bl_stat, ratqs_inter_
    … …  
    614615        IF ( ok_ice_supersat ) THEN
    615616          cf_ancien = 0.
    616           rvc_ancien = 0.
     617          qvc_ancien = 0.
     618        ENDIF
     619        IF ( ok_plane_contrail ) THEN
     620          cfa_ancien = 0.
     621          qta_ancien = 0.
    617622        ENDIF
    618623        rain_fall=0.

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

    r5589 r5601  
    7171
    7272!**********************************************************************************
    73 SUBROUTINE contrails_formation_evolution( &
    74       dtime, pplay, temp, qsat, qsatl, gamma_cond, rcont_seri, flight_dist, &
    75       cldfra, qvc, pdf_loc, pdf_scale, pdf_alpha, &
    76       Tcritcont, qcritcont, potcontfraP, potcontfraNP, contfra, &
    77       dcontfra_cir, dcf_avi, dqvc_avi, dqi_avi &
    78       )
    79 
    80 USE lmdz_lscp_ini, ONLY: RCPD, EPS_W, RTT
     73SUBROUTINE contrails_formation( &
     74      klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, flight_dist, &
     75      cldfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &
     76      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
     77      dcfa_ini, dqia_ini, dqta_ini)
     78
     79USE lmdz_lscp_ini, ONLY: RPI, RCPD, EPS_W, RTT
    8180USE lmdz_lscp_ini, ONLY: EI_H2O_aviation, qheat_fuel_aviation, prop_efficiency_aviation
    82 USE lmdz_lscp_ini, ONLY: linear_contrails_lifetime
    83 USE lmdz_lscp_ini, ONLY: eps
     81USE lmdz_lscp_ini, ONLY: eps, temp_nowater
    8482
    8583USE lmdz_lscp_tools, ONLY: GAMMAINC, calc_qsat_ecmwf
    … …  
    9088! Input
    9189!
    92 REAL, INTENT(IN)  :: dtime        ! time step [s]
    93 REAL, INTENT(IN)  :: pplay        ! layer pressure [Pa]
    94 REAL, INTENT(IN)  :: temp         ! temperature [K]
    95 REAL, INTENT(IN)  :: qsat         ! saturation specific humidity [kg/kg]
    96 REAL, INTENT(IN)  :: qsatl        ! saturation specific humidity w.r.t. liquid [kg/kg]
    97 REAL, INTENT(IN)  :: gamma_cond   ! condensation threshold w.r.t. qsat [-]
    98 REAL, INTENT(IN)  :: rcont_seri   ! ratio of contrails fraction to total cloud fraction [-]
    99 REAL, INTENT(IN) :: flight_dist  ! aviation distance flown concentration [m/s/m3]
    100 REAL, INTENT(IN) :: cldfra       ! cloud fraction [-]
    101 REAL, INTENT(IN)  :: qvc          ! gridbox-mean vapor in the cloud [kg/kg]
    102 REAL, INTENT(IN)  :: pdf_loc      ! location parameter of the clear sky PDF [%]
    103 REAL, INTENT(IN)  :: pdf_scale    ! scale parameter of the clear sky PDF [%]
    104 REAL, INTENT(IN)  :: pdf_alpha    ! shape parameter of the clear sky PDF [-]
     90INTEGER,  INTENT(IN)                   :: klon         ! number of horizontal grid points
     91REAL,     INTENT(IN)                   :: dtime        ! time step [s]
     92REAL,     INTENT(IN) , DIMENSION(klon) :: pplay        ! layer pressure [Pa]
     93REAL,     INTENT(IN) , DIMENSION(klon) :: temp         ! temperature [K]
     94REAL,     INTENT(IN) , DIMENSION(klon) :: qsat         ! saturation specific humidity [kg/kg]
     95REAL,     INTENT(IN) , DIMENSION(klon) :: qsatl        ! saturation specific humidity w.r.t. liquid [kg/kg]
     96REAL,     INTENT(IN) , DIMENSION(klon) :: gamma_cond   ! condensation threshold w.r.t. qsat [-]
     97REAL,     INTENT(IN) , DIMENSION(klon) :: flight_dist  ! aviation distance flown concentration [m/s/m3]
     98REAL,     INTENT(IN) , DIMENSION(klon) :: cldfra       ! cloud fraction [-]
     99REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_loc      ! location parameter of the clear sky PDF [%]
     100REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_scale    ! scale parameter of the clear sky PDF [%]
     101REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_alpha    ! shape parameter of the clear sky PDF [-]
     102LOGICAL,  INTENT(IN) , DIMENSION(klon) :: keepgoing    ! .TRUE. if a new condensation loop should be computed
    105103!
    106104! Output
    107105!
    108 REAL, INTENT(OUT) :: Tcritcont    ! critical temperature for contrail formation [K]
    109 REAL, INTENT(OUT) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
    110 REAL, INTENT(OUT) :: potcontfraP  ! potential persistent contrail fraction [-]
    111 REAL, INTENT(OUT) :: potcontfraNP ! potential non-persistent contrail fraction [-]
    112 REAL, INTENT(OUT) :: contfra      ! linear contrail fraction [-]
    113 REAL, INTENT(OUT) :: dcontfra_cir ! linear contrail fraction to cirrus cloud fraction tendency [s-1]
    114 REAL, INTENT(OUT) :: dcf_avi      ! cloud fraction tendency because of aviation [s-1]
    115 REAL, INTENT(OUT) :: dqvc_avi     ! specific ice content tendency because of aviation [kg/kg/s]
    116 REAL, INTENT(OUT) :: dqi_avi      ! specific cloud water vapor tendency because of aviation [kg/kg/s]
     106REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
     107REAL,     INTENT(INOUT), DIMENSION(klon) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
     108REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
     109REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-]
     110REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfa_ini     ! contrails cloud fraction tendency bec ause of initial formation [s-1]
     111REAL,     INTENT(INOUT), DIMENSION(klon) :: dqia_ini     ! contrails ice specific humidity tende ncy because of initial formation [kg/kg/s]
     112REAL,     INTENT(INOUT), DIMENSION(klon) :: dqta_ini     ! contrails total specific humidity ten dency because of initial formation [kg/kg/s]
    117113!
    118114! Local
    119115!
     116INTEGER :: i
    120117! for Schmidt-Appleman criteria
    121 REAL, DIMENSION(1) :: ztemp, zpplay, qzero, zqsatl, zdqsatl
    122 REAL :: Gcont, qsatl_crit, psatl_crit, pcrit
     118REAL, DIMENSION(klon) :: qzero
     119REAL, DIMENSION(klon) :: qsatl_crit, dqsatl_crit
     120REAL :: Gcont, psatl_crit, pcrit
    123121REAL :: pdf_x, pdf_y, pdf_e2, pdf_e3
    124122REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat, pdf_fra_above_qnuc
    … …  
    131129qzero(:) = 0.
    132130
     131DO i = 1, klon
     132  IF ( keepgoing(i) ) THEN
     133    Tcritcont(i)    = 0.
     134    qcritcont(i)    = 0.
     135    potcontfraP(i)  = 0.
     136    potcontfraNP(i) = 0.
     137  ENDIF
     138ENDDO
     139
    133140!---------------------------------
    134 !--  SCHIMDT-APPLEMAN CRITERIA  --
     141!--  SCHMIDT-APPLEMAN CRITERIA  --
    135142!---------------------------------
    136143!--Revised by Schumann (1996) and Rap et al. (2010)
    137144
    138 !--Gcont is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
    139 !--in Pa / K. See Rap et al. (2010) in JGR.
    140 !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
    141 Gcont = EI_H2O_aviation * RCPD * pplay &
    142        / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
    143 !--critical temperature below which no liquid contrail can form in exhaust
    144 !--noted T_LM in Schumann (1996), their Eq. 31 in Appendix 2
    145 !--in Kelvins
    146 Tcritcont = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) &
    147        + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2
    148 
    149 ztemp(1) = Tcritcont
    150 zpplay(1) = pplay
    151 CALL calc_qsat_ecmwf(1, ztemp, qzero, zpplay, RTT, 1, .FALSE., zqsatl, zdqsatl)
    152 qsatl_crit = zqsatl(1)
    153 
    154 psatl_crit = qsatl_crit / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit ) * pplay
    155 pcrit = Gcont * ( temp - Tcritcont ) + psatl_crit
    156 qcritcont = EPS_W * pcrit / ( pplay - ( 1. - EPS_W ) * pcrit )
    157 
    158 
    159 IF ( ( ( 1. - cldfra ) .GT. eps ) .AND. ( temp .LT. Tcritcont ) ) THEN
    160 
    161   pdf_x = qcritcont / qsatl * 100.
    162   pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
    163   pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
    164   pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
    165   pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
    166   pdf_fra_above_qcritcont = EXP( - pdf_y ) * ( 1. - cldfra )
    167   pdf_q_above_qcritcont = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qcritcont ) * qsatl / 100.
    168 
    169   pdf_x = gamma_cond * qsat / qsatl * 100.
    170   pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
    171   pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
    172   pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
    173   pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
    174   pdf_fra_above_qnuc = EXP( - pdf_y ) * ( 1. - cldfra )
    175   pdf_q_above_qnuc = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qnuc ) * qsatl / 100.
    176 
    177   pdf_x = qsat / qsatl * 100.
    178   pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
    179   pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
    180   pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
    181   pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
    182   pdf_fra_above_qsat = EXP( - pdf_y ) * ( 1. - cldfra )
    183   pdf_q_above_qsat = ( pdf_e3 * ( 1. - cldfra ) + pdf_loc * pdf_fra_above_qsat ) * qsatl / 100.
    184 
    185   potcontfraNP = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
    186   potcontfraP = MAX(0., MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
    187                     pdf_fra_above_qcritcont - pdf_fra_above_qnuc))
    188   qpotcontP = MAX(0., MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
    189                   pdf_q_above_qcritcont - pdf_q_above_qnuc))
    190 
    191 ELSE
    192 
    193   potcontfraNP = 0.
    194   potcontfraP = 0.
    195 
    196 ENDIF ! temp .LT. Tcritcont
    197 
    198 !--Convert existing contrail fraction into "natural" cirrus cloud fraction
    199 contfra = rcont_seri * cldfra
    200 dcontfra_cir = contfra * ( EXP( - dtime / linear_contrails_lifetime ) - 1. )
    201 contfra = contfra + dcontfra_cir
    202 
    203 !--Add a source of contrails from aviation
    204 dcf_avi = 0.
    205 dqi_avi = 0.
    206 dqvc_avi = 0.
    207 IF ( potcontfraP .GT. eps ) THEN
    208   contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
    209   contfra_new = MIN(1., flight_dist * dtime * contrail_cross_section)
    210   dcf_avi = potcontfraP * contfra_new
    211   IF ( cldfra .GT. eps ) THEN
    212     dqvc_avi = dcf_avi * qvc / cldfra
    213   ELSE
    214     dqvc_avi = dcf_avi * qsat
    215   ENDIF
    216   dqi_avi = dcf_avi * qpotcontP / potcontfraP - dqvc_avi
    217 
    218   !--Add tendencies
    219   contfra = contfra + dcf_avi
    220 ENDIF
    221 
    222 END SUBROUTINE contrails_formation_evolution
     145DO i = 1, klon
     146  !--Gcont is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute
     147  !--in Pa / K. See Rap et al. (2010) in JGR.
     148  !--kg H2O/kg fuel * J kg air-1 K-1 * Pa / (kg H2O / kg air * J kg fuel-1)
     149  Gcont = EI_H2O_aviation * RCPD * pplay(i) &
     150         / ( EPS_W * qheat_fuel_aviation * ( 1. - prop_efficiency_aviation ) )
     151  !--critical temperature below which no liquid contrail can form in exhaust
     152  !--noted T_LM in Schumann (1996), their Eq. 31 in Appendix 2
     153  !--in Kelvins
     154  Tcritcont(i) = 226.69 + 9.43 * LOG( MAX(Gcont, 0.1) - 0.053 ) &
     155         + 0.72 * LOG( MAX(Gcont, 0.1) - 0.053 )**2
     156ENDDO
     157
     158CALL calc_qsat_ecmwf(klon, Tcritcont, qzero, pplay, RTT, 1, .FALSE., qsatl_crit, dqsatl_crit)
     159
     160DO i = 1, klon
     161  IF ( keepgoing(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
     162   
     163    psatl_crit = qsatl_crit(i) / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit(i) ) * pplay(i)
     164    pcrit = Gcont * ( temp(i) - Tcritcont(i) ) + psatl_crit
     165    qcritcont(i) = EPS_W * pcrit / ( pplay(i) - ( 1. - EPS_W ) * pcrit )
     166
     167
     168    IF ( ( ( 1. - cldfra(i) ) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN
     169   
     170      pdf_x = qcritcont(i) / qsatl(i) * 100.
     171      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     172      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
     173      pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     174      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     175      pdf_fra_above_qcritcont = EXP( - pdf_y ) * ( 1. - cldfra(i) )
     176      pdf_q_above_qcritcont = ( pdf_e3 * ( 1. - cldfra(i) ) &
     177          + pdf_loc(i) * pdf_fra_above_qcritcont ) * qsatl(i) / 100.
     178   
     179      pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
     180      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     181      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
     182      pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     183      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     184      pdf_fra_above_qnuc = EXP( - pdf_y ) * ( 1. - cldfra(i) )
     185      pdf_q_above_qnuc = ( pdf_e3 * ( 1. - cldfra(i) ) &
     186          + pdf_loc(i) * pdf_fra_above_qnuc ) * qsatl(i) / 100.
     187   
     188      pdf_x = qsat(i) / qsatl(i) * 100.
     189      pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     190      pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
     191      pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     192      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     193      pdf_fra_above_qsat = EXP( - pdf_y ) * ( 1. - cldfra(i) )
     194      pdf_q_above_qsat = ( pdf_e3 * ( 1. - cldfra(i) ) &
     195          + pdf_loc(i) * pdf_fra_above_qsat ) * qsatl(i) / 100.
     196   
     197      potcontfraNP(i) = MAX(0., pdf_fra_above_qcritcont - pdf_fra_above_qsat)
     198      potcontfraP(i) = MAX(0., MIN(pdf_fra_above_qsat - pdf_fra_above_qnuc, &
     199                        pdf_fra_above_qcritcont - pdf_fra_above_qnuc))
     200      qpotcontP = MAX(0., MIN(pdf_q_above_qsat - pdf_q_above_qnuc, &
     201                      pdf_q_above_qcritcont - pdf_q_above_qnuc))
     202
     203    ENDIF ! temp .LT. Tcritcont
     204   
     205    !--Add a source of contrails from aviation
     206    IF ( potcontfraP(i) .GT. eps ) THEN
     207      contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
     208      contfra_new = MIN(1., flight_dist(i) * dtime * contrail_cross_section)
     209      dcfa_ini(i) = potcontfraP(i) * contfra_new
     210      dqta_ini(i) = qpotcontP * contfra_new
     211      dqia_ini(i) = dqta_ini(i) - qsat(i) * dcfa_ini(i)
     212    ENDIF
     213
     214  ENDIF ! keepgoing
     215ENDDO
     216
     217END SUBROUTINE contrails_formation
    223218
    224219
    225220!**********************************************************************************
    226 SUBROUTINE contrails_mixing( &
    227       dtime, pplay, shear, pbl_eps, cell_area, dz, temp, qtot, qsat, &
    228       subfra, qsub, issrfra, qissr, cldfra, qcld, qvc, rcont_seri, &
    229       dcf_mix, dqvc_mix, dqi_mix, dqt_mix, dcf_mix_issr, dqt_mix_issr &
    230       )
    231 
    232 !----------------------------------------------------------------------
    233 ! This subroutine calculates the mixing of contrails in their environment.
    234 ! Authors: Audran Borella, Etienne Vignon, Olivier Boucher
    235 ! December 2024
    236 !----------------------------------------------------------------------
    237 
    238 USE lmdz_lscp_ini,        ONLY: RPI, eps, ok_unadjusted_clouds
    239 USE lmdz_lscp_ini,        ONLY: aspect_ratio_contrails, coef_mixing_contrails, &
    240                                 coef_shear_contrails, chi_mixing_contrails
    241 
    242 IMPLICIT NONE
    243 
    244 !
    245 ! Input
    246 !
    247 REAL, INTENT(IN)    :: dtime        ! time step [s]
    248 !
    249 REAL, INTENT(IN)    :: pplay        ! layer pressure [Pa]
    250 REAL, INTENT(IN)    :: shear        ! vertical shear [s-1]
    251 REAL, INTENT(IN)    :: pbl_eps      ! TKE dissipation[m2/s3]
    252 REAL, INTENT(IN)    :: cell_area    ! cell area [m2]
    253 REAL, INTENT(IN)    :: dz           ! cell width [m]
    254 REAL, INTENT(IN)    :: temp         ! temperature [K]
    255 REAL, INTENT(IN)    :: qtot         ! total specific humidity (without precip) [kg/kg]
    256 REAL, INTENT(IN)    :: qsat         ! saturation specific humidity [kg/kg]
    257 REAL, INTENT(IN)    :: subfra       ! subsaturated clear sky fraction [-]
    258 REAL, INTENT(IN)    :: qsub         ! gridbox-mean subsaturated clear sky specific water [kg/kg]
    259 REAL, INTENT(IN)    :: issrfra      ! ISSR fraction [-]
    260 REAL, INTENT(IN)    :: qissr        ! gridbox-mean ISSR specific water [kg/kg]
    261 REAL, INTENT(IN)    :: cldfra       ! cloud fraction [-]
    262 REAL, INTENT(IN)    :: qcld         ! gridbox-mean cloudy specific total water [kg/kg]
    263 REAL, INTENT(IN)    :: qvc          ! gridbox-mean cloudy specific water vapor [kg/kg]
    264 REAL, INTENT(IN)    :: rcont_seri   ! ratio of contrails fraction to total cloud fraction [-]
    265 !
    266 !  Input/Output
    267 !
    268 REAL, INTENT(INOUT) :: dcf_mix      ! cloud fraction tendency because of cloud mixing [s-1]
    269 REAL, INTENT(INOUT) :: dqvc_mix     ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
    270 REAL, INTENT(INOUT) :: dqi_mix      ! specific ice content tendency because of cloud mixing [kg/kg/s]
    271 REAL, INTENT(INOUT) :: dqt_mix      ! total water tendency because of cloud mixing [kg/kg/s]
    272 REAL, INTENT(INOUT) :: dcf_mix_issr ! cloud fraction tendency because of cloud mixing in ISSR [s-1]
    273 REAL, INTENT(INOUT) :: dqt_mix_issr ! total water tendency because of cloud mixing in ISSR [kg/kg/s]
    274 !
    275 ! Local
    276 !
    277 REAL :: dqt_mix_sub_cont, dqt_mix_issr_cont
    278 REAL :: dcf_mix_sub_cont, dcf_mix_issr_cont
    279 REAL :: dqvc_mix_sub_cont, dqvc_mix_issr_cont
    280 REAL :: dcf_mix_cont, dqvc_mix_cont, dqi_mix_cont, dqt_mix_cont
    281 REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
    282 REAL :: envfra_mix, cldfra_mix
    283 REAL :: L_shear, shear_fra
    284 REAL :: sigma_mix, issrfra_mix, subfra_mix
    285 REAL :: qvapincld, qvapinmix, qvapincld_new, qiceincld
    286 
    287 !--Initialisation
    288 dcf_mix_sub_cont   = 0.
    289 dqt_mix_sub_cont   = 0.
    290 dqvc_mix_sub_cont  = 0.
    291 dcf_mix_issr_cont  = 0.
    292 dqt_mix_issr_cont  = 0.
    293 dqvc_mix_issr_cont = 0.
    294 
    295 !-- PART 1 - TURBULENT DIFFUSION
    296 
    297 !--Clouds within the mesh are assumed to be ellipses. The length of the
    298 !--semi-major axis is a and the length of the semi-minor axis is b.
    299 !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
    300 !--In particular, it is 0 if cldfra = 1.
    301 !--clouds_perim is the total perimeter of the clouds within the mesh,
    302 !--not considering interfaces with other meshes (only the interfaces with clear
    303 !--sky are taken into account).
    304 !--
    305 !--The area of each cloud is A = a * b * RPI,
    306 !--and the perimeter of each cloud is
    307 !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
    308 !--
    309 !--With cell_area the area of the cell, we have:
    310 !-- cldfra = A * N_cld_mix / cell_area
    311 !-- clouds_perim = P * N_cld_mix
    312 !--
    313 !--We assume that the ratio between b and a is a function of
    314 !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
    315 !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
    316 !--are spherical.
    317 !-- b / a = bovera = MAX(0.1, cldfra)
    318 bovera = aspect_ratio_contrails
    319 !--P / a is a function of b / a only, that we can calculate
    320 !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
    321 Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
    322 !--The clouds perimeter is imposed using the formula from Morcrette 2012,
    323 !--based on observations.
    324 !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
    325 !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
    326 !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
    327 !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
    328 a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - cldfra ) / bovera
    329 !--and finally,
    330 !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
    331 N_cld_mix = coef_mixing_contrails * cldfra * ( 1. - cldfra ) * cell_area &
    332           / Povera / a_mix
    333 
    334 !--The time required for turbulent diffusion to homogenize a region of size
    335 !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
    336 !--We compute L_mix and assume that the cloud is mixed over this length
    337 L_mix = SQRT( dtime**3 * pbl_eps )
    338 !--The mixing length cannot be greater than the semi-minor axis. In this case,
    339 !--the entire cloud is mixed.
    340 L_mix = MIN(L_mix, a_mix * bovera)
    341 
    342 !--The fraction of clear sky mixed is
    343 !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
    344 envfra_mix = N_cld_mix * RPI / cell_area &
    345            * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
    346 !--The fraction of cloudy sky mixed is
    347 !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
    348 cldfra_mix = N_cld_mix * RPI / cell_area &
    349            * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
    350 
    351 
    352 !-- PART 2 - SHEARING
    353 
    354 !--The clouds are then sheared. We keep the shape and number
    355 !--assumptions from before. The clouds are sheared with a random orientation
    356 !--of the wind, on average we assume that the wind and the semi-major axis
    357 !--make a 45 degrees angle. Moreover, the contrails only mix
    358 !--along their semi-minor axis (b), because it is easier to compute.
    359 !--With this, the clouds increase in size along b only, by a factor
    360 !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind)
    361 L_shear = coef_shear_contrails * shear * dz * dtime
    362 !--therefore, the fraction of clear sky mixed is
    363 !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area
    364 !--and the fraction of cloud mixed is
    365 !-- N_cld_mix * ( a * b - a * (b - L_shear * SQRT(2.) / 2.) ) * RPI / 2. / cell_area
    366 !--(note that they are equal)
    367 shear_fra = RPI * SQRT(2.) / 2. * L_shear * a_mix / 2. * N_cld_mix / cell_area
    368 !--and the environment and cloud mixed fractions are the same,
    369 !--which we add to the previous calculated mixed fractions.
    370 !--We therefore assume that the sheared clouds and the turbulent
    371 !--mixed clouds are different.
    372 envfra_mix = envfra_mix + shear_fra
    373 cldfra_mix = cldfra_mix + shear_fra
    374 
    375 
    376 !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
    377 
    378 !--The environment fraction is allocated to subsaturated sky or supersaturated sky,
    379 !--according to the factor sigma_mix. This is computed as the ratio of the
    380 !--subsaturated sky fraction to the environment fraction, corrected by a factor
    381 !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
    382 !--supersaturated sky is favoured. Physically, this means that it is more likely
    383 !--to have supersaturated sky around the cloud than subsaturated sky.
    384 sigma_mix = subfra / ( subfra + chi_mixing_contrails * issrfra )
    385 subfra_mix = MIN( sigma_mix * envfra_mix, subfra )
    386 issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra )
    387 cldfra_mix = MIN( cldfra_mix, cldfra )
    388 
    389 !--First, we mix the subsaturated sky (subfra_mix) and the cloud close
    390 !--to this fraction (sigma_mix * cldfra_mix).
    391 IF ( subfra .GT. eps ) THEN
    392   !--We compute the total humidity in the mixed air, which
    393   !--can be either sub- or supersaturated.
    394   qvapinmix = ( qsub * subfra_mix / subfra &
    395             + qcld * cldfra_mix * sigma_mix / cldfra ) &
    396             / ( subfra_mix + cldfra_mix * sigma_mix )
    397 
    398   IF ( ok_unadjusted_clouds ) THEN
    399     qiceincld = ( qcld - qvc ) * cldfra_mix * sigma_mix / cldfra &
    400               / ( subfra_mix + cldfra_mix * sigma_mix )
    401     qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( &
    402         qvapinmix, qiceincld, temp, qsat, pplay, dtime)
    403     IF ( qvapincld_new .EQ. 0. ) THEN
    404       !--If all the ice has been sublimated, we sublimate
    405       !--completely the cloud and do not activate the sublimation
    406       !--process
    407       !--Tendencies and diagnostics
    408       dcf_mix_sub_cont = - sigma_mix * cldfra_mix
    409       dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra
    410       dqvc_mix_sub_cont = dcf_mix_sub_cont * qvc / cldfra
    411     ELSE
    412       dcf_mix_sub_cont = subfra_mix
    413       dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra
    414       dqvc_mix_sub_cont = dcf_mix_sub_cont * qvapincld_new
    415     ENDIF
    416   ELSE
    417     IF ( qvapinmix .GT. qsat ) THEN
    418       !--If the mixed air is supersaturated, we condense the subsaturated
    419       !--region which was mixed.
    420       dcf_mix_sub_cont = subfra_mix
    421       dqt_mix_sub_cont = dcf_mix_sub_cont * qsub / subfra
    422       dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat
    423     ELSE
    424       !--Else, we sublimate the cloud which was mixed.
    425       dcf_mix_sub_cont = - sigma_mix * cldfra_mix
    426       dqt_mix_sub_cont = dcf_mix_sub_cont * qcld / cldfra
    427       dqvc_mix_sub_cont = dcf_mix_sub_cont * qsat
    428     ENDIF
    429   ENDIF ! ok_unadjusted_clouds
    430 ENDIF ! subfra .GT. eps
    431 
    432 !--We then mix the supersaturated sky (issrfra_mix) and the cloud,
    433 !--for which the mixed air is always supersatured, therefore
    434 !--the cloud necessarily expands
    435 IF ( issrfra .GT. eps ) THEN
    436 
    437   IF ( ok_unadjusted_clouds ) THEN
    438     qvapinmix = ( qissr * issrfra_mix / issrfra &
    439               + qcld * cldfra_mix * (1. - sigma_mix) / cldfra ) &
    440               / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
    441     qiceincld = ( qcld - qvc ) * cldfra_mix * (1. - sigma_mix) / cldfra &
    442               / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
    443     qvapincld_new = QVAPINCLD_DEPSUB_CONTRAILS( &
    444         qvapinmix, qiceincld, temp, qsat, pplay, dtime)
    445     dcf_mix_issr_cont = issrfra_mix
    446     dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra
    447     dqvc_mix_issr_cont = dcf_mix_issr_cont * qvapincld_new
    448   ELSE
    449     !--In this case, the additionnal vapor condenses
    450     dcf_mix_issr_cont = issrfra_mix
    451     dqt_mix_issr_cont = dcf_mix_issr_cont * qissr / issrfra
    452     dqvc_mix_issr_cont = dcf_mix_issr_cont * qsat
    453   ENDIF ! ok_unadjusted_clouds
    454 
    455 ENDIF ! issrfra .GT. eps
    456 
    457 !--Sum up the tendencies from subsaturated sky and supersaturated sky
    458 dcf_mix_cont  = dcf_mix_sub_cont  + dcf_mix_issr_cont
    459 dqt_mix_cont  = dqt_mix_sub_cont  + dqt_mix_issr_cont
    460 dqvc_mix_cont = dqvc_mix_sub_cont + dqvc_mix_issr_cont
    461 dqi_mix_cont  = dqt_mix_cont - dqvc_mix_cont
    462 
    463 !--Outputs
    464 !--The mixing tendencies are a linear combination of the calculation done for "classical" cirrus
    465 !--and contrails
    466 dcf_mix = ( 1. - rcont_seri ) * dcf_mix + rcont_seri * dcf_mix_cont
    467 dqvc_mix = ( 1. - rcont_seri ) * dqvc_mix + rcont_seri * dqvc_mix_cont
    468 dqi_mix = ( 1. - rcont_seri ) * dqi_mix + rcont_seri * dqi_mix_cont
    469 dqt_mix = ( 1. - rcont_seri ) * dqt_mix + rcont_seri * dqt_mix_cont
    470 dcf_mix_issr = ( 1. - rcont_seri ) * dcf_mix_issr + rcont_seri * dcf_mix_issr_cont
    471 dqt_mix_issr = ( 1. - rcont_seri ) * dqt_mix_issr + rcont_seri * dqt_mix_issr_cont
    472 
    473 END SUBROUTINE contrails_mixing
    474 
    475 
    476 !**********************************************************************************
    477 FUNCTION qvapincld_depsub_contrails( &
     221FUNCTION QVAPINCLD_DEPSUB_CONTRAILS( &
    478222    qvapincld, qiceincld, temp, qsat, pplay, dtime)
    479223
    … …  
    494238REAL :: qvapincld_depsub_contrails
    495239! local
    496 REAL :: qice_ratio
    497240REAL :: pres_sat, rho, kappa
    498241REAL :: air_thermal_conduct, water_vapor_diff
    … …  
    583326ELSE
    584327  !--If the cloud is initially subsaturated
    585   !--Exact explicit integration (qice exact, qvc explicit)
    586   !--The barrier is set so that the resulting vapor in cloud
    587   !--cannot be greater than qsat
    588   !--qice_ratio is the ratio between the new ice content and
    589   !--the old one, it is comprised between 0 and 1
    590   qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )
    591 
    592   IF ( qice_ratio .LT. 0. ) THEN
    593     !--The new vapor in cloud is set to 0 so that the
    594     !--sublimation process does not activate
    595     qvapincld_depsub_contrails = 0.
    596   ELSE
    597     !--Else, the sublimation process is activated with the
    598     !--diagnosed new cloud water vapor
    599     !--The new vapor in the cloud is increased with the
    600     !--sublimated ice
    601     qvapincld_depsub_contrails = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
    602     !--The new vapor in the cloud cannot be greater than qsat
    603     qvapincld_depsub_contrails = MIN(qvapincld_depsub_contrails, qsat)
    604   ENDIF ! qice_ratio .LT. 0.
     328  !--Exact explicit integration (qvc exact, qice explicit)
     329  !--Same but depo_coef_cirrus = 1
     330  qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
     331                * EXP( - dtime * qiceincld / kappa / rm_ice**2 )
    605332ENDIF ! qvapincld .GT. qsat
    606333
    607 END FUNCTION qvapincld_depsub_contrails
     334END FUNCTION QVAPINCLD_DEPSUB_CONTRAILS
    608335
    609336

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

    r5598 r5601  
    1212    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
    1313    !--AB contrails
    14     contfra, radocond_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
    15     pch_nocont, pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
     14    contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
     15    pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     16    xfiwp_nocont, xfiwc_nocont, reice_nocont)
    1617
    1718  ! Interface between the LMDZ physics monitor and the cloud properties calculation routines
    … …  
    9596  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
    9697  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
    97   REAL, INTENT(IN)  :: radocond_cont(klon, klev) ! contrails condensed water seen by radiation [kg/kg]
     98  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
    9899  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
    99100  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
    … …  
    101102  REAL, INTENT(OUT) :: xfiwc_nocont(klon, klev)  ! ice water content seen by radiation without contrails [kg/kg]
    102103  REAL, INTENT(OUT) :: pclc_nocont(klon, klev)   ! cloud fraction for radiation without contrails [-]
    103   REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [m]
     104  REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [-]
    104105  REAL, INTENT(OUT) :: pclemi_nocont(klon, klev) ! cloud emissivity without contrails [-]
    105   REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [m]
     106  REAL, INTENT(OUT) :: pcltau_cont(klon, klev)   ! contrails optical depth [-]
     107  REAL, INTENT(OUT) :: pclemi_cont(klon, klev)   ! contrails emissivity [-]
     108  REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [micronts]
    106109  !--AB
    107110
    … …  
    134137    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
    135138    !--AB for contrails
    136     contfra, radocond_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
    137     pch_nocont, pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
     139    contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
     140    pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     141    xfiwp_nocont, xfiwc_nocont, reice_nocont)
    138142  ELSE
    139143    CALL nuage (paprs, pplay, &

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

    r5598 r5601  
    1111    icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, &
    1212    !--AB contrails
    13     contfra, radocond_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
    14     pch_nocont, pct_cont, xfiwp_nocont, xfiwc_nocont, reice_nocont)
     13    contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
     14    pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     15    xfiwp_nocont, xfiwc_nocont, reice_nocont)
    1516
    1617  USE lmdz_cloud_optics_prop_ini , ONLY : flag_aerosol, ok_cdnc
    … …  
    119120  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
    120121  REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
    121   REAL, INTENT(IN)  :: radocond_cont(klon, klev) ! contrails condensed water seen by radiation [kg/kg]
     122  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
    122123  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
    123124  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
    … …  
    125126  REAL, INTENT(OUT) :: xfiwc_nocont(klon, klev)  ! ice water content seen by radiation without contrails [kg/kg]
    126127  REAL, INTENT(OUT) :: pclc_nocont(klon, klev)   ! cloud fraction for radiation without contrails [-]
    127   REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [m]
     128  REAL, INTENT(OUT) :: pcltau_nocont(klon, klev) ! cloud optical depth without contrails [-]
    128129  REAL, INTENT(OUT) :: pclemi_nocont(klon, klev) ! cloud emissivity without contrails [-]
    129   REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [m]
     130  REAL, INTENT(OUT) :: pcltau_cont(klon, klev)   ! contrails optical depth [-]
     131  REAL, INTENT(OUT) :: pclemi_cont(klon, klev)   ! contrails emissivity [-]
     132  REAL, INTENT(OUT) :: reice_nocont(klon, klev)  ! ice effective radius without contrails [microns]
    130133  !--AB
    131134
    … …  
    168171  REAL zflwp_var, zfiwp_var
    169172  REAL d_rei_dt
     173  REAL pclc_var
    170174
    171175
    … …  
    249253        DO i = 1, klon
    250254          pclc_nocont(i,k) = pclc(i, k) - contfra(i, k)
    251           xfiwc_nocont(i, k) = icefrac_optics(i, k) * radocond(i, k) - radocond_cont(i, k)
     255          xfiwc_nocont(i, k) = xfiwc(i, k) - qice_cont(i, k)
    252256        ENDDO
    253257      ENDDO
    … …  
    427431  ! --in the general case
    428432
     433  IF ( .NOT. ok_plane_contrail ) THEN
     434
    429435  DO k = 1, klev
    430436    DO i = 1, klon
    … …  
    447453        pcltau(i, k) = 0.0
    448454        pclemi(i, k) = 0.0
    449 
    450         !--AB
    451         IF ( ok_plane_contrail ) THEN
    452           reice_nocont(i,k) = 0.
    453           pclc_nocont(i,k) = 0.
    454           pcltau_nocont(i,k) = 0.
    455           pclemi_nocont(i,k) = 0.
    456         ENDIF
    457455
    458456      ELSE
    … …  
    534532        IF (zflwp_var==0.) rel = 1.
    535533        IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
    536         IF ( ok_plane_contrail ) THEN
    537           !--Diagnostics of clouds emissivity, optical depth and ice crystal radius
    538           !--without contrails
    539           IF ( pclc_nocont(i,k) .LE. seuil_neb ) THEN
    540             reice_nocont(i,k) = 0.
    541             pclc_nocont(i,k) = 0.
    542             pclc(i,k) = contfra(i,k)
    543             pcltau_nocont(i,k) = 0.
    544             pclemi_nocont(i, k) = 0.
    545           ELSE
    546             reice_nocont(i,k) = rei
    547             pcltau_nocont(i,k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/rei)
    548             k_ice = k_ice0 + 1.0/rei
    549             pclemi_nocont(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
    550           ENDIF
    551 
    552           !--If contrails are activated, rei is a weighted average between the natural
    553           !--rei and the contrails rei, with the weights being the fraction of natural
    554           !--vs contrail cirrus in the gridbox
    555           !--Beware, re_ice_crystals_contrails is in m, while rei is in microns
    556           rei = ( rei * pclc_nocont(i,k) &
    557               + re_ice_crystals_contrails * 1.E6 * contfra(i,k) ) / pclc(i,k)
    558         ENDIF
    559534
    560535        pcltau(i, k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/ &
    … …  
    576551      xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
    577552      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
    578       !--AB
    579       IF ( ok_plane_contrail ) THEN
    580           xfiwp_nocont(i) = xfiwp_nocont(i) + xfiwc_nocont(i, k)*rhodz(i, k)
     553
     554    ENDDO
     555  ENDDO
     556
     557  ELSE !--AB if contrails are activated
     558
     559  DO k = 1, klev
     560    DO i = 1, klon
     561
     562      IF (pclc(i,k)<=seuil_neb) THEN
     563
     564        ! effective cloud droplet radius (microns) for liquid water clouds:
     565        ! For output diagnostics cloud droplet effective radius [um]
     566        ! we multiply here with f * xl (fraction of liquid water
     567        ! clouds in the grid cell) to avoid problems in the averaging of the
     568        ! output.
     569        ! In the output of IOIPSL, derive the REAL cloud droplet
     570        ! effective radius as re/fl
     571
     572        fl(i, k) = seuil_neb*(1.-icefrac_optics(i,k))
     573        re(i, k) = rad_chaud(i, k)*fl(i, k)
     574        rel = 0.
     575        rei = 0.
     576        pclc(i, k) = 0.0
     577        pcltau(i, k) = 0.0
     578        pclemi(i, k) = 0.0
     579
     580        !--AB contrails
     581        reice_nocont(i,k) = 0.
     582        pclc_nocont(i,k) = 0.
     583        pcltau_cont(i,k) = 0.
     584        pclemi_cont(i,k) = 0.
     585        pcltau_nocont(i,k) = 0.
     586        pclemi_nocont(i,k) = 0.
     587
     588      ELSE
     589
     590        IF ( pclc_nocont(i,k) .GT. 0.01 * seuil_neb ) THEN
     591          ! -- liquid/ice cloud water paths:
     592
     593          zflwp_var = 1000.*xflwc(i, k)/pclc(i, k)*rhodz(i, k)
     594          zfiwp_var = 1000.*xfiwc_nocont(i, k)/pclc_nocont(i, k)*rhodz(i, k)
     595
     596          ! effective cloud droplet radius (microns) for liquid water clouds:
     597          ! For output diagnostics cloud droplet effective radius [um]
     598          ! we multiply here with f Effective radius of cloud droplet at top of cloud (m)* xl (fraction of liquid water
     599          ! clouds in the grid cell) to avoid problems in the averaging of the
     600          ! output.
     601          ! In the output of IOIPSL, derive the REAL cloud droplet
     602          ! effective radius as re/fl
     603
     604          fl(i, k) = pclc_nocont(i, k)*(1.-icefrac_optics(i,k))
     605          re(i, k) = rad_chaud(i, k)*fl(i, k)
     606
     607          rel = rad_chaud(i, k)
     608
     609          ! Calculation of ice cloud effective radius in micron
     610
     611          IF (iflag_rei .EQ. 2) THEN
     612              ! in-cloud ice water content in g/m3
     613              iwc = xfiwc_nocont(i, k) / pclc_nocont(i,k) * zrho(i,k) * 1000.
     614              ! this formula is a simplified version of the Sun 2001 one (as in the IFS model,
     615              ! and which is activated for iflag_rei = 1).
     616              ! In particular, the term in temperature**2 has been simplified.
     617              ! The new coefs are tunable, and are by default set so that the results fit well
     618              ! to the Sun 2001 formula
     619              ! By default, rei_coef = 2.4 and rei_min_temp = 175.
     620              ! The ranges of these parameters are determined so that the RMSE between this
     621              ! formula and the one from Sun 2001 is less than 4 times the minimum RMSE
     622              ! The ranges are [1.9, 2.9] for rei_coef and [160., 185.] for rei_min_temp
     623              dei = rei_coef * (iwc**0.2445) * ( temp(i,k) - rei_min_temp )
     624              ! we clip the results
     625              !deimin = 20.
     626              deimax = 155.
     627              !dei = MIN(MAX(dei, deimin), deimax)
     628              dei = MIN(dei, deimax)
     629              ! formula to convert effective diameter to effective radius
     630              rei = 3. * SQRT(3.) / 8. * dei
     631              rei = MAX(rei, rei_min)
     632             
     633          ELSEIF (iflag_rei .EQ. 1) THEN
     634              ! when we account for precipitation in the radiation scheme,
     635              ! we use the rei formula from Sun and Rikkus 1999 with a revision
     636              ! from Sun 2001 (as in the IFS model)
     637              iwc = xfiwc_nocont(i, k) / pclc_nocont(i,k) * zrho(i,k) * 1000. !in cloud ice water content in g/m3
     638              dei=(1.2351+0.0105*(temp(i,k)-273.15))*(45.8966*(iwc**0.2214) + &
     639                 &0.7957*(iwc**0.2535)*(temp(i,k)-83.15))
     640              !deimax=155.0
     641              !deimin=20.+40*cos(abs(latitude_deg(i))/180.*RPI)
     642              !Etienne: deimax and deimin controled by rei_max and rei_min in physiq.def
     643              deimax=rei_max*2.0
     644              deimin=2.0*rei_min+40*cos(abs(latitude_deg(i))/180.*RPI)
     645              dei=min(dei,deimax)
     646              dei=max(dei,deimin)
     647              rei=3.*sqrt(3.)/8.*dei
     648       
     649          ELSE
     650              ! Default
     651              ! for ice clouds: as a function of the ambiant temperature
     652              ! [formula used by Iacobellis and Somerville (2000), with an
     653              ! asymptotical value of 3.5 microns at T<-81.4 C added to be
     654              ! consistent with observations of Heymsfield et al. 1986]:
     655              ! 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as
     656              ! rei_max=61.29
     657              tc = temp(i, k) - 273.15
     658              rei = d_rei_dt*tc + rei_max
     659              IF (tc<=-81.4) rei = rei_min
     660          ENDIF
     661          ! -- cloud optical thickness :
     662          ! [for liquid clouds, traditional formula,
     663          ! for ice clouds, Ebert & Curry (1992)]
     664
     665
     666          IF (zflwp_var==0.) rel = 1.
     667          IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
     668
     669          !--Diagnostics of clouds emissivity, optical depth and ice crystal radius
     670          !--without contrails
     671          reice_nocont(i,k) = rei
     672          pcltau_nocont(i,k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/rei)
     673          ! -- cloud infrared emissivity:
     674          ! [the broadband infrared absorption coefficient is PARAMETERized
     675          ! as a function of the effective cld droplet radius]
     676          ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
     677          k_ice = k_ice0 + 1.0/rei
     678          pclemi_nocont(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
     679
     680        ELSE
     681
     682          rel = 1.
     683          rei = 1.
     684
     685          !--Diagnostics of clouds emissivity, optical depth and ice crystal radius
     686          !--without contrails
     687          reice_nocont(i,k) = 1.
     688          pclc_nocont(i,k) = 0.
     689          pclc(i,k) = contfra(i,k)
     690          pcltau_nocont(i, k) = 0.
     691          pclemi_nocont(i, k) = 0.
     692
     693        ENDIF
     694
     695
     696        IF ( contfra(i,k) .GT. 0.01 * seuil_neb ) THEN
     697          !--Everything is the same but with contrails
     698          zfiwp_var = 1000.*qice_cont(i, k)/contfra(i, k)*rhodz(i, k)
     699          pclc_var = contfra(i,k)
     700
     701          pcltau_cont(i, k) = zfiwp_var*(3.448E-03+2.431/(re_ice_crystals_contrails*1.E6))
     702          ! -- cloud infrared emissivity:
     703          ! [the broadband infrared absorption coefficient is PARAMETERized
     704          ! as a function of the effective cld droplet radius]
     705          ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
     706          k_ice = k_ice0 + 1.0/(re_ice_crystals_contrails*1.E6)
     707          pclemi_cont(i, k) = 1.0 - exp(-df*k_ice*zfiwp_var)
     708
     709        ELSE
     710          pclc_var = 0.
     711          pcltau_cont(i, k) = 0.
     712          pclemi_cont(i, k) = 0.
     713        ENDIF
     714
     715        rei = ( re_ice_crystals_contrails*1.E6 * pclc_var &
     716            + reice_nocont(i, k) * pclc_nocont(i,k) ) / ( pclc_var + pclc_nocont(i,k) )
     717
     718        zflwp_var = 1000.*xflwc(i, k)/pclc(i, k)*rhodz(i, k)
     719        zfiwp_var = 1000.*xfiwc(i, k)/pclc(i, k)*rhodz(i, k)
     720
     721        pcltau(i,k) = 3.0/2.0*(zflwp_var/rel) + zfiwp_var*(3.448E-03+2.431/rei)
     722        ! -- cloud infrared emissivity:
     723        ! [the broadband infrared absorption coefficient is PARAMETERized
     724        ! as a function of the effective cld droplet radius]
     725        ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
     726        k_ice = k_ice0 + 1.0/rei
     727        pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
     728
    581729      ENDIF
    582730
     731      reice(i, k) = rei
     732
     733      xflwp(i) = xflwp(i) + xflwc(i, k)*rhodz(i, k)
     734      xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
     735      xfiwp_nocont(i) = xfiwp_nocont(i) + xfiwc_nocont(i, k)*rhodz(i, k)
     736
    583737    ENDDO
    584738  ENDDO
     739
     740  ENDIF ! ok_plane_contrail
    585741
    586742  ! --if cloud droplet radius is fixed, then pcldtaupi=pcltau

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

    r5598 r5601  
    1919     tke, tke_dissip,                                   &
    2020     cell_area, stratomask,                             &
    21      cf_seri, rvc_seri, u_seri, v_seri,                 &
     21     cf_seri, qvc_seri, u_seri, v_seri,                 &
    2222     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
    2323     dcf_sub, dcf_con, dcf_mix,                         &
    2424     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
    2525     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
    26      rcont_seri, flight_dist, flight_h2o,               &
    27      contfra, radocond_cont, Tcritcont, qcritcont,      &
    28      potcontfraP, potcontfraNP, dcontfra_cir, dcf_avi,  &
    29      dqi_avi, dqvc_avi, cloudth_sth,cloudth_senv,       &
    30      cloudth_sigmath,cloudth_sigmaenv,                  &
     26     cfa_seri, qta_seri, flight_dist, flight_h2o,       &
     27     qice_cont, Tcritcont,                              &
     28     qcritcont, potcontfraP, potcontfraNP,              &
     29     cloudth_sth,                                       &
     30     cloudth_senv, cloudth_sigmath, cloudth_sigmaenv,   &
    3131     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
    3232     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
    … …  
    126126USE phys_local_var_mod, ONLY : issrfra100to150, issrfra150to200, issrfra200to250
    127127USE phys_local_var_mod, ONLY : issrfra250to300, issrfra300to400, issrfra400to500
     128USE phys_local_var_mod, ONLY : dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub
     129USE phys_local_var_mod, ONLY : dcfa_cir, dqta_cir, dcfa_mix, dqia_mix, dqta_mix
    128130
    129131IMPLICIT NONE
    … …  
    172174  !--------------------------------------------------
    173175  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
    174   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
     176  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qvc_seri         ! cloudy water vapor [kg/kg]
    175177  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
    176178  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
    … …  
    180182  ! INPUT/OUTPUT aviation
    181183  !--------------------------------------------------
    182   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rcont_seri       ! ratio of contrails fraction to total cloud fraction [-]
     184  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfa_seri         ! linear contrails fraction [-]
     185  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qta_seri         ! linear contrails total specific humidity [kg/kg]
    183186  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
    184187  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
    … …  
    238241  ! for contrails and aviation
    239242
    240   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: contfra        !--linear contrail fraction [-]
    241   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond_cont  !--condensed water in linear contrails used in the radiation scheme [kg/kg]
     243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_cont      !--condensed water in linear contrails used in the radiation scheme [kg/kg]
    242244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
    243245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
    244246  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: potcontfraP    !--potential persistent contrail fraction [-]
    245247  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: potcontfraNP   !--potential non-persistent contrail fraction [-]
    246   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcontfra_cir   !--linear contrail fraction to cirrus cloud fraction tendency [s-1]
    247   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
    248   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
    249   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
    250248
    251249
    … …  
    279277  ! LOCAL VARIABLES:
    280278  !----------------
    281   REAL, DIMENSION(klon) :: cldfra_in, qvc_in, qliq_in, qice_in
     279  REAL, DIMENSION(klon) :: qliq_in, qice_in
    282280  REAL, DIMENSION(klon,klev) :: ctot
    283281  REAL, DIMENSION(klon,klev) :: ctot_vol
    … …  
    324322  REAL :: delta_z
    325323  ! for contrails
    326   REAL, DIMENSION(klon) :: zqcont
     324  REAL, DIMENSION(klon) :: contfra, qcont
    327325  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrail)
    328326  ! Constants used for calculating ratios that are advected (using a parent-child
    … …  
    432430dqvc_con(:,:)   = 0.
    433431dqvc_mix(:,:)   = 0.
    434 contfra(:,:)    = 0.
    435 radocond_cont(:,:)= 0.
    436 Tcritcont(:,:)  = missing_val
    437 qcritcont(:,:)  = missing_val
    438 potcontfraP(:,:)= 0.
    439 potcontfraNP(:,:)= 0.
    440 dcontfra_cir(:,:)= 0.
    441 dcf_avi(:,:)    = 0.
    442 dqi_avi(:,:)    = 0.
    443 dqvc_avi(:,:)   = 0.
    444432qvc(:)          = 0.
    445433shear(:)        = 0.
    … …  
    491479DO k = klev, 1, -1
    492480
    493     !--Initialisation of advected properties
    494     cldfra_in(:) = cf_seri(:,k)
    495     qvc_in(:) = rvc_seri(:,k) * qt(:,k)
    496     qliq_in(:) = ql_seri(:,k)
    497     qice_in(:) = qi_seri(:,k)
    498 
    499481    IF (k.LE.klev-1) THEN
    500482       iftop=.false.
    … …  
    512494        zt(i)=temp(i,k)
    513495        zq(i)=qt(i,k)
     496        qliq_in(i) = ql_seri(i,k)
     497        qice_in(i) = qi_seri(i,k)
    514498        IF (.not. iftop) THEN
    515499           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
    … …  
    534518                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
    535519                        zqvapclr, zqupnew, zqice_sedim, &
    536                         cldfra_in, qvc_in, qliq_in, qice_in, &
     520                        cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, &
    537521                        zrfl, zrflclr, zrflcld, &
    538522                        zifl, ziflclr, ziflcld, &
    … …  
    759743                        klon, dtime, missing_val, &
    760744                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
    761                         cldfra_in, qvc_in, qliq_in, qice_in, &
     745                        cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, &
    762746                        shear, tke_dissip(:,k), cell_area, stratomask(:,k), &
    763747                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
    … …  
    766750                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
    767751                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
    768                         rcont_seri(:,k), flight_dist(:,k), flight_h2o(:,k), contfra(:,k), &
    769                         Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
    770                         dcontfra_cir(:,k), dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
     752                        cfa_seri(:,k), qta_seri(:,k), flight_dist(:,k), flight_h2o(:,k), &
     753                        contfra, qcont, Tcritcont(:,k), qcritcont(:,k), &
     754                        potcontfraP(:,k), potcontfraNP(:,k), &
     755                        dcfa_ini(:,k), dqia_ini(:,k), dqta_ini(:,k), &
     756                        dcfa_sub(:,k), dqia_sub(:,k), dqta_sub(:,k), &
     757                        dcfa_cir(:,k), dqta_cir(:,k), &
     758                        dcfa_mix(:,k), dqia_mix(:,k), dqta_mix(:,k))
    771759
    772760
    … …  
    963951      !--Contrails do not precipitate. We remove then from the variables temporarily
    964952      DO i = 1, klon
    965         IF (rneb(i,k) .GT. eps) THEN
    966           zqcont(i) = zcond(i) * zfice(i) * contfra(i,k) / rneb(i,k)
    967         ELSE
    968           zqcont(i) = 0.
    969         ENDIF
    970         rneb(i,k) = rneb(i,k) - contfra(i,k)
    971         zoliqi(i) = zoliqi(i) - zqcont(i)
     953        rneb(i,k) = rneb(i,k) - contfra(i)
     954        zoliqi(i) = zoliqi(i) - ( qcont(i) - zqs(i) * contfra(i) )
    972955      ENDDO
    973956    ENDIF
    … …  
    1009992      !--Contrails are reintroduced in the variables
    1010993      DO i = 1, klon
    1011         rneb(i,k) = rneb(i,k) + contfra(i,k)
    1012         zoliqi(i) = zoliqi(i) + zqcont(i)
     994        rneb(i,k) = rneb(i,k) + contfra(i)
     995        zoliqi(i) = zoliqi(i) + ( qcont(i) - zqs(i) * contfra(i) )
    1013996      ENDDO
    1014997    ENDIF
    … …  
    11041087   
    11051088    !--AB Write diagnostics and tracers for ice supersaturation
     1089    IF ( ok_plane_contrail ) THEN
     1090      cfa_seri(:,k) = contfra(:)
     1091      qta_seri(:,k) = qcont(:)
     1092      qice_cont(:,k) = qcont(:) - zqs(:) * contfra(:)
     1093    ENDIF
     1094
     1095    IF ( ok_ice_supersat .AND. .NOT. ok_unadjusted_clouds) qvc(:) = zqs(:) * rneb(:,k)
     1096
    11061097    IF ( ok_ice_supersat ) THEN
    11071098      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
    … …  
    11201111        ENDIF
    11211112
    1122         IF ( .NOT. ok_unadjusted_clouds ) THEN
    1123           qvc(i) = zqs(i) * rneb(i,k)
    1124         ENDIF
    1125         IF ( zq(i) .GT. min_qParent ) THEN
    1126           rvc_seri(i,k) = qvc(i) / zq(i)
    1127         ELSE
    1128           rvc_seri(i,k) = min_ratio
    1129         ENDIF
    1130         !--The MIN barrier is NEEDED because of:
    1131         !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
    1132         !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
    1133         !--The MAX barrier is a safeguard that should not be activated
    1134         rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
    1135         IF ( ok_plane_contrail ) THEN
    1136           IF ( rneb(i,k) .GT. min_qParent ) THEN
    1137             rcont_seri(i,k) = contfra(i,k) / rneb(i,k)
    1138           ELSE
    1139             rcont_seri(i,k) = min_ratio
    1140           ENDIF
    1141           !--This barrier should never be activated
    1142           rcont_seri(i,k) = MIN(MAX(rcont_seri(i,k), 0.), 1.)
    1143           radocond_cont(i,k) = radocond(i,k) * rcont_seri(i,k)
    1144         ENDIF
     1113        qvc_seri(i,k) = qvc(i)
    11451114
    11461115        !--Diagnostics

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

    r5589 r5601  
    9999      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
    100100      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
    101       rcont_seri, flight_dist, flight_h2o, contfra, &
     101      contfra_in, qcont_in, flight_dist, flight_h2o, contfra, qcont, &
    102102      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    103       dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi)
     103      dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub, &
     104      dcfa_cir, dqta_cir, dcfa_mix, dqia_mix, dqta_mix)
    104105
    105106!----------------------------------------------------------------------
    … …  
    116117!----------------------------------------------------------------------
    117118
    118 USE lmdz_lscp_tools,      ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
    119 USE lmdz_lscp_ini,        ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
    120 USE lmdz_lscp_ini,        ONLY: eps, temp_nowater, ok_weibull_warm_clouds
    121 USE lmdz_lscp_ini,        ONLY: ok_unadjusted_clouds
    122 USE lmdz_lscp_ini,        ONLY: ok_plane_contrail, ok_no_issr_strato
    123 
    124 USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp, &
    125                                 beta_pdf_lscp, temp_thresh_pdf_lscp, &
    126                                 std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, &
    127                                 coef_mixing_lscp, coef_shear_lscp, chi_mixing_lscp
    128 
    129 USE lmdz_aviation,        ONLY: contrails_mixing, contrails_formation_evolution
     119USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
     120USE lmdz_lscp_ini,   ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
     121USE lmdz_lscp_ini,   ONLY: eps, temp_nowater, ok_weibull_warm_clouds
     122USE lmdz_lscp_ini,   ONLY: ok_unadjusted_clouds, ok_no_issr_strato
     123
     124USE lmdz_lscp_ini,   ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp
     125USE lmdz_lscp_ini,   ONLY: beta_pdf_lscp, temp_thresh_pdf_lscp
     126USE lmdz_lscp_ini,   ONLY: std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp
     127USE lmdz_lscp_ini,   ONLY: coef_mixing_lscp, coef_shear_lscp, chi_mixing_lscp
     128
     129USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_contrails
     130USE lmdz_lscp_ini,   ONLY: coef_mixing_contrails, coef_shear_contrails
     131USE lmdz_lscp_ini,   ONLY: linear_contrails_lifetime
     132USE lmdz_aviation,   ONLY: contrails_formation
    130133
    131134IMPLICIT NONE
    … …  
    158161! Input for aviation
    159162!
    160 REAL,     INTENT(INOUT), DIMENSION(klon) :: rcont_seri    ! ratio of contrails fraction to total cloud fraction [-]
    161 REAL,     INTENT(IN),    DIMENSION(klon) :: flight_dist   ! aviation distance flown concentration [m/s/m3]
    162 REAL,     INTENT(IN),    DIMENSION(klon) :: flight_h2o    ! aviation emitted H2O concentration [kgH2O/s/m3]
     163REAL,     INTENT(IN)   , DIMENSION(klon) :: contfra_in    ! input linear contrails fraction [-]
     164REAL,     INTENT(IN)   , DIMENSION(klon) :: qcont_in      ! input linear contrails total specific humidity [kg/kg]
     165REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_dist   ! aviation distance flown concentration [m/s/m3]
     166REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_h2o    ! aviation emitted H2O concentration [kgH2O/s/m3]
    163167!
    164168!  Output
    … …  
    194198!
    195199REAL,     INTENT(INOUT), DIMENSION(klon) :: contfra      ! linear contrail fraction [-]
     200REAL,     INTENT(INOUT), DIMENSION(klon) :: qcont        ! linear contrail specific humidity [kg/kg]
    196201REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
    197202REAL,     INTENT(INOUT), DIMENSION(klon) :: qcritcont    ! critical specific humidity for contrail formation [kg/kg]
    198203REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraP  ! potential persistent contrail fraction [-]
    199204REAL,     INTENT(INOUT), DIMENSION(klon) :: potcontfraNP ! potential non-persistent contrail fraction [-]
    200 REAL,     INTENT(INOUT), DIMENSION(klon) :: dcontfra_cir ! linear contrail fraction to cirrus cloud fraction tendency [s-1]
    201 REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_avi      ! cloud fraction tendency because of aviation [s-1]
    202 REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_avi      ! specific ice content tendency because of aviation [kg/kg/s]
    203 REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_avi     ! specific cloud water vapor tendency because of aviation [kg/kg/s]
     205REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfa_ini     ! contrails cloud fraction tendency because of initial formation [s-1]
     206REAL,     INTENT(INOUT), DIMENSION(klon) :: dqia_ini     ! contrails ice specific humidity tendency because of initial formation [kg/kg/s]
     207REAL,     INTENT(INOUT), DIMENSION(klon) :: dqta_ini     ! contrails total specific humidity tendency because of initial formation [kg/kg/s]
     208REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfa_sub     ! contrails cloud fraction tendency because of sublimation [s-1]
     209REAL,     INTENT(INOUT), DIMENSION(klon) :: dqia_sub     ! contrails ice specific humidity tendency because of sublimation [kg/kg/s]
     210REAL,     INTENT(INOUT), DIMENSION(klon) :: dqta_sub     ! contrails total specific humidity tendency because of sublimation [kg/kg/s]
     211REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfa_cir     ! contrails cloud fraction tendency because of conversion in cirrus [s-1]
     212REAL,     INTENT(INOUT), DIMENSION(klon) :: dqta_cir     ! contrails total specific humidity tendency because of conversion in cirrus [kg/kg/s]
     213REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfa_mix     ! contrails cloud fraction tendency because of mixing [s-1]
     214REAL,     INTENT(INOUT), DIMENSION(klon) :: dqia_mix     ! contrails ice specific humidity tendency because of mixing [kg/kg/s]
     215REAL,     INTENT(INOUT), DIMENSION(klon) :: dqta_mix     ! contrails total specific humidity tendency because of mixing [kg/kg/s]
    204216!
    205217! Local
    … …  
    208220LOGICAL :: ok_warm_cloud
    209221REAL, DIMENSION(klon) :: qcld, qzero
     222REAL, DIMENSION(klon) :: clrfra, qclr
     223REAL, DIMENSION(klon) :: clrfra_pdf, qclr_pdf
    210224!
    211225! for lognormal
    … …  
    214228!
    215229! for unadjusted clouds
    216 REAL :: qvapincld, qvapincld_new
    217 REAL :: qiceincld
     230REAL :: qvapincld, qiceincld, qvapincld_new
    218231!
    219232! for sublimation
    220 REAL :: pdf_alpha
    221233REAL :: dqt_sub
    222234!
    223235! for condensation
    224236REAL, DIMENSION(klon) :: qsatl, dqsatl
    225 REAL :: clrfra, qclr, sl_clr, rhl_clr
    226 REAL :: pdf_ratqs, pdf_skew, pdf_scale, pdf_loc
     237REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_loc
     238REAL :: rhl_clr
     239REAL :: pdf_ratqs, pdf_skew
    227240REAL :: pdf_x, pdf_y, pdf_T
    228241REAL :: pdf_e3, pdf_e4
    … …  
    230243!
    231244! for mixing
    232 REAL, DIMENSION(klon) :: subfra, qsub
    233 REAL :: dqt_mix_sub, dqt_mix_issr
    234 REAL :: dcf_mix_sub, dcf_mix_issr
    235 REAL :: dqvc_mix_sub, dqvc_mix_issr
    236 REAL :: dqt_mix
    237245REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
    238 REAL :: envfra_mix, cldfra_mix
     246REAL :: clrfra_mix, sigma_mix
    239247REAL :: L_shear, shear_fra
    240 REAL :: sigma_mix, issrfra_mix, subfra_mix, qvapinmix
     248REAL :: qvapinclr_lim
     249REAL :: pdf_fra_above_lim, pdf_q_above_lim
     250REAL :: pdf_fra_below_lim, pdf_q_below_lim
    241251!
    242252! for cell properties
    243 REAL :: rho, rhodz, dz
     253REAL, DIMENSION(klon) :: dz
    244254
    245255qzero(:) = 0.
    … …  
    261271    issrfra(i)  = 0.
    262272    qissr(i)    = 0.
     273    contfra(i)  = 0.
     274    qcont(i)    = 0.
    263275
    264276    !--If the temperature is higher than the threshold below which
    … …  
    298310
    299311      !--Initialisation
    300       IF ( temp(i) .GT. temp_nowater ) THEN
    301         !--If the air mass is warm (liquid water can exist),
    302         !--all the memory is lost and the scheme becomes statistical,
    303         !--i.e., the sublimation and mixing processes are deactivated,
    304         !--and the condensation process is slightly adapted
    305         !--This can happen only if ok_weibull_warm_clouds = .TRUE.
    306         cldfra(i) = MAX(0., MIN(1., cldfra_in(i)))
    307         qcld(i)   = MAX(0., MIN(qtot(i), qliq_in(i) + qice_in(i) + qvc_in(i)))
    308         qvc(i)    = MAX(0., MIN(qcld(i), qvc_in(i)))
    309         ok_warm_cloud = .TRUE.
    310       ELSE
    311         !--The following barriers ensure that the traced cloud properties
    312         !--are consistent. In some rare cases, i.e. the cloud water vapor
    313         !--can be greater than the total water in the gridbox
    314         cldfra(i) = MAX(0., MIN(1., cldfra_in(i)))
    315         qcld(i)   = MAX(0., MIN(qtot(i), qvc_in(i) + qice_in(i)))
    316         qvc(i)    = MAX(0., MIN(qcld(i), qvc_in(i)))
    317         ok_warm_cloud = .FALSE.
    318       ENDIF
     312      !--If the air mass is warm (liquid water can exist),
     313      !--all the memory is lost and the scheme becomes statistical,
     314      !--i.e., the sublimation and mixing processes are deactivated,
     315      !--and the condensation process is slightly adapted
     316      !--This can happen only if ok_weibull_warm_clouds = .TRUE.
     317      ok_warm_cloud = ( temp(i) .GT. temp_nowater )
     318
     319      !--The following barriers ensure that the traced cloud properties
     320      !--are consistent. In some rare cases, i.e. the cloud water vapor
     321      !--can be greater than the total water in the gridbox
     322      cldfra(i) = MAX(0., MIN(1., cldfra_in(i)))
     323      qcld(i)   = MAX(0., MIN(qtot(i), qliq_in(i) + qice_in(i) + qvc_in(i)))
     324      qvc(i)    = MAX(0., MIN(qcld(i), qvc_in(i)))
     325
     326      !--We init clear fraction properties
     327      clrfra(i) = 1. - cldfra(i)
     328      qclr(i) = qtot(i) - qcld(i)
    319329
    320330      dcf_sub(i)  = 0.
    … …  
    332342      !--Initialisation of the cell properties
    333343      !--Dry density [kg/m3]
    334       rho = pplay(i) / temp(i) / RD
     344      !rho = pplay(i) / temp(i) / RD
    335345      !--Dry air mass [kg/m2]
    336       rhodz = ( paprsdn(i) - paprsup(i) ) / RG
     346      !rhodz = ( paprsdn(i) - paprsup(i) ) / RG
    337347      !--Cell thickness [m]
    338       dz = rhodz / rho
     348      !dz = rhodz / rho
     349      dz(i) = ( paprsdn(i) - paprsup(i) ) / pplay(i) * temp(i) * RD / RG
     350
     351      !--If contrails are activated
     352      IF ( ok_plane_contrail ) THEN
     353        contfra(i) = contfra_in(i)
     354        qcont(i) = qcont_in(i)
     355
     356        dcfa_ini(i) = 0.
     357        dqia_ini(i) = 0.
     358        dqta_ini(i) = 0.
     359        dcfa_sub(i) = 0.
     360        dqia_sub(i) = 0.
     361        dqta_sub(i) = 0.
     362        dcfa_cir(i) = 0.
     363        dqta_cir(i) = 0.
     364        dcfa_mix(i) = 0.
     365        dqia_mix(i) = 0.
     366        dqta_mix(i) = 0.
     367      ENDIF
     368
     369
     370      !----------------------------------------------------------------------
     371      !--    SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CONTRAIL    --
     372      !----------------------------------------------------------------------
     373
     374      !--If there is a contrail
     375      IF ( contfra(i) .GT. eps ) THEN
     376        !--The contrail is always adjusted to saturation
     377        qiceincld = ( qcont(i) / contfra(i) - qsat(i) )
     378       
     379        !--If the ice water content is too low, the cloud is purely sublimated
     380        IF ( qiceincld .LT. eps ) THEN
     381          dcfa_sub(i) = - contfra(i)
     382          dqia_sub(i) = - qiceincld * contfra(i)
     383          dqta_sub(i) = - qcont(i)
     384          contfra(i) = 0.
     385          qcont(i) = 0.
     386        ELSE
     387          !--We remove contrails from the main class
     388          cldfra(i) = cldfra(i) - contfra(i)
     389          qcld(i) = qcld(i) - qcont(i)
     390          qvc(i) = qvc(i) - qsat(i) * contfra(i)
     391        ENDIF   ! qiceincld .LT. eps
     392      ENDIF ! contfra(i) .GT. eps
    339393
    340394
    … …  
    345399      !--If there is a cloud
    346400      IF ( cldfra(i) .GT. eps ) THEN
    347        
     401
    348402        qvapincld = qvc(i) / cldfra(i)
    349403        qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
    … …  
    375429
    376430          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    377             CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
    378                                         pplay(i), dtime, qvapincld_new)
     431            qvapincld_new = QVAPINCLD_DEPSUB( &
     432                qvapincld, qiceincld, temp(i), qsat(i), pplay(i), dtime)
    379433            IF ( qvapincld_new .EQ. 0. ) THEN
    380434              !--If all the ice has been sublimated, we sublimate
    381               !--completely the cloud and do not activate the sublimation
     435              !--completely the cloud and do not activate the dissipation
    382436              !--process
    383437              !--Tendencies and diagnostics
    … …  
    443497      !--This section relies on a distribution of water in the clear-sky region of
    444498      !--the mesh.
     499      clrfra_pdf(i) = 1. - cldfra(i) - contfra(i)
     500      qclr_pdf(i) = qtot(i) - qcld(i) - qcont(i)
    445501
    446502      !--If there is a clear-sky region
    447       IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
    448 
    449         !--Water quantity in the clear-sky + potential liquid cloud (gridbox average)
    450         qclr = qtot(i) - qcld(i)
     503      IF ( clrfra_pdf(i) .GT. eps ) THEN
    451504
    452505        !--New PDF
    453         rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100.
     506        rhl_clr = qclr_pdf(i) / clrfra_pdf(i) / qsatl(i) * 100.
    454507
    455508        !--Calculation of the properties of the PDF
    … …  
    463516        !       * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
    464517        !--  tuning coef * (cell area**0.25) * (function of temperature)
    465         pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 &
     518        pdf_e1 = beta_pdf_lscp * ( clrfra_pdf(i) * cell_area(i) )**0.25 &
    466519                * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
    467520        IF ( rhl_clr .GT. 50. ) THEN
    … …  
    471524        ENDIF
    472525        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
    473         pdf_alpha = EXP( MIN(1000., rhl_clr) / 100. ) * pdf_e3
     526        pdf_alpha(i) = EXP( MIN(1000., rhl_clr) / 100. ) * pdf_e3
    474527       
    475528        IF ( ok_warm_cloud ) THEN
    … …  
    480533        ENDIF
    481534
    482         pdf_e2 = GAMMA(1. + 1. / pdf_alpha)
    483         pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 ))
    484         pdf_loc = rhl_clr - pdf_scale * pdf_e2
     535        pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
     536        pdf_scale(i) = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha(i)) - pdf_e2**2 ))
     537        pdf_loc(i) = rhl_clr - pdf_scale(i) * pdf_e2
    485538
    486539        !--Calculation of the newly condensed water and fraction (pronostic)
    … …  
    489542
    490543        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
    491         pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
    492         pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
    493         pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
     544        pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     545        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     546        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
    494547        cf_cond = EXP( - pdf_y )
    495         qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100.
     548        qt_cond = ( pdf_e3 + pdf_loc(i) * cf_cond ) * qsatl(i) / 100.
    496549
    497550        IF ( cf_cond .GT. eps ) THEN
    498551
    499           dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond
    500           dqt_con = ( 1. - cldfra(i) ) * qt_cond
    501          
    502           !--Barriers
    503           dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i))
    504           dqt_con = MIN(dqt_con, qclr)
     552          dcf_con(i) = clrfra_pdf(i) * cf_cond
     553          dqt_con = clrfra_pdf(i) * qt_cond
    505554
    506555          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    … …  
    510559            qvapincld = gamma_cond(i) * qsat(i)
    511560            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
    512             CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), &
    513                                         pplay(i), dtime / 2., qvapincld_new)
     561            qvapincld_new = QVAPINCLD_DEPSUB( &
     562                qvapincld, qiceincld, temp(i), qsat(i), pplay(i), dtime / 2.)
    514563            qvapincld = qvapincld_new
    515564          ELSE
    … …  
    527576          qcld(i)   = MIN(qtot(i), qcld(i) + dqt_con)
    528577          qvc(i)    = MIN(qcld(i), qvc(i) + dqvc_con(i))
     578          clrfra(i) = MAX(0., clrfra(i) - dcf_con(i))
     579          qclr(i)   = MAX(0., qclr(i) - dqt_con)
    529580
    530581        ENDIF ! ok_warm_cloud, cf_cond .GT. eps
    … …  
    533584        !--and lower than gamma_cond * qsat, which is the ice supersaturated region
    534585        pdf_x = qsat(i) / qsatl(i) * 100.
    535         pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha
    536         pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y )
    537         pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2
    538         issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) )
    539         qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
     586        pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     587        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     588        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     589        issrfra(i) = EXP( - pdf_y ) * clrfra_pdf(i)
     590        qissr(i) = ( pdf_e3 * clrfra_pdf(i) + pdf_loc(i) * issrfra(i) ) * qsatl(i) / 100.
    540591
    541592        issrfra(i) = issrfra(i) - dcf_con(i)
    … …  
    543594
    544595      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
    545 
    546       !--Calculation of the subsaturated clear sky fraction and water
    547       subfra(i) = 1. - cldfra(i) - issrfra(i)
    548       qsub(i) = qtot(i) - qcld(i) - qissr(i)
    549596
    550597
    … …  
    556603      !--but does not cover the entire mesh.
    557604      !
    558       IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN
    559         !--Initialisation
    560         dcf_mix_sub   = 0.
    561         dqt_mix_sub   = 0.
    562         dqvc_mix_sub  = 0.
    563         dcf_mix_issr  = 0.
    564         dqt_mix_issr  = 0.
    565         dqvc_mix_issr = 0.
     605      IF ( ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN
    566606
    567607        !-- PART 1 - TURBULENT DIFFUSION
    … …  
    614654        !--The fraction of clear sky mixed is
    615655        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
    616         envfra_mix = N_cld_mix * RPI / cell_area(i) &
     656        clrfra_mix = N_cld_mix * RPI / cell_area(i) &
    617657                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
    618         !--The fraction of cloudy sky mixed is
    619         !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
    620         cldfra_mix = N_cld_mix * RPI / cell_area(i) &
    621                    * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
    622658
    623659
    … …  
    628664        !--semi-major axis (a_mix), on the entire cell heigh dz.
    629665        !--The increase in size is
    630         L_shear = coef_shear_lscp * shear(i) * dz * dtime
     666        L_shear = coef_shear_lscp * shear(i) * dz(i) * dtime
    631667        !--therefore, the fraction of clear sky mixed is
    632668        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
    … …  
    639675        !--We therefore assume that the sheared clouds and the turbulent
    640676        !--mixed clouds are different.
    641         envfra_mix = envfra_mix + shear_fra
    642         cldfra_mix = cldfra_mix + shear_fra
     677        clrfra_mix = clrfra_mix + shear_fra
    643678
    644679
    645680        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
    646681
    647         !--The environment fraction is allocated to subsaturated sky or supersaturated sky,
    648         !--according to the factor sigma_mix. This is computed as the ratio of the
    649         !--subsaturated sky fraction to the environment fraction, corrected by a factor
    650         !--chi_mixing_lscp for the supersaturated part. If chi is greater than 1, the
    651         !--supersaturated sky is favoured. Physically, this means that it is more likely
    652         !--to have supersaturated sky around the cloud than subsaturated sky.
    653         sigma_mix = subfra(i) / ( subfra(i) + chi_mixing_lscp * issrfra(i) )
    654         subfra_mix = MIN( sigma_mix * envfra_mix, subfra(i) )
    655         issrfra_mix = MIN( ( 1. - sigma_mix ) * envfra_mix, issrfra(i) )
    656         cldfra_mix = MIN( cldfra_mix, cldfra(i) )
    657 
    658         !--First, we mix the subsaturated sky (subfra_mix) and the cloud close
    659         !--to this fraction (sigma_mix * cldfra_mix).
    660         IF ( subfra(i) .GT. eps ) THEN
    661           !--We compute the total humidity in the mixed air, which
    662           !--can be either sub- or supersaturated.
    663           qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &
    664                     + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &
    665                     / ( subfra_mix + cldfra_mix * sigma_mix )
    666 
    667           IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    668             qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) &
    669                       / ( subfra_mix + cldfra_mix * sigma_mix )
    670             CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
    671                                         pplay(i), dtime, qvapincld_new)
    672             IF ( qvapincld_new .EQ. 0. ) THEN
    673               !--If all the ice has been sublimated, we sublimate
    674               !--completely the cloud and do not activate the sublimation
    675               !--process
    676               !--Tendencies and diagnostics
    677               dcf_mix_sub = - sigma_mix * cldfra_mix
    678               dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
    679               dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i)
    680             ELSE
    681               dcf_mix_sub = subfra_mix
    682               dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
    683               dqvc_mix_sub = dcf_mix_sub * qvapincld_new
    684             ENDIF
    685           ELSE
    686             IF ( qvapinmix .GT. qsat(i) ) THEN
    687               !--If the mixed air is supersaturated, we condense the subsaturated
    688               !--region which was mixed.
    689               dcf_mix_sub = subfra_mix
    690               dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i)
    691               dqvc_mix_sub = dcf_mix_sub * qsat(i)
    692             ELSE
    693               !--Else, we sublimate the cloud which was mixed.
    694               dcf_mix_sub = - sigma_mix * cldfra_mix
    695               dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i)
    696               dqvc_mix_sub = dcf_mix_sub * qsat(i)
    697             ENDIF
    698           ENDIF ! ok_unadjusted_clouds
    699         ENDIF ! subfra .GT. eps
    700    
    701         !--We then mix the supersaturated sky (issrfra_mix) and the cloud,
    702         !--for which the mixed air is always supersatured, therefore
    703         !--the cloud necessarily expands
    704         IF ( issrfra(i) .GT. eps ) THEN
    705 
    706           IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    707             qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) &
    708                       + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) &
    709                       / ( issrfra_mix + cldfra_mix * (1. -  sigma_mix) )
    710             qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) &
    711                       / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) )
    712             CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), &
    713                                         pplay(i), dtime, qvapincld_new)
    714             dcf_mix_issr = issrfra_mix
    715             dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
    716             dqvc_mix_issr = dcf_mix_issr * qvapincld_new
    717           ELSE
    718             !--In this case, the additionnal vapor condenses
    719             dcf_mix_issr = issrfra_mix
    720             dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i)
    721             dqvc_mix_issr = dcf_mix_issr * qsat(i)
    722           ENDIF ! ok_unadjusted_clouds
    723 
    724         ENDIF ! issrfra .GT. eps
    725 
    726         !--Sum up the tendencies from subsaturated sky and supersaturated sky
    727         dcf_mix(i)  = dcf_mix_sub  + dcf_mix_issr
    728         dqt_mix     = dqt_mix_sub  + dqt_mix_issr
    729         dqvc_mix(i) = dqvc_mix_sub + dqvc_mix_issr
    730         dqi_mix(i)  = dqt_mix - dqvc_mix(i)
    731 
    732         IF ( ok_plane_contrail .AND. ( rcont_seri(i) .GT. eps ) ) THEN
    733           CALL contrails_mixing( &
    734               dtime, pplay(i), shear(i), pbl_eps(i), cell_area(i), dz, &
    735               temp(i), qtot(i), qsat(i), subfra(i), qsub(i), issrfra(i), qissr(i), &
    736               cldfra(i), qcld(i), qvc(i), rcont_seri(i), &
    737               dcf_mix(i), dqvc_mix(i), dqi_mix(i), dqt_mix, dcf_mix_issr, dqt_mix_issr)
     682        clrfra_mix = MIN( clrfra_mix, clrfra(i) )
     683
     684        !--We compute the limit vapor in clear sky where the mixed cloud could not
     685        !--survive if all the ice crystals were sublimated. Note that here we assume,
     686        !--for growth or reduction of the cloud, that the mixed cloud is adjusted
     687        !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a
     688        !--diagnostic, and if the cloud size is increased, we add the new vapor to the
     689        !--cloud's vapor without condensing or sublimating ice crystals
     690        qvapinclr_lim = ( ( cldfra(i) + clrfra_mix ) * qsat(i) - qcld(i) ) / clrfra_mix
     691
     692        IF ( qvapinclr_lim .LT. 0. ) THEN
     693          !--Whatever we do, the cloud will increase in size
     694          dcf_mix(i)  = clrfra_mix
     695          dqvc_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
     696          dqi_mix(i)  = 0.
     697        ELSE
     698          !--We then calculate the clear sky part where the humidity is lower than
     699          !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim
     700          !--This is the clear-sky PDF calculated in the condensation section. Note
     701          !--that if we are here, we necessarily went through the condensation part
     702          !--because the clear sky fraction can only be reduced by condensation.
     703          !--Thus the `pdf_xxx` variables are well defined.
     704          pdf_x = qvapinclr_lim / qsatl(i) * 100.
     705          pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     706          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     707          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     708          pdf_fra_above_lim = EXP( - pdf_y ) * clrfra_pdf(i)
     709          pdf_q_above_lim = ( pdf_e3 * clrfra_pdf(i) &
     710                          + pdf_loc(i) * pdf_fra_above_lim ) * qsatl(i) / 100.
     711
     712          pdf_fra_below_lim = clrfra_pdf(i) - pdf_fra_above_lim
     713          pdf_q_below_lim = qclr_pdf(i) - pdf_q_above_lim
     714          !--We remove the already condensated part (it is well defined)
     715          pdf_fra_above_lim = pdf_fra_above_lim - dcf_con(i)
     716          pdf_q_above_lim = pdf_q_above_lim - dqvc_con(i) - dqi_con(i)
     717
     718          !--sigma_mix quantifies
     719          sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim )
     720
     721          IF ( pdf_fra_above_lim .GT. eps ) THEN
     722            dcf_mix(i)  = clrfra_mix * sigma_mix
     723            dqvc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
     724            dqi_mix(i)  = 0.
     725          ENDIF
     726
     727          IF ( pdf_fra_below_lim .GT. eps ) THEN
     728            dcf_mix(i)  = dcf_mix(i)  - cldfra(i) * ( 1. - sigma_mix )
     729            dqvc_mix(i) = dqvc_mix(i) - cldfra(i) * ( 1. - sigma_mix ) &
     730                                                  * qvc(i) / cldfra(i)
     731            dqi_mix(i)  = dqi_mix(i)  - cldfra(i) * ( 1. - sigma_mix ) &
     732                                                  * ( qcld(i) - qvc(i) ) / cldfra(i)
     733          ENDIF
    738734        ENDIF
    739735
    740736        !--Add tendencies
    741         issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr)
    742         qissr(i)   = MAX(0., qissr(i) - dqt_mix_issr)
    743737        cldfra(i)  = MAX(0., MIN(1., cldfra(i) + dcf_mix(i)))
    744         qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix))
     738        qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqvc_mix(i) + dqi_mix(i)))
    745739        qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))
    746 
    747       ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) )
    748 
    749 
    750       !----------------------------------------
    751       !--       CONTRAILS AND AVIATION       --
    752       !----------------------------------------
    753 
    754       IF ( ok_plane_contrail ) THEN
    755         CALL contrails_formation_evolution( &
    756             dtime, pplay(i), temp(i), qsat(i), qsatl(i), gamma_cond(i), &
    757             rcont_seri(i), flight_dist(i), cldfra(i), qvc(i), &
    758             pdf_loc, pdf_scale, pdf_alpha, &
    759             Tcritcont(i), qcritcont(i), potcontfraP(i), potcontfraNP(i), contfra(i), &
    760             dcontfra_cir(i), dcf_avi(i), dqvc_avi(i), dqi_avi(i) &
    761             )
     740        clrfra(i)  = MAX(0., MIN(1., clrfra(i) - dcf_mix(i)))
     741        qclr(i)    = MAX(0., MIN(qtot(i), qclr(i) - dqvc_mix(i) - dqi_mix(i)))
     742     
     743      ENDIF ! ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )
     744
     745      !--------------------------------------
     746      !--        CONTRAIL MIXING           --
     747      !--------------------------------------
     748
     749      IF ( ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps ) ) THEN
     750
     751        !-- PART 1 - TURBULENT DIFFUSION
     752       
     753        !--Clouds within the mesh are assumed to be ellipses. The length of the
     754        !--semi-major axis is a and the length of the semi-minor axis is b.
     755        !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer.
     756        !--In particular, it is 0 if cldfra = 1.
     757        !--clouds_perim is the total perimeter of the clouds within the mesh,
     758        !--not considering interfaces with other meshes (only the interfaces with clear
     759        !--sky are taken into account).
     760        !--
     761        !--The area of each cloud is A = a * b * RPI,
     762        !--and the perimeter of each cloud is
     763        !-- P ~= RPI * ( 3 * (a + b) - SQRT( (3 * a + b) * (a + 3 * b) ) )
     764        !--
     765        !--With cell_area the area of the cell, we have:
     766        !-- cldfra = A * N_cld_mix / cell_area
     767        !-- clouds_perim = P * N_cld_mix
     768        !--
     769        !--We assume that the ratio between b and a is a function of
     770        !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
     771        !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
     772        !--are spherical.
     773        !-- b / a = bovera = MAX(0.1, cldfra)
     774        bovera = aspect_ratio_contrails
     775        !--P / a is a function of b / a only, that we can calculate
     776        !-- P / a = RPI * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) )
     777        Povera = RPI * ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )
     778     
     779        !--The clouds perimeter is imposed using the formula from Morcrette 2012,
     780        !--based on observations.
     781        !-- clouds_perim / cell_area = N_cld_mix * ( P / a * a ) / cell_area = coef_mix_lscp * cldfra * ( 1. - cldfra )
     782        !--With cldfra = a * ( b / a * a ) * RPI * N_cld_mix / cell_area, we have:
     783        !-- cldfra = a * b / a * RPI / (P / a) * coef_mix_lscp * cldfra * ( 1. - cldfra )
     784        !-- a = (P / a) / ( coef_mix_lscp * RPI * ( 1. - cldfra ) * (b / a) )
     785        a_mix = Povera / coef_mixing_contrails / RPI / ( 1. - contfra(i) ) / bovera
     786        !--and finally,
     787        !-- N_cld_mix = coef_mix_lscp * cldfra * ( 1. - cldfra ) * cell_area / ( P / a * a )
     788        N_cld_mix = coef_mixing_contrails * contfra(i) * ( 1. - contfra(i) ) * cell_area(i) &
     789                  / Povera / a_mix
     790       
     791        !--The time required for turbulent diffusion to homogenize a region of size
     792        !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)
     793        !--We compute L_mix and assume that the cloud is mixed over this length
     794        L_mix = SQRT( dtime**3 * pbl_eps(i) )
     795        !--The mixing length cannot be greater than the semi-minor axis. In this case,
     796        !--the entire cloud is mixed.
     797        L_mix = MIN(L_mix, a_mix * bovera)
     798       
     799        !--The fraction of clear sky mixed is
     800        !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area
     801        clrfra_mix = N_cld_mix * RPI / cell_area(i) &
     802                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
     803       
     804       
     805        !-- PART 2 - SHEARING
     806       
     807        !--The clouds are then sheared. We keep the shape and number
     808        !--assumptions from before. The clouds are sheared with a random orientation
     809        !--of the wind, on average we assume that the wind and the semi-major axis
     810        !--make a 45 degrees angle. Moreover, the contrails only mix
     811        !--along their semi-minor axis (b), because it is easier to compute.
     812        !--With this, the clouds increase in size along b only, by a factor
     813        !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind)
     814        L_shear = coef_shear_contrails * shear(i) * dz(i) * dtime
     815        !--therefore, the fraction of clear sky mixed is
     816        !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area
     817        !--and the fraction of cloud mixed is
     818        !-- N_cld_mix * ( a * b - a * (b - L_shear * SQRT(2.) / 2.) ) * RPI / 2. / cell_area
     819        !--(note that they are equal)
     820        shear_fra = RPI * L_shear * a_mix * SQRT(2.) / 2. / 2. * N_cld_mix / cell_area(i)
     821        !--and the environment and cloud mixed fractions are the same,
     822        !--which we add to the previous calculated mixed fractions.
     823        !--We therefore assume that the sheared clouds and the turbulent
     824        !--mixed clouds are different.
     825        clrfra_mix = clrfra_mix + shear_fra
     826       
     827       
     828        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
     829       
     830        clrfra_mix = MIN( clrfra_mix, clrfra(i) )
     831       
     832        !--We compute the limit vapor in clear sky where the mixed cloud could not
     833        !--survive if all the ice crystals were sublimated. Note that here we assume,
     834        !--for growth or reduction of the cloud, that the mixed cloud is adjusted
     835        !--to saturation, ie the vapor in the mixed cloud is qsat. This is only a
     836        !--diagnostic, and if the cloud size is increased, we add the new vapor to the
     837        !--cloud's vapor without condensing or sublimating ice crystals
     838        qvapinclr_lim = ( ( contfra(i) + clrfra_mix ) * qsat(i) - qcont(i) ) / clrfra_mix
     839
     840        IF ( qvapinclr_lim .LT. 0. ) THEN
     841          !--Whatever we do, the cloud will increase in size
     842          dcfa_mix(i) = clrfra_mix
     843          dqta_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
     844          dqia_mix(i) = 0.
     845        ELSE
     846          !--We then calculate the clear sky part where the humidity is lower than
     847          !--qvapinclr_lim, and the part where it is higher than qvapinclr_lim
     848          !--This is the clear-sky PDF calculated in the condensation section. Note
     849          !--that if we are here, we necessarily went through the condensation part
     850          !--because the clear sky fraction can only be reduced by condensation.
     851          !--Thus the `pdf_xxx` variables are well defined.
     852          pdf_x = qvapinclr_lim / qsatl(i) * 100.
     853          pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     854          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     855          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     856          pdf_fra_above_lim = EXP( - pdf_y ) * clrfra_pdf(i)
     857          pdf_q_above_lim = ( pdf_e3 * clrfra_pdf(i) &
     858                          + pdf_loc(i) * pdf_fra_above_lim ) * qsatl(i) / 100.
     859
     860          pdf_fra_below_lim = clrfra_pdf(i) - pdf_fra_above_lim
     861          pdf_q_below_lim = qclr_pdf(i) - pdf_q_above_lim
     862          !--We remove the already condensated part (it is well defined)
     863          pdf_fra_above_lim = pdf_fra_above_lim - dcf_con(i)
     864          pdf_q_above_lim = pdf_q_above_lim - dqvc_con(i) - dqi_con(i)
     865
     866          !--sigma_mix quantifies
     867          sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim )
     868
     869          IF ( pdf_fra_above_lim .GT. eps ) THEN
     870            dcfa_mix(i) = clrfra_mix * sigma_mix
     871            dqta_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
     872            dqia_mix(i) = 0.
     873          ENDIF
     874
     875          IF ( pdf_fra_below_lim .GT. eps ) THEN
     876            dcfa_mix(i) = dcfa_mix(i) - contfra(i) * ( 1. - sigma_mix )
     877            dqta_mix(i) = dqta_mix(i) - contfra(i) * ( 1. - sigma_mix ) &
     878                                      * qcont(i) / contfra(i)
     879            dqia_mix(i) = dqia_mix(i) - contfra(i) * ( 1. - sigma_mix ) &
     880                                      * ( qcont(i) / contfra(i) - qsat(i) )
     881          ENDIF
     882        ENDIF
    762883
    763884        !--Add tendencies
    764         issrfra(i) = MAX(0., issrfra(i) - dcf_avi(i))
    765         qissr(i)   = MAX(0., qissr(i) - dqvc_avi(i) - dqi_avi(i))
    766         cldfra(i)  = MIN(1., cldfra(i) + dcf_avi(i))
    767         qcld(i)    = MIN(qtot(i), qcld(i) + dqvc_avi(i) + dqi_avi(i))
    768         qvc(i)     = MIN(qcld(i), qvc(i) + dqvc_avi(i))
    769       ENDIF
    770 
    771       !-------------------------------------------
    772       !--       FINAL BARRIERS AND OUTPUTS      --
    773       !-------------------------------------------
    774 
    775       IF ( cldfra(i) .LT. eps ) THEN
    776         !--If the cloud is too small, it is sublimated.
    777         cldfra(i) = 0.
    778         qcld(i)   = 0.
    779         qvc(i)    = 0.
    780         qincld(i) = qsat(i)
    781         IF ( ok_plane_contrail ) contfra(i) = 0.
    782       ELSE
    783         qincld(i) = qcld(i) / cldfra(i)
    784       ENDIF ! cldfra .LT. eps
     885        contfra(i) = MAX(0., MIN(1., contfra(i) + dcfa_mix(i)))
     886        qcont(i)   = MAX(0., MIN(qtot(i), qcont(i) + dqta_mix(i)))
     887        clrfra(i)  = MAX(0., MIN(1., clrfra(i) - dcfa_mix(i)))
     888        qclr(i)    = MAX(0., MIN(qtot(i), qclr(i) - dqta_mix(i)))
     889
     890      ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
     891
     892      !--We put back contrails in the clouds class
     893      cldfra(i) = cldfra(i) + contfra(i)
     894      qcld(i) = qcld(i) + qcont(i)
     895      qvc(i) = qvc(i) + qsat(i) * contfra(i)
     896
    785897
    786898      !--Diagnostics
    … …  
    796908      dqvc_con(i) = dqvc_con(i) / dtime
    797909      dqvc_mix(i) = dqvc_mix(i) / dtime
    798       IF ( ok_plane_contrail ) THEN
    799         dcontfra_cir(i) = dcontfra_cir(i) / dtime
    800         dcf_avi(i)  = dcf_avi(i)  / dtime
    801         dqi_avi(i)  = dqi_avi(i)  / dtime
    802         dqvc_avi(i) = dqvc_avi(i) / dtime
     910
     911      !-------------------------------------------
     912      !--       FINAL BARRIERS AND OUTPUTS      --
     913      !-------------------------------------------
     914
     915      IF ( cldfra(i) .LT. eps ) THEN
     916        !--If the cloud is too small, it is sublimated.
     917        cldfra(i) = 0.
     918        qcld(i)   = 0.
     919        qvc(i)    = 0.
     920        qincld(i) = qsat(i)
     921      ELSE
     922        qincld(i) = qcld(i) / cldfra(i)
     923      ENDIF ! cldfra .LT. eps
     924
     925    ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
     926
     927  ENDIF   ! end keepgoing
     928
     929ENDDO     ! end loop on i
     930
     931
     932!----------------------------------------
     933!--       CONTRAILS AND AVIATION       --
     934!----------------------------------------
     935IF ( ok_plane_contrail ) THEN
     936
     937  CALL contrails_formation( &
     938      klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, &
     939      flight_dist, cldfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &
     940      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
     941      dcfa_ini, dqia_ini, dqta_ini)
     942
     943  DO i = 1, klon
     944    IF ( keepgoing(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
     945
     946      !--Convert existing contrail fraction into "natural" cirrus cloud fraction
     947      dcfa_cir(i) = contfra(i) * ( EXP( - dtime / linear_contrails_lifetime ) - 1. )
     948      dqta_cir(i) = qcont(i)   * ( EXP( - dtime / linear_contrails_lifetime ) - 1. )
     949
     950      !--Add tendencies
     951      issrfra(i) = MAX(0., issrfra(i) - dcfa_ini(i))
     952      qissr(i)   = MAX(0., qissr(i) - dqta_ini(i))
     953      cldfra(i)  = MAX(0., MIN(1., cldfra(i) + dcfa_ini(i)))
     954      qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqta_ini(i)))
     955      qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dcfa_ini(i) * qsat(i)))
     956      contfra(i) = MAX(0., MIN(cldfra(i), contfra(i) + dcfa_cir(i) + dcfa_ini(i)))
     957      qcont(i)   = MAX(0., MIN(qcld(i), qcont(i) + dqta_cir(i) + dqta_ini(i)))
     958
     959      !--Diagnostics
     960      dcfa_ini(i) = dcfa_ini(i) / dtime
     961      dqia_ini(i) = dqia_ini(i) / dtime
     962      dqta_ini(i) = dqta_ini(i) / dtime
     963      dcfa_sub(i) = dcfa_sub(i) / dtime
     964      dqia_sub(i) = dqta_sub(i) / dtime
     965      dqta_sub(i) = dqta_sub(i) / dtime
     966      dcfa_cir(i) = dcfa_cir(i) / dtime
     967      dqta_cir(i) = dqta_cir(i) / dtime
     968      dcfa_mix(i) = dcfa_mix(i) / dtime
     969      dqia_mix(i) = dqia_mix(i) / dtime
     970      dqta_mix(i) = dqta_mix(i) / dtime
     971
     972      !-------------------------------------------
     973      !--       FINAL BARRIERS AND OUTPUTS      --
     974      !-------------------------------------------
     975
     976      IF ( cldfra(i) .LT. eps ) THEN
     977        !--If the cloud is too small, it is sublimated.
     978        cldfra(i) = 0.
     979        contfra(i)= 0.
     980        qcld(i)   = 0.
     981        qvc(i)    = 0.
     982        qincld(i) = qsat(i)
     983      ELSE
     984        qincld(i) = qcld(i) / cldfra(i)
     985      ENDIF ! cldfra .LT. eps
     986
     987      IF ( contfra(i) .LT. eps ) THEN
     988        contfra(i) = 0.
     989        qcont(i) = 0.
    803990      ENDIF
    804991
    805     ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
    806 
    807   ENDIF   ! end keepgoing
    808 
    809 ENDDO     ! end loop on i
     992    ENDIF ! keepgoing
     993  ENDDO
     994ENDIF ! ok_plane_contrail
     995
    810996
    811997END SUBROUTINE condensation_ice_supersat
    812998!**********************************************************************************
    813999
    814 SUBROUTINE deposition_sublimation( &
    815     qvapincld, qiceincld, temp, qsat, pplay, dtime, &
    816     qvapincld_new)
    817 
    818 USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, RD, EPS_W
    819 USE lmdz_lscp_ini,        ONLY: eps
    820 USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
    821 
    822 REAL,     INTENT(IN) :: qvapincld     !
    823 REAL,     INTENT(IN) :: qiceincld     !
    824 REAL,     INTENT(IN) :: temp          !
    825 REAL,     INTENT(IN) :: qsat          !
    826 REAL,     INTENT(IN) :: pplay         !
    827 REAL,     INTENT(IN) :: dtime         ! time step [s]
    828 REAL,     INTENT(OUT):: qvapincld_new !
    829 
    830 !
    831 ! for unadjusted clouds
    832 REAL :: qice_ratio
     1000FUNCTION QVAPINCLD_DEPSUB( &
     1001    qvapincld, qiceincld, temp, qsat, pplay, dtime)
     1002
     1003USE lmdz_lscp_ini, ONLY: RV, RLSTT, RTT, RD, EPS_W
     1004USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
     1005USE lmdz_lscp_ini, ONLY: eps
     1006
     1007IMPLICIT NONE
     1008
     1009! inputs
     1010REAL :: qvapincld     !
     1011REAL :: qiceincld     !
     1012REAL :: temp          !
     1013REAL :: qsat          !
     1014REAL :: pplay         !
     1015REAL :: dtime         ! time step [s]
     1016! outpu
     1017REAL :: qvapincld_depsub
     1018
     1019!
     1020! local
    8331021REAL :: pres_sat, rho, kappa
    8341022REAL :: air_thermal_conduct, water_vapor_diff
    … …  
    9391127  !--If the cloud is initially supersaturated
    9401128  !--Exact explicit integration (qvc exact, qice explicit)
    941   qvapincld_new = qsat + ( qvapincld - qsat ) &
    942                 * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
     1129  qvapincld_depsub = qsat + ( qvapincld - qsat ) &
     1130                   * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
    9431131ELSE
    9441132  !--If the cloud is initially subsaturated
    945   !--Exact explicit integration (qice exact, qvc explicit)
    946   !--The barrier is set so that the resulting vapor in cloud
    947   !--cannot be greater than qsat
    948   !--qice_ratio is the ratio between the new ice content and
    949   !--the old one, it is comprised between 0 and 1
    950   qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) )
    951 
    952   IF ( qice_ratio .LT. 0. ) THEN
    953     !--The new vapor in cloud is set to 0 so that the
    954     !--sublimation process does not activate
    955     qvapincld_new = 0.
    956   ELSE
    957     !--Else, the sublimation process is activated with the
    958     !--diagnosed new cloud water vapor
    959     !--The new vapor in the cloud is increased with the
    960     !--sublimated ice
    961     qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 )
    962     !--The new vapor in the cloud cannot be greater than qsat
    963     qvapincld_new = MIN(qvapincld_new, qsat)
    964   ENDIF ! qice_ratio .LT. 0.
     1133  !--Exact explicit integration (qvc exact, qice explicit)
     1134  !--Same but depo_coef_cirrus = 1
     1135  qvapincld_depsub = qsat + ( qvapincld - qsat ) &
     1136                   * EXP( - dtime * qiceincld / kappa / rm_ice**2 )
    9651137ENDIF ! qvapincld .GT. qsat
    9661138
    967 END SUBROUTINE deposition_sublimation
     1139END FUNCTION QVAPINCLD_DEPSUB
    9681140
    9691141END MODULE lmdz_lscp_condensation

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

    r5536 r5601  
    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, rvc_ancien, rcont_ancien, radpas, radsol, rain_fall, ratqs, &
     25       cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, radpas, radsol, rain_fall, ratqs, &
    2626       rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
    2727       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
    … …  
    414414  IF ( ok_ice_supersat ) THEN
    415415    ancien_ok=ancien_ok.AND.phyetat0_get(cf_ancien,"CFANCIEN","CFANCIEN",0.)
    416     ancien_ok=ancien_ok.AND.phyetat0_get(rvc_ancien,"RVCANCIEN","RVCANCIEN",0.)
     416    ancien_ok=ancien_ok.AND.phyetat0_get(qvc_ancien,"QVCANCIEN","QVCANCIEN",0.)
    417417  ELSE
    418418    cf_ancien(:,:)=0.
    419     rvc_ancien(:,:)=0.
     419    qvc_ancien(:,:)=0.
    420420  ENDIF
    421421
    422422  ! cas specifique des variables de l'aviation
    423423  IF ( ok_plane_contrail ) THEN
    424     ancien_ok=ancien_ok.AND.phyetat0_get(rcont_ancien,"RCONTANCIEN","RCONTANCIEN",0.)
    425   ELSE
    426     rcont_ancien(:,:)=0.
     424    ancien_ok=ancien_ok.AND.phyetat0_get(cfa_ancien,"CFAANCIEN","CFAANCIEN",0.)
     425    ancien_ok=ancien_ok.AND.phyetat0_get(cfa_ancien,"QTAANCIEN","QTAANCIEN",0.)
     426  ELSE
     427    cfa_ancien(:,:)=0.
     428    qta_ancien(:,:)=0.
    427429  ENDIF
    428430
    … …  
    449451  IF ( ok_ice_supersat ) THEN
    450452    IF ( (maxval(cf_ancien).EQ.minval(cf_ancien))     .OR. &
    451          (maxval(rvc_ancien).EQ.minval(rvc_ancien)) ) THEN
     453         (maxval(qvc_ancien).EQ.minval(qvc_ancien)) ) THEN
    452454       ancien_ok=.false.
    453455     ENDIF
    … …  
    455457
    456458  IF ( ok_plane_contrail ) THEN
    457     IF ( maxval(rcont_ancien).EQ.minval(rcont_ancien) ) THEN
     459    IF ( ( maxval(cfa_ancien).EQ.minval(cfa_ancien) ) .OR. &
     460         ( maxval(qta_ancien).EQ.minval(qta_ancien) ) ) THEN
    458461       ancien_ok=.false.
    459462     ENDIF

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

    r5536 r5601  
    2323                                prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
    2424                                ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
    25                                 rvc_ancien, rcont_ancien, u_ancien, v_ancien,&
     25                                qvc_ancien, cfa_ancien, qta_ancien,          &
     26                                u_ancien, v_ancien,                          &
    2627                                clwcon, rnebcon, ratqs, pbl_tke,             &
    2728                                wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
    … …  
    253254    IF ( ok_ice_supersat ) THEN
    254255      CALL put_field(pass,"CFANCIEN", "CFANCIEN", cf_ancien)
    255       CALL put_field(pass,"RVCANCIEN", "RVCANCIEN", rvc_ancien)
     256      CALL put_field(pass,"QVCANCIEN", "QVCANCIEN", qvc_ancien)
    256257    ENDIF
    257258
    258259    IF ( ok_plane_contrail ) THEN
    259       CALL put_field(pass,"RCONTANCIEN", "RCONTANCIEN", rcont_ancien)
     260      CALL put_field(pass,"CFAANCIEN", "CFAANCIEN", cfa_ancien)
     261      CALL put_field(pass,"QTAANCIEN", "QTAANCIEN", qta_ancien)
    260262    ENDIF
    261263

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

    r5598 r5601  
    2323      REAL, SAVE, ALLOCATABLE :: u_seri(:,:), v_seri(:,:)
    2424      !$OMP THREADPRIVATE(u_seri, v_seri)
    25       REAL, SAVE, ALLOCATABLE :: cf_seri(:,:), rvc_seri(:,:)
    26       !$OMP THREADPRIVATE(cf_seri, rvc_seri)
     25      REAL, SAVE, ALLOCATABLE :: cf_seri(:,:), qvc_seri(:,:)
     26      !$OMP THREADPRIVATE(cf_seri, qvc_seri)
    2727      REAL, SAVE, ALLOCATABLE :: l_mixmin(:,:,:),l_mix(:,:,:),wprime(:,:,:)
    2828      !$OMP THREADPRIVATE(l_mixmin, l_mix, wprime)
    … …  
    4343      REAL, SAVE, ALLOCATABLE :: d_u_dyn(:,:), d_v_dyn(:,:)
    4444      !$OMP THREADPRIVATE(d_u_dyn, d_v_dyn)
    45       REAL, SAVE, ALLOCATABLE :: d_cf_dyn(:,:), d_rvc_dyn(:,:)
    46       !$OMP THREADPRIVATE(d_cf_dyn, d_rvc_dyn)
     45      REAL, SAVE, ALLOCATABLE :: d_cf_dyn(:,:), d_qvc_dyn(:,:)
     46      !$OMP THREADPRIVATE(d_cf_dyn, d_qvc_dyn)
    4747      REAL, SAVE, ALLOCATABLE :: d_tr_dyn(:,:,:)
    4848      !$OMP THREADPRIVATE(d_tr_dyn)
    … …  
    649649      REAL, SAVE, ALLOCATABLE :: ql_seri_lscp(:,:)
    650650      !$OMP THREADPRIVATE(ql_seri_lscp)
    651       REAL, SAVE, ALLOCATABLE :: ratio_ql_qtot(:,:)
    652       !$OMP THREADPRIVATE(ratio_ql_qtot)
    653651      REAL, SAVE, ALLOCATABLE :: qi_seri_lscp(:,:)
    654652      !$OMP THREADPRIVATE(qi_seri_lscp)
    655       REAL, SAVE, ALLOCATABLE :: ratio_qi_qtot(:,:)
    656       !$OMP THREADPRIVATE(ratio_qi_qtot)
    657653      REAL, SAVE, ALLOCATABLE :: dcf_sub(:,:), dcf_con(:,:), dcf_mix(:,:)
    658654      !$OMP THREADPRIVATE(dcf_sub, dcf_con, dcf_mix)
    … …  
    671667      REAL, SAVE, ALLOCATABLE :: d_q_avi(:,:)
    672668      !$OMP THREADPRIVATE(d_q_avi)
    673       REAL, SAVE, ALLOCATABLE :: rcont_seri(:,:), d_rcont_dyn(:,:)
    674       !$OMP THREADPRIVATE(rcont_seri, d_rcont_dyn)
     669      REAL, SAVE, ALLOCATABLE :: cfa_seri(:,:), d_cfa_dyn(:,:)
     670      !$OMP THREADPRIVATE(cfa_seri, d_cfa_dyn)
     671      REAL, SAVE, ALLOCATABLE :: qta_seri(:,:), d_qta_dyn(:,:)
     672      !$OMP THREADPRIVATE(qta_seri, d_qta_dyn)
    675673      REAL, SAVE, ALLOCATABLE :: flight_dist(:,:), flight_h2o(:,:)
    676674      !$OMP THREADPRIVATE(flight_dist, flight_h2o)
    … …  
    679677      REAL, SAVE, ALLOCATABLE :: potcontfraP(:,:), potcontfraNP(:,:)
    680678      !$OMP THREADPRIVATE(potcontfraP, potcontfraNP)
    681       REAL, SAVE, ALLOCATABLE :: contfra(:,:), radocond_cont(:,:), dcontfra_cir(:,:)
    682       !$OMP THREADPRIVATE(contfra, radocond_cont, dcontfra_cir)
    683       REAL, SAVE, ALLOCATABLE :: dcf_avi(:,:), dqi_avi(:,:), dqvc_avi(:,:)
    684       !$OMP THREADPRIVATE(dcf_avi, dqi_avi, dqvc_avi)
     679      REAL, SAVE, ALLOCATABLE :: qice_cont(:,:), dcfa_cir(:,:), dqta_cir(:,:)
     680      !$OMP THREADPRIVATE(qice_cont, dcfa_cir, dqta_cir)
     681      REAL, SAVE, ALLOCATABLE :: dcfa_ini(:,:), dqia_ini(:,:), dqta_ini(:,:)
     682      !$OMP THREADPRIVATE(dcfa_ini, dqia_ini, dqta_ini)
     683      REAL, SAVE, ALLOCATABLE :: dcfa_sub(:,:), dqia_sub(:,:), dqta_sub(:,:)
     684      !$OMP THREADPRIVATE(dcfa_sub, dqia_sub, dqta_sub)
     685      REAL, SAVE, ALLOCATABLE :: dcfa_mix(:,:), dqia_mix(:,:), dqta_mix(:,:)
     686      !$OMP THREADPRIVATE(dcfa_mix, dqia_mix, dqta_mix)
    685687      REAL, SAVE, ALLOCATABLE :: cldfra_nocont(:,:), cldtau_nocont(:,:), cldemi_nocont(:,:)
    686688      !$OMP THREADPRIVATE(cldfra_nocont, cldtau_nocont, cldemi_nocont)
    687       REAL, SAVE, ALLOCATABLE :: cldh_nocont(:), contcov(:)
    688       !$OMP THREADPRIVATE(cldh_nocont, contcov)
     689      REAL, SAVE, ALLOCATABLE :: cldh_nocont(:), contcov(:), conttau(:,:), contemi(:,:)
     690      !$OMP THREADPRIVATE(cldh_nocont, contcov, conttau, contemi)
    689691      REAL, SAVE, ALLOCATABLE :: fiwp_nocont(:), fiwc_nocont(:,:), ref_ice_nocont(:,:)
    690692      !$OMP THREADPRIVATE(fiwp_nocont, fiwc_nocont, ref_ice_nocont)
    … …  
    851853! SN
    852854      ALLOCATE(u_seri(klon,klev),v_seri(klon,klev))
    853       ALLOCATE(cf_seri(klon,klev),rvc_seri(klon,klev))
     855      ALLOCATE(cf_seri(klon,klev),qvc_seri(klon,klev))
    854856      ALLOCATE(l_mixmin(klon,klev+1,nbsrf),l_mix(klon,klev+1,nbsrf),wprime(klon,klev+1,nbsrf))
    855857      ALLOCATE(pbl_eps(klon,klev+1,nbsrf+1))
    … …  
    864866      ALLOCATE(d_q_dyn2d(klon),d_ql_dyn2d(klon),d_qs_dyn2d(klon), d_qbs_dyn2d(klon))
    865867      ALLOCATE(d_u_dyn(klon,klev),d_v_dyn(klon,klev))
    866       ALLOCATE(d_cf_dyn(klon,klev),d_rvc_dyn(klon,klev))
     868      ALLOCATE(d_cf_dyn(klon,klev),d_qvc_dyn(klon,klev))
    867869      ALLOCATE(d_tr_dyn(klon,klev,nbtr))                   !RomP
    868870      ALLOCATE(d_t_con(klon,klev),d_q_con(klon,klev),d_q_con_zmasse(klon,klev))
    … …  
    12311233      ALLOCATE(qsub(klon,klev), qissr(klon,klev), qcld(klon,klev))
    12321234      ALLOCATE(subfra(klon,klev), issrfra(klon,klev))
    1233       ALLOCATE(gamma_cond(klon,klev), ql_seri_lscp(klon,klev), ratio_ql_qtot(klon,klev))
    1234       ALLOCATE(qi_seri_lscp(klon,klev), ratio_qi_qtot(klon,klev))
     1235      ALLOCATE(gamma_cond(klon,klev), ql_seri_lscp(klon,klev), qi_seri_lscp(klon,klev))
    12351236      ALLOCATE(dcf_sub(klon,klev), dcf_con(klon,klev), dcf_mix(klon,klev))
    12361237      ALLOCATE(dqi_adj(klon,klev), dqi_sub(klon,klev), dqi_con(klon,klev), dqi_mix(klon,klev))
    … …  
    12421243!-- LSCP - aviation and contrails variables
    12431244      ALLOCATE(d_q_avi(klon,klev))
    1244       ALLOCATE(rcont_seri(klon,klev), d_rcont_dyn(klon,klev))
     1245      ALLOCATE(cfa_seri(klon,klev), d_cfa_dyn(klon,klev))
     1246      ALLOCATE(qta_seri(klon,klev), d_qta_dyn(klon,klev))
    12451247      ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
    12461248      ALLOCATE(Tcritcont(klon,klev), qcritcont(klon,klev))
    12471249      ALLOCATE(potcontfraP(klon,klev), potcontfraNP(klon,klev))
    1248       ALLOCATE(contfra(klon,klev), radocond_cont(klon,klev), dcontfra_cir(klon,klev))
    1249       ALLOCATE(dcf_avi(klon,klev), dqi_avi(klon,klev), dqvc_avi(klon,klev))
     1250      ALLOCATE(qice_cont(klon,klev), dcfa_cir(klon,klev), dqta_cir(klon,klev))
     1251      ALLOCATE(dcfa_ini(klon,klev), dqia_ini(klon,klev), dqta_ini(klon,klev))
     1252      ALLOCATE(dcfa_sub(klon,klev), dqia_sub(klon,klev), dqta_sub(klon,klev))
     1253      ALLOCATE(dcfa_mix(klon,klev), dqia_mix(klon,klev), dqta_mix(klon,klev))
    12501254      ALLOCATE(cldfra_nocont(klon,klev), cldtau_nocont(klon,klev), cldemi_nocont(klon,klev))
    1251       ALLOCATE(cldh_nocont(klon), contcov(klon))
     1255      ALLOCATE(cldh_nocont(klon), contcov(klon), conttau(klon,klev), contemi(klon,klev))
    12521256      ALLOCATE(fiwp_nocont(klon), fiwc_nocont(klon,klev), ref_ice_nocont(klon,klev))
    12531257      ALLOCATE(topsw_nocont(klon), solsw_nocont(klon))
    … …  
    13241328! SN
    13251329      DEALLOCATE(u_seri,v_seri)
    1326       DEALLOCATE(cf_seri,rvc_seri)
     1330      DEALLOCATE(cf_seri,qvc_seri)
    13271331      DEALLOCATE(l_mixmin,l_mix,wprime)
    13281332      DEALLOCATE(tke_shear,tke_buoy,tke_trans)
    … …  
    13341338      DEALLOCATE(d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, d_qbs_dyn2d)
    13351339      DEALLOCATE(d_u_dyn,d_v_dyn)
    1336       DEALLOCATE(d_cf_dyn,d_rvc_dyn)
     1340      DEALLOCATE(d_cf_dyn,d_qvc_dyn)
    13371341      DEALLOCATE(d_tr_dyn)                      !RomP
    13381342      DEALLOCATE(d_t_con,d_q_con,d_q_con_zmasse)
    … …  
    16481652      DEALLOCATE(qsub, qissr, qcld)
    16491653      DEALLOCATE(subfra, issrfra)
    1650       DEALLOCATE(gamma_cond, ql_seri_lscp, ratio_ql_qtot, qi_seri_lscp, ratio_qi_qtot)
     1654      DEALLOCATE(gamma_cond, ql_seri_lscp, qi_seri_lscp)
    16511655      DEALLOCATE(dcf_sub, dcf_con, dcf_mix)
    16521656      DEALLOCATE(dqi_adj, dqi_sub, dqi_con, dqi_mix)
    … …  
    16571661
    16581662!-- LSCP - aviation and contrails variables
    1659       DEALLOCATE(d_q_avi, rcont_seri, d_rcont_dyn, flight_dist, flight_h2o)
     1663      DEALLOCATE(d_q_avi, cfa_seri, d_cfa_dyn, qta_seri, d_qta_dyn, flight_dist, flight_h2o)
    16601664      DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
    1661       DEALLOCATE(contfra, radocond_cont, dcontfra_cir)
    1662       DEALLOCATE(dcf_avi, dqi_avi, dqvc_avi)
    1663       DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont)
     1665      DEALLOCATE(qice_cont, dcfa_cir, dqta_cir, dcfa_ini, dqia_ini, dqta_ini)
     1666      DEALLOCATE(dcfa_sub, dqia_sub, dqta_sub, dcfa_mix, dqia_mix, dqta_mix)
     1667      DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont, conttau, contemi)
    16641668      DEALLOCATE(cldh_nocont, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)
    16651669      DEALLOCATE(topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont)

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

    r5598 r5601  
    21302130  TYPE(ctrl_out), SAVE :: o_dcfdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    21312131    'dcfdyn', 'Dynamics cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
    2132   TYPE(ctrl_out), SAVE :: o_rvcseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2133     'rvcseri', 'Cloudy water vapor to total water vapor ratio', '-', (/ ('', i=1, 10) /))
    2134   TYPE(ctrl_out), SAVE :: o_drvcdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2135     'drvcdyn', 'Dynamics cloudy water vapor to total water vapor ratio tendency', 's-1', (/ ('', i=1, 10) /))
     2132  TYPE(ctrl_out), SAVE :: o_qvcseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2133    'qvcseri', 'Cloudy water vapor to total water vapor ratio', '-', (/ ('', i=1, 10) /))
     2134  TYPE(ctrl_out), SAVE :: o_dqvcdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2135    'dqvcdyn', 'Dynamics cloudy water vapor to total water vapor ratio tendency', 's-1', (/ ('', i=1, 10) /))
    21362136  TYPE(ctrl_out), SAVE :: o_qsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    21372137    'qsub', 'Subsaturated clear sky total water', 'kg/kg', (/ ('', i=1, 10) /))
    … …  
    21882188  TYPE(ctrl_out), SAVE :: o_dqavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    21892189    'dqavi', 'Water vapor emissions from aviation tendency', 'kg/kg/s', (/ ('',i=1,10) /))
    2190   TYPE(ctrl_out), SAVE :: o_rcontseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2191     'rcontseri', 'Contrails fraction to total cloud fraction ratio', '-', (/ ('',i=1,10) /))
    2192   TYPE(ctrl_out), SAVE :: o_drcontdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2193     'drcontdyn', 'Dynamics contrails fraction to total cloud fraction ratio tendency', 's-1', (/ ('',i=1,10) /))
     2190  TYPE(ctrl_out), SAVE :: o_cfaseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2191    'cfaseri', 'Linear contrail fraction', '-', (/ ('',i=1,10) /))
     2192  TYPE(ctrl_out), SAVE :: o_dcfadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2193    'dcfadyn', 'Dynamics linear contrail fraction tendency', 's-1', (/ ('',i=1,10) /))
     2194  TYPE(ctrl_out), SAVE :: o_qtaseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2195    'qtaseri', 'Linear contrail total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
     2196  TYPE(ctrl_out), SAVE :: o_dqtadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2197    'dqtadyn', 'Dynamics linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
    21942198  TYPE(ctrl_out), SAVE :: o_Tcritcont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    21952199    'Tcritcont', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
    … …  
    22002204  TYPE(ctrl_out), SAVE :: o_potcontfraNP = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    22012205    'potcontfraNP', 'Potential non-persistent contrail fraction', '-', (/ ('', i=1,10)/))
    2202   TYPE(ctrl_out), SAVE :: o_contfra = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2203     'contfra', 'Linear contrail fraction', '-', (/ ('', i=1,10)/))
    22042206  TYPE(ctrl_out), SAVE :: o_qice_cont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2205     'qice_cont', 'Linear contrails ice specific humidity seen by radiation', 'kg/kg', (/ ('', i=1, 10) /))
    2206   TYPE(ctrl_out), SAVE :: o_dcontfracir = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    2207     'dcontfracir', 'Linear contrail fraction to cirrus cloud fraction tendency', '-', (/ ('', i=1,10)/))
    2208   TYPE(ctrl_out), SAVE :: o_dcfavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2209     'dcfavi', 'Aviation cloud fraction tendency', 's-1', (/ ('', i=1, 10) /))
    2210   TYPE(ctrl_out), SAVE :: o_dqiavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2211     'dqiavi', 'Aviation ice tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    2212   TYPE(ctrl_out), SAVE :: o_dqvcavi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2213     'dqvcavi', 'Aviation cloudy water vapor tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2207    'qice_cont', 'Linear contrails ice specific humidity', 'kg/kg', (/ ('', i=1, 10) /))
     2208  TYPE(ctrl_out), SAVE :: o_dcfaini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2209    'dcfaini', 'Initial formation linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2210  TYPE(ctrl_out), SAVE :: o_dqiaini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2211    'dqiaini', 'Initial formation linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2212  TYPE(ctrl_out), SAVE :: o_dqtaini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2213    'dqtaini', 'Initial formation linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2214  TYPE(ctrl_out), SAVE :: o_dcfacir = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2215    'dcfacir', 'Conversion to cirrus linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2216  TYPE(ctrl_out), SAVE :: o_dqtacir = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2217    'dqtacir', 'Conversion to cirrus linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2218  TYPE(ctrl_out), SAVE :: o_dcfasub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2219    'dcfasub', 'Sublimation linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2220  TYPE(ctrl_out), SAVE :: o_dqiasub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2221    'dqiasub', 'Sublimation linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2222  TYPE(ctrl_out), SAVE :: o_dqtasub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2223    'dqtasub', 'Sublimation linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2224  TYPE(ctrl_out), SAVE :: o_dcfamix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2225    'dcfamix', 'Mixing linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2226  TYPE(ctrl_out), SAVE :: o_dqiamix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2227    'dqiamix', 'Mixing linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2228  TYPE(ctrl_out), SAVE :: o_dqtamix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2229    'dqtamix', 'Mixing linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    22142230  TYPE(ctrl_out), SAVE :: o_flight_dist = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    22152231    'flightdist', 'Aviation flown distance', 'm/s/m^3', (/ ('', i=1, 10) /))
    … …  
    22262242  TYPE(ctrl_out), SAVE :: o_contcov = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    22272243    'contcov', 'Total contrails cover', '-', (/ ('', i=1, 10) /))
     2244  TYPE(ctrl_out), SAVE :: o_conttau = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2245    'conttau', 'Contrails optical thickness', '1', (/ ('', i=1, 10) /))
     2246  TYPE(ctrl_out), SAVE :: o_contemi = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2247    'contemi', 'Contrails optical emissivity', '1', (/ ('', i=1, 10) /))
    22282248  TYPE(ctrl_out), SAVE :: o_iwp_nocont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    22292249    'iwp_nocont', 'Cloud ice water path w/o contrails', 'kg/m2', (/ ('', i=1, 10) /))

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

    r5598 r5601  
    222222         o_col_O3_strato, o_col_O3_tropo,                 &
    223223!-- LSCP - condensation and ice supersaturation variables
    224          o_cfseri, o_dcfdyn, o_rvcseri, o_drvcdyn, &
     224         o_cfseri, o_dcfdyn, o_qvcseri, o_dqvcdyn, &
    225225         o_qsub, o_qissr, o_qcld, o_subfra, o_issrfra, o_gammacond, &
    226226         o_dcfsub, o_dcfcon, o_dcfmix, o_dqiadj, o_dqisub, o_dqicon, o_dqimix, &
    … …  
    229229         o_issrfra250to300, o_issrfra300to400, o_issrfra400to500, &
    230230!-- LSCP - aviation variables
    231          o_rcontseri, o_drcontdyn, o_dqavi, o_contfra, o_qice_cont, &
     231         o_cfaseri, o_dcfadyn, o_qtaseri, o_dqtadyn, o_dqavi, &
    232232         o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
    233          o_dcontfracir, o_dcfavi, o_dqiavi, o_dqvcavi, o_flight_dist, o_flight_h2o, &
     233         o_flight_dist, o_flight_h2o, o_qice_cont, o_dcfacir, o_dqtacir, &
     234         o_dcfaini, o_dqiaini, o_dqtaini, o_dcfasub, o_dqiasub, o_dqtasub, &
     235         o_dcfamix, o_dqiamix, o_dqtamix, &
    234236         o_cldfra_nocont, o_cldtau_nocont, o_cldemi_nocont, o_cldh_nocont, &
    235          o_contcov, o_iwp_nocont, o_iwc_nocont, o_ref_ice_nocont, &
     237         o_contcov, o_conttau, o_contemi, o_iwp_nocont, o_iwc_nocont, o_ref_ice_nocont, &
    236238         o_tops_nocont, o_topl_nocont, o_sols_nocont, o_soll_nocont, o_nettop_nocont, &
    237239!--interactive CO2
    … …  
    350352         wdtrainA, wdtrainS, wdtrainM, n2, s2, strig, zcong, zlcl_th, proba_notrig, &
    351353         random_notrig, &
    352          cf_seri, d_cf_dyn, rvc_seri, d_rvc_dyn, &
     354         cf_seri, d_cf_dyn, qvc_seri, d_qvc_dyn, &
    353355         qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
    354356         dcf_sub, dcf_con, dcf_mix, &
    … …  
    358360         issrfra100to150, issrfra150to200, issrfra200to250, &
    359361         issrfra250to300, issrfra300to400, issrfra400to500, &
    360          rcont_seri, d_rcont_dyn, d_q_avi, contfra, radocond_cont, &
     362         cfa_seri, d_cfa_dyn, qta_seri, d_qta_dyn, d_q_avi, &
    361363         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    362          dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
     364         flight_dist, flight_h2o, qice_cont, dcfa_cir, dqta_cir, &
     365         dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub, &
     366         dcfa_mix, dqia_mix, dqta_mix, &
    363367         cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
    364          contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
     368         contcov, conttau, contemi, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
    365369         topsw_nocont, toplw_nocont, solsw_nocont, sollw_nocont, &
    366370         alp_bl_det, alp_bl_fluct_m, alp_bl_conv, &
    … …  
    21182122         CALL histwrite_phy(o_cfseri, cf_seri)
    21192123         CALL histwrite_phy(o_dcfdyn, d_cf_dyn)
    2120          CALL histwrite_phy(o_rvcseri, rvc_seri)
    2121          CALL histwrite_phy(o_drvcdyn, d_rvc_dyn)
     2124         CALL histwrite_phy(o_qvcseri, qvc_seri)
     2125         CALL histwrite_phy(o_dqvcdyn, d_qvc_dyn)
    21222126         CALL histwrite_phy(o_qsub, qsub)
    21232127         CALL histwrite_phy(o_qissr, qissr)
    … …  
    21482152!-- LSCP - aviation variables
    21492153       IF (ok_plane_contrail) THEN
    2150          CALL histwrite_phy(o_rcontseri, rcont_seri)
    2151          CALL histwrite_phy(o_drcontdyn, d_rcont_dyn)
     2154         CALL histwrite_phy(o_cfaseri, cfa_seri)
     2155         CALL histwrite_phy(o_dcfadyn, d_cfa_dyn)
     2156         CALL histwrite_phy(o_qtaseri, qta_seri)
     2157         CALL histwrite_phy(o_dqtadyn, d_qta_dyn)
    21522158         CALL histwrite_phy(o_flight_dist, flight_dist)
    21532159         CALL histwrite_phy(o_Tcritcont, Tcritcont)
    … …  
    21552161         CALL histwrite_phy(o_potcontfraP, potcontfraP)
    21562162         CALL histwrite_phy(o_potcontfraNP, potcontfraNP)
    2157          CALL histwrite_phy(o_contfra, contfra)
    2158          CALL histwrite_phy(o_qice_cont, radocond_cont)
    2159          CALL histwrite_phy(o_dcontfracir, dcontfra_cir)
    2160          CALL histwrite_phy(o_dcfavi, dcf_avi)
    2161          CALL histwrite_phy(o_dqiavi, dqi_avi)
    2162          CALL histwrite_phy(o_dqvcavi, dqvc_avi)
     2163         CALL histwrite_phy(o_qice_cont, qice_cont)
     2164         CALL histwrite_phy(o_dcfacir, dcfa_cir)
     2165         CALL histwrite_phy(o_dqtacir, dqta_cir)
     2166         CALL histwrite_phy(o_dcfaini, dcfa_ini)
     2167         CALL histwrite_phy(o_dqiaini, dqia_ini)
     2168         CALL histwrite_phy(o_dqtaini, dqta_ini)
     2169         CALL histwrite_phy(o_dcfasub, dcfa_sub)
     2170         CALL histwrite_phy(o_dqiasub, dqia_sub)
     2171         CALL histwrite_phy(o_dqtasub, dqta_sub)
     2172         CALL histwrite_phy(o_dcfamix, dcfa_mix)
     2173         CALL histwrite_phy(o_dqiamix, dqia_mix)
     2174         CALL histwrite_phy(o_dqtamix, dqta_mix)
    21632175         CALL histwrite_phy(o_cldfra_nocont, cldfra_nocont)
    21642176         CALL histwrite_phy(o_cldtau_nocont, cldtau_nocont)
    … …  
    21662178         CALL histwrite_phy(o_cldh_nocont, cldh_nocont)
    21672179         CALL histwrite_phy(o_contcov, contcov)
     2180         CALL histwrite_phy(o_conttau, conttau)
     2181         CALL histwrite_phy(o_contemi, contemi)
    21682182         CALL histwrite_phy(o_iwp_nocont, fiwp_nocont)
    21692183         CALL histwrite_phy(o_iwc_nocont, fiwc_nocont)

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

    r5575 r5601  
    9292      REAL, ALLOCATABLE, SAVE :: u_ancien(:,:), v_ancien(:,:)
    9393!$OMP THREADPRIVATE(u_ancien, v_ancien)
    94       REAL, ALLOCATABLE, SAVE :: cf_ancien(:,:), rvc_ancien(:,:), rcont_ancien(:,:)
    95 !$OMP THREADPRIVATE(cf_ancien, rvc_ancien, rcont_ancien)
     94      REAL, ALLOCATABLE, SAVE :: cf_ancien(:,:), qvc_ancien(:,:)
     95!$OMP THREADPRIVATE(cf_ancien, qvc_ancien)
     96      REAL, ALLOCATABLE, SAVE :: cfa_ancien(:,:), qta_ancien(:,:)
     97!$OMP THREADPRIVATE(cfa_ancien, qta_ancien)
    9698!!! RomP >>>
    9799      REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
    … …  
    587589      ALLOCATE(prw_ancien(klon), prlw_ancien(klon), prsw_ancien(klon), prbsw_ancien(klon))
    588590      ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
    589       ALLOCATE(cf_ancien(klon,klev), rvc_ancien(klon,klev), rcont_ancien(klon,klev))
     591      ALLOCATE(cf_ancien(klon,klev), qvc_ancien(klon,klev))
     592      ALLOCATE(cfa_ancien(klon,klev), qta_ancien(klon,klev))
    590593!!! Rom P >>>
    591594      ALLOCATE(tr_ancien(klon,klev,nbtr))
    … …  
    814817      DEALLOCATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
    815818      DEALLOCATE(u_ancien, v_ancien)
    816       DEALLOCATE(cf_ancien, rvc_ancien, rcont_ancien)
     819      DEALLOCATE(cf_ancien, qvc_ancien, cfa_ancien, qta_ancien)
    817820      DEALLOCATE(tr_ancien)                           !RomP
    818821      DEALLOCATE(ratqs, pbl_tke,coefh,coefm)

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

    r5598 r5601  
    146146       ! Variables locales pour effectuer les appels en serie
    147147       t_seri,q_seri,ql_seri,qs_seri,qbs_seri, &
    148        u_seri,v_seri,cf_seri,rvc_seri,tr_seri, &
     148       u_seri,v_seri,cf_seri,qvc_seri,tr_seri, &
    149149       rhcl, &
    150150       ! Dynamic tendencies (diagnostics)
    151151       d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_qbs_dyn, &
    152        d_u_dyn,d_v_dyn,d_cf_dyn,d_rvc_dyn,d_tr_dyn, &
     152       d_u_dyn,d_v_dyn,d_cf_dyn,d_qvc_dyn,d_tr_dyn, &
    153153       d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d,d_qbs_dyn2d, &
    154154       ! Physic tendencies
    … …  
    324324       !-- LSCP - condensation and ice supersaturation variables
    325325       qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
    326        ql_seri_lscp, ratio_ql_qtot, qi_seri_lscp, ratio_qi_qtot, &
     326       ql_seri_lscp, qi_seri_lscp, &
    327327       dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
    328328       dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    329329       !-- LSCP - aviation and contrails variables
    330        d_q_avi, rcont_seri, d_rcont_dyn, flight_dist, flight_h2o, &
    331        contfra, radocond_cont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    332        dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, &
     330       cfa_seri, qta_seri, d_cfa_dyn, d_qta_dyn, &
     331       d_q_avi, flight_dist, flight_h2o, &
     332       qice_cont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    333333       cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
    334        contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
     334       conttau, contemi, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
    335335       topsw_nocont, solsw_nocont, toplw_nocont, sollw_nocont, &
    336336       !
    … …  
    518518    !
    519519    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
    520     INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, irvc, ircont
    521 !$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, irvc, ircont)
     520    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, iqvc, icfa, iqta
     521!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, icfa, iqta)
    522522    !
    523523    !
    … …  
    13641364       ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
    13651365       icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
    1366        irvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
    1367        ircont = strIdx(tracers(:)%name, addPhase('H2O', 'a'))
     1366       iqvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
     1367       icfa = strIdx(tracers(:)%name, addPhase('H2O', 'a'))
     1368       iqta = strIdx(tracers(:)%name, addPhase('H2O', 'q'))
    13681369!       CALL init_etat0_limit_unstruct
    13691370!       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
    … …  
    23582359                    sollwdown(:))
    23592360
    2360       !--Init for LSCP - condensation
    2361       ratio_ql_qtot(:,:) = 0.
    2362       ratio_qi_qtot(:,:) = 0.
    2363 
    2364 
    23652361    ENDIF
    23662362    !
    … …  
    24602456          qbs_seri(i,k)= 0.
    24612457          cf_seri(i,k) = 0.
    2462           rvc_seri(i,k)= 0.
    2463           rcont_seri(i,k)= 0.
     2458          qvc_seri(i,k)= 0.
     2459          cfa_seri(i,k)= 0.
     2460          qta_seri(i,k)= 0.
    24642461          !CR: ATTENTION, on rajoute la variable glace
    24652462          IF (nqo.EQ.2) THEN             !--vapour and liquid only
    … …  
    24712468             IF (ok_ice_supersat) THEN
    24722469               cf_seri(i,k) = qx(i,k,icf)
    2473                rvc_seri(i,k) = qx(i,k,irvc)
     2470               qvc_seri(i,k) = qx(i,k,iqvc)
    24742471             ENDIF
    24752472             IF (ok_plane_contrail) THEN
    2476                rcont_seri(i,k) = qx(i,k,ircont)
     2473               cfa_seri(i,k) = qx(i,k,icfa)
     2474               qta_seri(i,k) = qx(i,k,iqta)
    24772475             ENDIF
    24782476             IF (ok_bs) THEN
    … …  
    25522550       d_qbs_dyn(:,:)= (qbs_seri(:,:)-qbs_ancien(:,:))/phys_tstep
    25532551       d_cf_dyn(:,:) = (cf_seri(:,:)-cf_ancien(:,:))/phys_tstep
    2554        d_rvc_dyn(:,:)= (rvc_seri(:,:)-rvc_ancien(:,:))/phys_tstep
    2555        d_rcont_dyn(:,:)= (rcont_seri(:,:)-rcont_ancien(:,:))/phys_tstep
     2552       d_qvc_dyn(:,:)= (qvc_seri(:,:)-qvc_ancien(:,:))/phys_tstep
     2553       d_cfa_dyn(:,:)= (cfa_seri(:,:)-cfa_ancien(:,:))/phys_tstep
     2554       d_qta_dyn(:,:)= (qta_seri(:,:)-qta_ancien(:,:))/phys_tstep
    25562555       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
    25572556       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
    … …  
    25742573       d_qbs_dyn(:,:)= 0.0
    25752574       d_cf_dyn(:,:) = 0.0
    2576        d_rvc_dyn(:,:)= 0.0
    2577        d_rcont_dyn(:,:)= 0.0
     2575       d_qvc_dyn(:,:)= 0.0
     2576       d_cfa_dyn(:,:)= 0.0
     2577       d_qta_dyn(:,:)= 0.0
    25782578       d_q_dyn2d(:)  = 0.0
    25792579       d_ql_dyn2d(:) = 0.0
    … …  
    26932693
    26942694    !-- Needed for LSCP - condensation and ice supersaturation
    2695     IF (ok_ice_supersat) THEN
     2695    IF (ok_plane_contrail.AND.ok_ice_supersat) THEN
    26962696      DO k = 1, klev
    26972697        DO i = 1, klon
    26982698          IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
    2699             ratio_ql_qtot(i,k) = ql_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    2700             ratio_qi_qtot(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    2701             rvc_seri(i,k) = rvc_seri(i,k) * q_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2699            ql_seri_lscp(i,k) = ql_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2700            qi_seri_lscp(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2701            qvc_seri(i,k) = qvc_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2702            qta_seri(i,k) = qta_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27022703          ELSE
    2703             ratio_ql_qtot(i,k) = 0.
    2704             ratio_qi_qtot(i,k) = 0.
    2705             rvc_seri(i,k) = 0.
     2704            ql_seri_lscp(i,k) = 0.
     2705            qi_seri_lscp(i,k) = 0.
     2706            qvc_seri(i,k) = 0.
     2707            qta_seri(i,k) = 0.
     2708          ENDIF
     2709        ENDDO
     2710      ENDDO
     2711    ELSEIF (ok_ice_supersat) THEN
     2712      DO k = 1, klev
     2713        DO i = 1, klon
     2714          IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
     2715            ql_seri_lscp(i,k) = ql_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2716            qi_seri_lscp(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2717            qvc_seri(i,k) = qvc_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2718          ELSE
     2719            ql_seri_lscp(i,k) = 0.
     2720            qi_seri_lscp(i,k) = 0.
     2721            qvc_seri(i,k) = 0.
    27062722          ENDIF
    27072723        ENDDO
    … …  
    27112727        DO i = 1, klon
    27122728          IF ( ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) ) .GT. 0. ) THEN
    2713             ratio_ql_qtot(i,k) = ql_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    2714             ratio_qi_qtot(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2729            ql_seri_lscp(i,k) = ql_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2730            qi_seri_lscp(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27152731          ELSE
    2716             ratio_ql_qtot(i,k) = 0.
    2717             ratio_qi_qtot(i,k) = 0.
     2732            ql_seri_lscp(i,k) = 0.
     2733            qi_seri_lscp(i,k) = 0.
    27182734          ENDIF
    27192735        ENDDO
    … …  
    39073923    IF (ok_plane_contrail) CALL vertical_interpolation_aviation(klon, klev, paprs, pplay, t_seri, flight_dist, flight_h2o)
    39083924
     3925    IF (ok_ice_supersat) THEN
     3926      DO k = 1, klev
     3927        DO i = 1, klon
     3928          qvc_seri(i,k) = qvc_seri(i,k) * q_seri(i,k)
     3929        ENDDO
     3930      ENDDO
     3931    ENDIF
     3932    IF (ok_plane_contrail) THEN
     3933      DO k = 1, klev
     3934        DO i = 1, klon
     3935          qta_seri(i,k) = qta_seri(i,k) * q_seri(i,k)
     3936        ENDDO
     3937      ENDDO
     3938    ENDIF
    39093939    DO k = 1, klev
    39103940      DO i = 1, klon
    3911         ql_seri_lscp(i,k) = ratio_ql_qtot(i,k) * q_seri(i,k)
    3912         qi_seri_lscp(i,k) = ratio_qi_qtot(i,k) * q_seri(i,k)
     3941        ql_seri_lscp(i,k) = ql_seri_lscp(i,k) * q_seri(i,k)
     3942        qi_seri_lscp(i,k) = qi_seri_lscp(i,k) * q_seri(i,k)
    39133943      ENDDO
    39143944    ENDDO
    … …  
    39253955         pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave), &
    39263956         cell_area, stratomask, &
    3927          cf_seri, rvc_seri, u_seri, v_seri, &
     3957         cf_seri, qvc_seri, u_seri, v_seri, &
    39283958         qsub, qissr, qcld, subfra, issrfra, gamma_cond,  &
    39293959         dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
    39303960         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    3931          rcont_seri, flight_dist, flight_h2o, &
    3932          contfra, radocond_cont, Tcritcont, qcritcont, potcontfraP, &
    3933          potcontfraNP, dcontfra_cir, dcf_avi, dqi_avi, dqvc_avi, &
     3961         cfa_seri, qta_seri, flight_dist, flight_h2o, &
     3962         qice_cont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    39343963         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
    39353964         qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
    … …  
    44834512               zfice, dNovrN, ptconv, rnebcon, clwcon, &
    44844513               !--AB contrails
    4485                contfra, radocond_cont, cldfra_nocont, cldtau_nocont, cldemi_nocont, &
    4486                cldh_nocont, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)
     4514               cfa_seri, qice_cont, cldfra_nocont, cldtau_nocont, cldemi_nocont, &
     4515               conttau, contemi, cldh_nocont, contcov, &
     4516               fiwp_nocont, fiwc_nocont, ref_ice_nocont)
    44874517
    44884518       !
    … …  
    56815711          IF (nqo.ge.5 .and. ok_ice_supersat) THEN
    56825712             d_qx(i,k,icf) = ( cf_seri(i,k) - qx(i,k,icf) ) / phys_tstep
    5683              d_qx(i,k,irvc) = ( rvc_seri(i,k) - qx(i,k,irvc) ) / phys_tstep
     5713             d_qx(i,k,iqvc) = ( qvc_seri(i,k) - qx(i,k,iqvc) ) / phys_tstep
    56845714             IF (nqo.ge.6 .and. ok_plane_contrail) THEN
    5685                d_qx(i,k,ircont) = ( rcont_seri(i,k) - qx(i,k,ircont) ) / phys_tstep
     5715               d_qx(i,k,icfa) = ( cfa_seri(i,k) - qx(i,k,icfa) ) / phys_tstep
     5716               d_qx(i,k,iqta) = ( qta_seri(i,k) - qx(i,k,iqta) ) / phys_tstep
    56865717             ENDIF
    56875718          ENDIF
    … …  
    57205751    qbs_ancien(:,:)= qbs_seri(:,:)
    57215752    cf_ancien(:,:) = cf_seri(:,:)
    5722     rvc_ancien(:,:)= rvc_seri(:,:)
    5723     rcont_ancien(:,:)= rcont_seri(:,:)
     5753    qvc_ancien(:,:)= qvc_seri(:,:)
     5754    cfa_ancien(:,:)= cfa_seri(:,:)
     5755    qta_ancien(:,:)= qta_seri(:,:)
    57245756    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
    57255757    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
Note: See TracChangeset for help on using the changeset viewer.