Changeset 5609


Ignore:
Timestamp:
Apr 13, 2025, 7:10:19 PM (8 weeks ago)
Author:
aborella
Message:
  • changed treatment of prognostic variables for prognostic clouds
  • adapted sedimentation and autoconversion for prognostic cirrus clouds
  • cloud mixing, ice sedimentation and ISSR diagnosis are now consistent with the water vapor PDF
  • simplified assumptions for ice crystals deposition / sublimation
  • first version of the coupling between prognostic cirrus clouds and deep convection
  • added persistent contrail cirrus clouds in radiative diagnostics
Location:
LMDZ6/branches/contrails
Files:
23 edited

Legend:

Unmodified
Added
Removed

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

    r5601 r5609  
    678678        <field id="tke_max_sic"    long_name="Max Turb. Kinetic Energy sic"    unit="m2/s2" operation="maximum"/>
    679679       
    680         <field id="qrain_lsc"    long_name="LS specific rain content"    unit="kg/kg" />
    681         <field id="qsnow_lsc"    long_name="LS specific snow content"    unit="kg/kg" />
    682         <field id="dqreva"    long_name="LS rain tendency due to evaporation"    unit="kg/kg/s" />
    683         <field id="dqssub"    long_name="LS snow tendency due to sublimation"    unit="kg/kg/s" />
    684         <field id="dqrauto"    long_name="LS rain tendency due to autoconversion"    unit="kg/kg/s" />
    685         <field id="dqrcol"    long_name="LS rain tendency due to collection"    unit="kg/kg/s" />
    686         <field id="dqrmelt"    long_name="LS rain tendency due to melting"    unit="kg/kg/s" />
    687         <field id="dqrfreez"    long_name="LS rain tendency due to freezing"    unit="kg/kg/s" />
    688         <field id="dqsauto"    long_name="LS snow tendency due to autoconversion"    unit="kg/kg/s" />
    689         <field id="dqsagg"    long_name="LS snow tendency due to aggregation"    unit="kg/kg/s" />
    690         <field id="dqsrim"    long_name="LS snow tendency due to riming"    unit="kg/kg/s" />
    691         <field id="dqsmelt"    long_name="LS snow tendency due to melting"    unit="kg/kg/s" />
    692         <field id="dqsfreez"    long_name="LS snow tendency due to freezing"    unit="kg/kg/s" />
    693 
    694680        <field id="l_mix_ter"    long_name="PBL mixing length ter"    unit="m" />
    695681        <field id="l_mix_lic"    long_name="PBL mixing length lic"    unit="m" />
    … …  
    909895        <field id="dcfcon"     long_name="Condensation cloud fraction tendency"    unit="s-1" />
    910896        <field id="dcfmix"     long_name="Cloud mixing cloud fraction tendency"    unit="s-1" />
     897        <field id="dcfsed"     long_name="Sedimentation cloud fraction tendency"    unit="s-1" />
    911898        <field id="dqiadj"     long_name="Temperature adjustment ice tendency"    unit="kg/kg/s" />
    912899        <field id="dqisub"     long_name="Sublimation ice tendency"    unit="kg/kg/s" />
    913900        <field id="dqicon"     long_name="Condensation ice tendency"    unit="kg/kg/s" />
    914901        <field id="dqimix"     long_name="Cloud mixing ice tendency"    unit="kg/kg/s" />
     902        <field id="dqised"     long_name="Sedimentation ice tendency"    unit="kg/kg/s" />
    915903        <field id="dqvcadj"    long_name="Temperature adjustment cloudy water vapor tendency"    unit="kg/kg/s" />
    916904        <field id="dqvcsub"    long_name="Sublimation cloudy water vapor tendency"    unit="kg/kg/s" />
    917905        <field id="dqvccon"    long_name="Condensation cloudy water vapor tendency"    unit="kg/kg/s" />
    918906        <field id="dqvcmix"    long_name="Cloud mixing cloudy water vapor tendency"    unit="kg/kg/s" />
     907        <field id="dqvcsed"    long_name="Sedimentation cloudy water vapor tendency"    unit="kg/kg/s" />
    919908        <field id="qsatl"      long_name="Saturation with respect to liquid"    unit="kg/kg" />
    920909        <field id="qsati"      long_name="Saturation with respect to ice"    unit="kg/kg" />
    921910
    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" />
    926             <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
    927             <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
     911        <field id="cfaseri"    long_name="Linear contrail fraction" unit="-" />
     912        <field id="dcfadyn"    long_name="Dynamics linear contrail fraction tendency" unit="s-1" />
     913        <field id="pcfseri"    long_name="Contrail induced cirrus fraction" unit="-" />
     914        <field id="dpcfdyn"    long_name="Dynamics contrail induced cirrus fraction tendency" unit="s-1" />
     915        <field id="qvaseri"    long_name="Linear contrail vapor specific humidity" unit="kg/kg" />
     916        <field id="dqvadyn"    long_name="Dynamics linear contrail vapor specific humidity tendency" unit="kg/kg/s" />
     917        <field id="qiaseri"    long_name="Linear contrail ice specific humidity" unit="kg/kg" />
     918        <field id="dqiadyn"    long_name="Dynamics linear contrail ice specific humidity tendency" unit="kg/kg/s" />
     919        <field id="Tcritcont"  long_name="Temperature threshold for contrail formation"    unit="K" />
     920        <field id="qcritcont"  long_name="Specific humidity threshold for contrail formation"    unit="kg/kg" />
    928921        <field id="potcontfraP"  long_name="Fraction with pontential persistent contrail"    unit="-" />
    929922        <field id="potcontfraNP" long_name="Fraction with potential non-persistent contrail"    unit="-" />
    930         <field id="qice_cont"  long_name="Contrails ice specific humidity"    unit="kg/kg" />
     923        <field id="qice_perscont" long_name="Contrail induced cirrus ice specific humidity"    unit="kg/kg" />
    931924        <field id="dcfaini"    long_name="Initial formation linear contrail fraction tendency"    unit="s-1" />
    932925        <field id="dqiaini"    long_name="Initial formation linear contrail ice specific humidity tendency"    unit="kg/kg/s" />

  • LMDZ6/branches/contrails/libf/dyn3d_common/infotrac.f90

    r5536 r5609  
    266266   !---------------------------------------------------------------------------------------------------------------------------
    267267   nqtrue = SIZE(tracers)                                                                               !--- "true" tracers
    268    nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%name)     == 'H2O')     !--- Water phases
    269    nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%gen0Name) == 'H2O')     !--- Passed to phytrac
     268!   nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%name)     == 'H2O')     !--- Water phases
     269   nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. &
     270       ((delPhase(tracers(:)%name)     == 'H2O') .OR. &    !--- Water phases
     271        (delPhase(tracers(:)%name)     == 'CLDFRA')))
     272!   nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%gen0Name) == 'H2O')     !--- Passed to phytrac
     273   nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. &
     274       ((delPhase(tracers(:)%gen0Name)     == 'H2O') .OR. &    !--- Passed to phytrac
     275        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
    270276   nqCO2  =      COUNT( [type_trac == 'inco', type_trac == 'co2i'] )
    271277IF (CPPKEY_INCA) THEN
    … …  
    335341      t1%iadv       = iad
    336342      t1%isAdvected = iad >= 0
    337       t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' .OR. t1%component /= 'lmdz' !=== MORE EXCEPTIONS ? CO2i, SURSAT CLOUD H2O
     343!      t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' .OR. t1%component /= 'lmdz' !=== MORE EXCEPTIONS ? CO2i, SURSAT CLOUD H2O
     344      t1%isInPhysics=((delPhase(t1%gen0Name) /= 'H2O') .AND. &
     345                      (delPhase(t1%gen0Name) /= 'CLDFRA')) .OR. t1%component /= 'lmdz'
    338346      ttr(iq)       = t1
    339347
    … …  
    388396
    389397   !--- Note: nqtottr can differ from nbtr when nmom/=0
    390    nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
    391    IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
    392       CALL abort_gcm(modname, 'problem with the computation of nqtottr', 1)
     398!   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
     399   nqtottr = nqtot - COUNT(tracers(:)%component == 'lmdz' .AND. &
     400       ((delPhase(tracers(:)%gen0Name)     == 'H2O') .OR. &
     401        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
     402!   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
     403!   IF(COUNT(tracers%iso_iName == 0) - COUNT(tracers(:)%component == 'lmdz' .AND. &
     404!       ((delPhase(tracers(:)%name)     == 'H2O') .OR. &
     405!        (delPhase(tracers(:)%name)     == 'CLDFRA') /= nqtottr) &
     406!      CALL abort_gcm(modname, 'problem with the computation of nqtottr', 1)
    393407
    394408   !=== DISPLAY THE RESULTS

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

    r5601 r5609  
    4444          solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
    4545          sollw,sollwdown, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
    46     sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
     46    sig1, ftsol, cwcon, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
    4747    zmax0,fevap, rnebcon,falb_dir, falb_dif, wake_fip,    agesno,  detr_therm, pbl_tke,  &
    4848    phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
    … …  
    5050    ale_bl, ale_bl_trig, alp_bl, &
    5151    ale_wake, ale_bl_stat, AWAKE_S, &
    52     cf_ancien, qvc_ancien, cfa_ancien, qta_ancien
     52    cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien
    5353
    5454  USE comconst_mod, ONLY: pi, dtvr
    … …  
    243243  cf_ancien = 0.
    244244  qvc_ancien = 0.
     245  cwcon = 0.
    245246  cfa_ancien = 0.
    246   qta_ancien = 0.
     247  pcf_ancien = 0.
     248  qva_ancien = 0.
     249  qia_ancien = 0.
    247250 
    248251  z0m(:,:)=0 ! ym missing 5th subsurface initialization

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

    r5601 r5609  
    6363    LOGICAL               :: isAdvected  = .FALSE.              !--- "true" tracers: iadv > 0.   COUNT(isAdvected )=nqtrue
    6464    LOGICAL               :: isInPhysics = .TRUE.               !--- "true" tracers: in tr_seri. COUNT(isInPhysics)=nqtottr
    65     INTEGER               :: iso_iGroup  = 0                    !--- Isotopes group index in isotopes(:)
     65    INTEGER               :: iso_iGroup  = -1                   !--- Isotopes group index in isotopes(:)
    6666    INTEGER               :: iso_iName   = 0                    !--- Isotope  name  index in isotopes(iso_iGroup)%trac(:)
    6767    INTEGER               :: iso_iZone   = 0                    !--- Isotope  zone  index in isotopes(iso_iGroup)%zone(:)
    … …  
    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   = 'vlibfcaq'    !--- Old phases for water (no separator)
    125   CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsbfcaq'    !--- Known phases initials
     124  CHARACTER(LEN=maxlen), PARAMETER :: old_phases   = 'vlib'     !--- Old phases for water (no separator)
     125  CHARACTER(LEN=maxlen), PARAMETER :: known_phases = 'glsb'     !--- 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', 'cldVapor ', 'aviContra', 'qtContra ']
     128                                = ['gaseous  ', 'liquid   ', 'solid    ','blownSnow']
    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

    r5601 r5609  
    260260    cf_ancien = 0.
    261261    qvc_ancien = 0.
     262    cwcon = 0.
    262263
    263264    wake_delta_pbl_TKE(:,:,:)=0

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

    r5601 r5609  
    1010      USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
    1111   USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
    12        clwcon, detr_therm, &
     12       cwcon, clwcon, detr_therm, &
    1313       qsol, fevap, z0m, z0h, agesno, &
    1414       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
    … …  
    2323       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
    2424       ql_ancien, qs_ancien, qbs_ancien, &
    25        cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, &
     25       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
    2626       prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
    2727       u10m,v10m,ale_wake,ale_bl_stat
    … …  
    874874          cf_ancien = 0.
    875875          qvc_ancien = 0.
     876          cwcon = 0.
    876877        ENDIF
    877878        IF ( ok_plane_contrail ) THEN
    878879          cfa_ancien = 0.
    879           qta_ancien = 0.
     880          pcf_ancien = 0.
     881          qva_ancien = 0.
     882          qia_ancien = 0.
    880883        ENDIF
    881884!jyg<

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

    r5601 r5609  
    66   USE ioipsl, only: ju2ymds, ymds2ju, ioconf_calendar,getin
    77   USE phys_state_var_mod, ONLY : phys_state_var_init, phys_state_var_end, &
    8        clwcon, detr_therm, &
     8       cwcon, clwcon, detr_therm, &
    99       qsol, fevap, z0m, z0h, agesno, &
    1010       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
    … …  
    1919       zstd, zthe, zval, ale_bl, ale_bl_trig, alp_bl, &
    2020       ql_ancien, qs_ancien, qbs_ancien, &
    21        cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, &
     21       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
    2222       prlw_ancien, prsw_ancien, prbsw_ancien, prw_ancien, &
    2323       u10m,v10m,ale_wake,ale_bl_stat, ratqs_inter_
    … …  
    616616          cf_ancien = 0.
    617617          qvc_ancien = 0.
     618          cwcon = 0.
    618619        ENDIF
    619620        IF ( ok_plane_contrail ) THEN
    620621          cfa_ancien = 0.
    621           qta_ancien = 0.
     622          pcf_ancien = 0.
     623          qva_ancien = 0.
     624          qia_ancien = 0.
    622625        ENDIF
    623626        rain_fall=0.
    … …  
    820823      ENDIF
    821824     
     825      ! calculation of potential temperature for the advection
     826      teta=temp*(pzero/play)**rkappa
     827
    822828      ! vertical tendencies computed as d X / d t = -W  d X / d z
    823829
    … …  
    829835        d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l-1))/(z_adv(l+1)-z_adv(l-1))
    830836        ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
    831         d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l-1))/(z_adv(l+1)-z_adv(l-1))
     837        d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l-1))/(z_adv(l+1)-z_adv(l-1)) / (pzero/play(l))**rkappa
    832838        d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l-1,:))/(z_adv(l+1)-z_adv(l-1))
    833839      enddo
    … …  
    842848          d_u_vert_adv(l)=-w_adv(l)*(u(l)-u(l-1))/(z_adv(l)-z_adv(l-1))
    843849          d_v_vert_adv(l)=-w_adv(l)*(v(l)-v(l-1))/(z_adv(l)-z_adv(l-1))
    844           d_t_vert_adv(l)=-w_adv(l)*(temp(l)-temp(l-1))/(z_adv(l)-z_adv(l-1))
     850          ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
     851          d_t_vert_adv(l)=-w_adv(l)*(teta(l)-teta(l-1))/(z_adv(l)-z_adv(l-1)) / (pzero/play(l))**rkappa
    845852          d_q_vert_adv(l,:)=-w_adv(l)*(q(l,:)-q(l-1,:))/(z_adv(l)-z_adv(l-1))
    846853        ELSE
    847854          d_u_vert_adv(l)=-w_adv(l)*(u(l+1)-u(l))/(z_adv(l+1)-z_adv(l))
    848855          d_v_vert_adv(l)=-w_adv(l)*(v(l+1)-v(l))/(z_adv(l+1)-z_adv(l))
    849           d_t_vert_adv(l)=-w_adv(l)*(temp(l+1)-temp(l))/(z_adv(l+1)-z_adv(l))
     856          ! d theta / dt = -W d theta / d z, transformed into d temp / d t dividing by (pzero/play(l))**rkappa
     857          d_t_vert_adv(l)=-w_adv(l)*(teta(l+1)-teta(l))/(z_adv(l+1)-z_adv(l)) / (pzero/play(l))**rkappa
    850858          d_q_vert_adv(l,:)=-w_adv(l)*(q(l+1,:)-q(l,:))/(z_adv(l+1)-z_adv(l))
    851859        ENDIF

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

    r5536 r5609  
    293293   !---------------------------------------------------------------------------------------------------------------------------
    294294   nqtrue = SIZE(tracers)                                                                               !--- "true" tracers
    295    nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%name)     == 'H2O')     !--- Water phases
    296    nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. delPhase(tracers(:)%gen0Name) == 'H2O')     !--- Passed to phytrac
     295   nqo    =      COUNT(tracers(:)%component == 'lmdz' .AND. &
     296       ((delPhase(tracers(:)%name)     == 'H2O') .OR. &    !--- Water phases
     297        (delPhase(tracers(:)%name)     == 'CLDFRA')))
     298   nbtr = nqtrue-COUNT(tracers(:)%component == 'lmdz' .AND. &
     299       ((delPhase(tracers(:)%gen0Name)     == 'H2O') .OR. &    !--- Passed to phytrac
     300        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
    297301   nqCO2  =      COUNT( [type_trac == 'inco', type_trac == 'co2i'] )
    298302IF (CPPKEY_INCA) THEN
    … …  
    350354      t1%longName   = t1%name; IF(iad > 0) t1%longName=TRIM(t1%name)//descrq(iad)
    351355      t1%isAdvected = iad >= 0
    352       t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' .OR. t1%component /= 'lmdz' !=== MORE EXCEPTIONS ? CO2i, SURSAT CLOUD H2O
     356      !t1%isInPhysics= delPhase(t1%gen0Name) /= 'H2O' .OR. t1%component /= 'lmdz' !=== MORE EXCEPTIONS ? CO2i, SURSAT CLOUD H2O
     357      t1%isInPhysics=((delPhase(t1%gen0Name) /= 'H2O') .AND. &
     358                      (delPhase(t1%gen0Name) /= 'CLDFRA')) .OR. t1%component /= 'lmdz'
    353359      ttr(iq)       = t1
    354360
    … …  
    397403
    398404   !--- Note: nqtottr can differ from nbtr when nmom/=0
    399    nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
    400    IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
    401       CALL abort_physic(modname, 'problem with the computation of nqtottr', 1)
     405!   nqtottr = nqtot - COUNT(delPhase(tracers%gen0Name) == 'H2O' .AND. tracers%component == 'lmdz')
     406   nqtottr = nqtot - COUNT(tracers(:)%component == 'lmdz' .AND. &
     407       ((delPhase(tracers(:)%gen0Name)     == 'H2O') .OR. &
     408        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
     409!   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
     410!   IF(COUNT(tracers%iso_iName == 0) - COUNT(tracers(:)%component == 'lmdz' .AND. &
     411!       ((delPhase(tracers(:)%name)     == 'H2O') .OR. &
     412!        (delPhase(tracers(:)%name)     == 'CLDFRA'))) /= nqtottr) &
     413!      CALL abort_physic(modname, 'problem with the computation of nqtottr', 1)
    402414
    403415   !=== DISPLAY THE RESULTS

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

    r5602 r5609  
    7373SUBROUTINE contrails_formation( &
    7474      klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, flight_dist, &
    75       clrfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &
     75      clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, keepgoing, pt_pron_clds, &
    7676      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    7777      dcfa_ini, dqia_ini, dqta_ini)
    … …  
    9797REAL,     INTENT(IN) , DIMENSION(klon) :: flight_dist  ! aviation distance flown concentration [m/s/m3]
    9898REAL,     INTENT(IN) , DIMENSION(klon) :: clrfra       ! clear sky fraction [-]
    99 REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_loc      ! location parameter of the clear sky PDF [%]
     99REAL,     INTENT(IN) , DIMENSION(klon) :: qclr         ! clear sky specific humidity [kg/kg]
    100100REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_scale    ! scale parameter of the clear sky PDF [%]
    101101REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_alpha    ! shape parameter of the clear sky PDF [-]
     102REAL,     INTENT(IN) , DIMENSION(klon) :: pdf_gamma    ! tmp parameter of the clear sky PDF [-]
    102103LOGICAL,  INTENT(IN) , DIMENSION(klon) :: keepgoing    ! .TRUE. if a new condensation loop should be computed
     104LOGICAL,  INTENT(IN) , DIMENSION(klon) :: pt_pron_clds  ! .TRUE. if clouds are prognostic in this mesh
    103105!
    104106! Output
    … …  
    119121REAL, DIMENSION(klon) :: qsatl_crit, dqsatl_crit
    120122REAL :: Gcont, psatl_crit, pcrit
    121 REAL :: pdf_x, pdf_y, pdf_e2, pdf_e3
    122 REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat, pdf_fra_above_qnuc
    123 REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat, pdf_q_above_qnuc
     123REAL :: rhl_clr, pdf_loc
     124REAL :: pdf_x, pdf_y, pdf_e3
     125REAL :: pdf_fra_above_qcritcont, pdf_fra_above_qsat
     126REAL :: pdf_q_above_qcritcont, pdf_q_above_qsat
    124127REAL :: qpotcontP
    125128!
    … …  
    159162
    160163DO i = 1, klon
    161   IF ( keepgoing(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
     164  IF ( keepgoing(i) .AND. pt_pron_clds(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
    162165   
    163166    psatl_crit = qsatl_crit(i) / ( EPS_W + ( 1. - EPS_W ) * qsatl_crit(i) ) * pplay(i)
    … …  
    168171    IF ( ( clrfra(i) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN
    169172   
     173      rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     174      pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
    170175      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))
     176      pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    173177      pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    174       pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     178      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
    175179      pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i)
    176180      pdf_q_above_qcritcont = ( pdf_e3 * clrfra(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))
     181          + pdf_loc * pdf_fra_above_qcritcont ) * qsatl(i) / 100.
     182   
     183      pdf_x = qsat(i) / qsatl(i) * 100.
     184      pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    182185      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 ) * clrfra(i)
    185       pdf_q_above_qnuc = ( pdf_e3 * clrfra(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
     186      pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
    193187      pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i)
    194188      pdf_q_above_qsat = ( pdf_e3 * clrfra(i) &
    195           + pdf_loc(i) * pdf_fra_above_qsat ) * qsatl(i) / 100.
     189          + pdf_loc * pdf_fra_above_qsat ) * qsatl(i) / 100.
    196190   
    197191      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))
     192      potcontfraP(i) = MIN(pdf_fra_above_qsat, pdf_fra_above_qcritcont)
     193      qpotcontP = MIN(pdf_q_above_qsat, pdf_q_above_qcritcont)
    202194
    203195    ENDIF ! temp .LT. Tcritcont
    … …  
    219211
    220212!**********************************************************************************
    221 FUNCTION QVAPINCLD_DEPSUB_CONTRAILS( &
    222     qvapincld, qiceincld, temp, qsat, pplay, dtime)
    223 
    224 USE lmdz_lscp_ini,        ONLY: RV, RLSTT, RTT, EPS_W
    225 USE lmdz_lscp_ini,        ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
    226 USE lmdz_lscp_ini,        ONLY: rm_ice_crystals_contrails
     213FUNCTION contrail_cross_section_onera()
     214
     215USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
    227216
    228217IMPLICIT NONE
    229218
    230 ! inputs
    231 REAL :: qvapincld     !
    232 REAL :: qiceincld     !
    233 REAL :: temp          !
    234 REAL :: qsat          !
    235 REAL :: pplay         !
    236 REAL :: dtime         ! time step [s]
    237 ! output
    238 REAL :: qvapincld_depsub_contrails
    239 ! local
    240 REAL :: pres_sat, rho, kappa
    241 REAL :: air_thermal_conduct, water_vapor_diff
    242 REAL :: rm_ice
    243 
    244 !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
    245 !--hypothesis is lost, and the vapor in the cloud is purely prognostic.
    246 !
    247 !--The deposition equation is
    248 !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
    249 !--from Lohmann et al. (2016), where
    250 !--alpha is the deposition coefficient [-]
    251 !--mi is the mass of one ice crystal [kg]
    252 !--C is the capacitance of an ice crystal [m]
    253 !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
    254 !--R_v is the specific gas constant for humid air [J/kg/K]
    255 !--T is the temperature [K]
    256 !--esi is the saturation pressure w.r.t. ice [Pa]
    257 !--Dv is the diffusivity of water vapor [m2/s]
    258 !--Ls is the specific latent heat of sublimation [J/kg/K]
    259 !--ka is the thermal conductivity of dry air [J/m/s/K]
    260 !
    261 !--alpha is a coefficient to take into account the fact that during deposition, a water
    262 !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
    263 !--coherent (with the same structure). It has no impact for sublimation.
    264 !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
    265 !--during deposition, and alpha = 1. during sublimation.
    266 !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
    267 !-- C = capa_cond_cirrus * rm_ice
    268 !
    269 !--We have qice = Nice * mi, where Nice is the ice crystal
    270 !--number concentration per kg of moist air
    271 !--HYPOTHESIS 1: the ice crystals are spherical, therefore
    272 !-- mi = 4/3 * pi * rm_ice**3 * rho_ice
    273 !--HYPOTHESIS 2: the ice crystals are monodisperse with the
    274 !--initial radius rm_ice_0.
    275 !--NB. this is notably different than the assumption
    276 !--of a distributed qice in the cloud made in the sublimation process;
    277 !--should it be consistent?
    278 !
    279 !--As the deposition process does not create new ice crystals,
    280 !--and because we assume a same rm_ice value for all crystals
    281 !--therefore the sublimation process does not destroy ice crystals
    282 !--(or, in a limit case, it destroys all ice crystals), then
    283 !--Nice is a constant during the sublimation/deposition process.
    284 !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
    285 !
    286 !--The deposition equation then reads:
    287 !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
    288 !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
    289 !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
    290 !--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
    291 !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
    292 !--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
    293 !--and we have
    294 !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
    295 !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
    296 !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
    297 !
    298 !--This system of equations can be resolved with an exact
    299 !--explicit numerical integration, having one variable resolved
    300 !--explicitly, the other exactly. The exactly resolved variable is
    301 !--the one decreasing, so it is qvc if the process is
    302 !--condensation, qi if it is sublimation.
    303 !
    304 !--kappa is computed as an initialisation constant, as it depends only
    305 !--on temperature and other pre-computed values
    306 pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
    307 !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
    308 air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
    309 !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
    310 water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
    311 kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
    312       * ( RV * temp / water_vapor_diff / pres_sat &
    313         + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
    314 !--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
    315 
    316 !--Here, the initial vapor in the cloud is qvapincld, and we compute
    317 !--the new vapor qvapincld_depsub_contrails
    318 
    319 rm_ice = rm_ice_crystals_contrails
    320 
    321 IF ( qvapincld .GE. qsat ) THEN
    322   !--If the cloud is initially supersaturated
    323   !--Exact explicit integration (qvc exact, qice explicit)
    324   qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) &
    325                 * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
    326 ELSE
    327   !--If the cloud is initially subsaturated
    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 )
    332   IF ( qvapincld_depsub_contrails .GE. ( qvapincld + qiceincld ) ) &
    333                 qvapincld_depsub_contrails = 0.
    334 ENDIF ! qvapincld .GT. qsat
    335 
    336 END FUNCTION QVAPINCLD_DEPSUB_CONTRAILS
    337 
    338 
    339 !**********************************************************************************
    340 FUNCTION contrail_cross_section_onera()
    341 
    342 USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
    343 
    344 IMPLICIT NONE
    345 
    346219!
    347220! Output
    … …  
    355228
    356229END FUNCTION contrail_cross_section_onera
     230
    357231
    358232SUBROUTINE read_aviation_emissions(klon, klev)

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

    r5601 r5609  
    1212    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
    1313    !--AB contrails
    14     contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
    15     pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     14    contfra, perscontfra, qice_cont, qice_perscont, pclc_nocont, &
     15    pcltau_nocont, pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
    1616    xfiwp_nocont, xfiwc_nocont, reice_nocont)
    1717
    … …  
    9595
    9696  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
    97   REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
    98   REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
     97  REAL, INTENT(IN)  :: contfra(klon, klev)       ! linear contrails fraction [-]
     98  REAL, INTENT(IN)  :: perscontfra(klon, klev)   ! contrail induced cirrus fraction [-]
     99  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! linear contrails condensed water [kg/kg]
     100  REAL, INTENT(IN)  :: qice_perscont(klon, klev) ! contrail induced cirrus condensed water [kg/kg]
    99101  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
    100102  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
    … …  
    137139    icefrac_optics, dNovrN, ptconv,rnebcon, ccwcon, &
    138140    !--AB for contrails
    139     contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
    140     pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     141    contfra, perscontfra, qice_cont, qice_perscont, pclc_nocont, pcltau_nocont, &
     142    pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
    141143    xfiwp_nocont, xfiwc_nocont, reice_nocont)
    142144  ELSE

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

    r5601 r5609  
    1111    icefrac_optics, dNovrN, ptconv, rnebcon, ccwcon, &
    1212    !--AB contrails
    13     contfra, qice_cont, pclc_nocont, pcltau_nocont, pclemi_nocont, &
    14     pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
     13    contfra, perscontfra, qice_cont, qice_perscont, pclc_nocont, pcltau_nocont, &
     14    pclemi_nocont, pcltau_cont, pclemi_cont, pch_nocont, pct_cont, &
    1515    xfiwp_nocont, xfiwc_nocont, reice_nocont)
    1616
    … …  
    119119
    120120  !--AB for contrails. All these are used / outputed only if ok_plane_contrail=y
    121   REAL, INTENT(IN)  :: contfra(klon, klev)       ! contrails fraction [-]
    122   REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! contrails condensed water [kg/kg]
     121  REAL, INTENT(IN)  :: contfra(klon, klev)       ! linear contrails fraction [-]
     122  REAL, INTENT(IN)  :: perscontfra(klon, klev)   ! contrail induced cirrus fraction [-]
     123  REAL, INTENT(IN)  :: qice_cont(klon, klev)     ! linear contrails condensed water [kg/kg]
     124  REAL, INTENT(IN)  :: qice_perscont(klon, klev) ! contrail induced cirrus condensed water [kg/kg]
    123125  REAL, INTENT(OUT) :: pch_nocont(klon)          ! 2D high cloud cover without contrails[-]
    124126  REAL, INTENT(OUT) :: pct_cont(klon)            ! 2D total contrails cover[-]
    … …  
    171173  REAL zflwp_var, zfiwp_var
    172174  REAL d_rei_dt
    173   REAL pclc_var
     175  REAL pclc_cont, pclc_perscont
     176  REAL rei_cont, rei_perscont
    174177
    175178
    … …  
    252255      DO k = 1, klev
    253256        DO i = 1, klon
    254           pclc_nocont(i,k) = pclc(i, k) - contfra(i, k)
    255           xfiwc_nocont(i, k) = xfiwc(i, k) - qice_cont(i, k)
     257          pclc_nocont(i,k) = pclc(i, k) - contfra(i, k) - perscontfra(i, k)
     258          xfiwc_nocont(i, k) = xfiwc(i, k) - qice_cont(i, k) - qice_perscont(i, k)
    256259        ENDDO
    257260      ENDDO
    … …  
    393396
    394397          IF (zfiwp_var==0. .OR. rei<=0.) rei = 1.
    395           IF ( ok_plane_contrail ) THEN
    396             !--If contrails are activated, rei is a weighted average between the natural
    397             !--rei and the contrails rei, with the weights being the fraction of natural
    398             !--vs contrail cirrus in the gridbox
    399             !--Beware, re_ice_crystals_contrails is in m, while rei is in microns
    400             rei = ( rei * pclc_nocont(i,k) &
    401                 + re_ice_crystals_contrails * 1.E6 * contfra(i,k) ) / pclc(i,k)
    402           ENDIF
     398
    403399          pcldtaupi(i, k) = 3.0/2.0*zflwp_var/rad_chaud_pi(i, k) + &
    404400            zfiwp_var*(3.448E-03+2.431/rei)
    … …  
    687683          reice_nocont(i,k) = 1.
    688684          pclc_nocont(i,k) = 0.
    689           pclc(i,k) = contfra(i,k)
     685          pclc(i,k) = contfra(i,k) + perscontfra(i,k)
    690686          pcltau_nocont(i, k) = 0.
    691687          pclemi_nocont(i, k) = 0.
    … …  
    695691
    696692        IF ( contfra(i,k) .GT. 0.01 * seuil_neb ) THEN
     693          pclc_cont = contfra(i,k)
     694          rei_cont = re_ice_crystals_contrails * 1.E6
     695        ELSE
     696          pclc_cont = 0.
     697          rei_cont = 1.
     698        ENDIF
     699
     700
     701        IF ( perscontfra(i,k) .GT. 0.01 * seuil_neb ) THEN
    697702          !--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))
     703          pclc_perscont = perscontfra(i,k)
     704
     705          tc = temp(i, k) - 273.15
     706          rei_perscont = d_rei_dt*tc + rei_max
     707          IF (tc<=-81.4) rei_perscont = rei_min
     708
     709        ELSE
     710          pclc_perscont = 0.
     711          rei_perscont = 1.
     712        ENDIF
     713
     714        IF ( MAX(contfra(i,k), perscontfra(i,k)) .GT. 0.01 * seuil_neb ) THEN
     715
     716          rei = ( rei_cont * pclc_cont + rei_perscont * pclc_perscont ) &
     717              / ( pclc_cont + pclc_perscont )
     718          zfiwp_var = 1000.*(xfiwc(i, k)-xfiwc_nocont(i, k))&
     719                    / (pclc(i, k)-pclc_nocont(i, k))*rhodz(i, k)
     720
     721          pcltau_cont(i, k) = zfiwp_var*(3.448E-03+2.431/rei)
    702722          ! -- cloud infrared emissivity:
    703723          ! [the broadband infrared absorption coefficient is PARAMETERized
    704724          ! as a function of the effective cld droplet radius]
    705725          ! Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
    706           k_ice = k_ice0 + 1.0/(re_ice_crystals_contrails*1.E6)
     726          k_ice = k_ice0 + 1.0/rei
    707727          pclemi_cont(i, k) = 1.0 - exp(-df*k_ice*zfiwp_var)
    708728
    709729        ELSE
    710           pclc_var = 0.
    711730          pcltau_cont(i, k) = 0.
    712731          pclemi_cont(i, k) = 0.
    713732        ENDIF
    714733
    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) )
     734
     735        rei = ( rei_cont * pclc_cont + rei_perscont * pclc_perscont &
     736            + reice_nocont(i, k) * pclc_nocont(i, k) ) &
     737            / ( pclc_cont + pclc_perscont + pclc_nocont(i,k) )
    717738
    718739        zflwp_var = 1000.*xflwc(i, k)/pclc(i, k)*rhodz(i, k)
    … …  
    861882    DO k = klev, 1, -1
    862883      DO i = 1, klon
    863         zclear(i) = zclear(i)*(1.-max(contfra(i,k),zcloud(i)))/(1.-min(real( &
    864           zcloud(i),kind=8),1.-zepsec))
     884        zclear(i) = zclear(i)*(1.-max(contfra(i,k)+perscontfra(i,k),zcloud(i)))/&
     885          (1.-min(real(zcloud(i),kind=8),1.-zepsec))
    865886        pct_cont(i) = 1. - zclear(i)
    866         zcloud(i) = contfra(i,k)
     887        zcloud(i) = contfra(i,k) + perscontfra(i,k)
    867888        IF (paprs(i,k)<prmhc) THEN
    868889          pch_nocont(i) = pch_nocont(i)*(1.-max(pclc_nocont(i,k),zcloudh(i)))/(1.-min(real( &

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

    r5602 r5609  
    88SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
    99     paprs, pplay, omega, temp, qt, ql_seri, qi_seri,   &
    10      ptconv, cldfracv, qcondcv, ratqs, sigma_qtherm,    &
     10     ratqs, sigma_qtherm, ptconv, cfcon_old, qvcon_old, &
     11     qccon_old, cfcon, qvcon, qccon,                    &
    1112     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
    1213     pfraclr, pfracld,                                  &
    … …  
    2425     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
    2526     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
    26      cfa_seri, qta_seri, flight_dist, flight_h2o,       &
    27      qice_cont, Tcritcont,                              &
     27     cfa_seri, pcf_seri, qva_seri, qia_seri,flight_dist,&
     28     flight_h2o, qice_perscont, Tcritcont,              &
    2829     qcritcont, potcontfraP, potcontfraNP,              &
    2930     cloudth_sth,                                       &
    … …  
    117118USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
    118119USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
    119 USE lmdz_lscp_ini, ONLY : t_glace_min, min_frac_th_cld
     120USE lmdz_lscp_ini, ONLY : t_glace_min, temp_nowater, min_frac_th_cld
    120121USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
    121122USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
    122123USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
    123 USE lmdz_lscp_ini, ONLY : ok_plane_contrail
     124USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato
     125USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_ice_sedim
    124126
    125127! Temporary call for Lamquin et al (2012) diagnostics
    … …  
    153155                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
    154156  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
    155   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cldfracv        ! cloud fraction from deep convection [-]
    156   REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qcondcv         ! in-cloud condensed specific humidity from deep convection [kg/kg]
     157  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon_old       ! cloud fraction from deep convection from previous timestep [-]
     158  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon_old       ! in-cloud vapor specific humidity from deep convection from previous timestep [kg/kg]
     159  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon_old       ! in-cloud condensed specific humidity from deep convection from previous timestep [kg/kg]
     160  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon           ! cloud fraction from deep convection [-]
     161  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon           ! in-cloud vapor specific humidity from deep convection [kg/kg]
     162  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon           ! in-cloud condensed specific humidity from deep convection [kg/kg]
    157163
    158164  !Inputs associated with thermal plumes
    … …  
    185191  !--------------------------------------------------
    186192  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfa_seri         ! linear contrails fraction [-]
    187   REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qta_seri         ! linear contrails total specific humidity [kg/kg]
     193  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: pcf_seri         ! contrails induced cirrus fraction [-]
     194  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qva_seri         ! linear contrails total specific humidity [kg/kg]
     195  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qia_seri         ! linear contrails total specific humidity [kg/kg]
    188196  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown within the mesh [m/s/mesh]
    189197  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! aviation H2O emitted within the mesh [kgH2O/s/mesh]
    … …  
    243251  ! for contrails and aviation
    244252
    245   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_cont      !--condensed water in linear contrails used in the radiation scheme [kg/kg]
     253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_perscont  !--condensed water in contrail induced cirrus used in the radiation scheme [kg/kg]
    246254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
    247255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
    … …  
    279287  ! LOCAL VARIABLES:
    280288  !----------------
    281   REAL, DIMENSION(klon) :: qliq_in, qice_in
     289  REAL, DIMENSION(klon) :: qliq_in, qice_in, qvc_in, cldfra_in
    282290  REAL, DIMENSION(klon,klev) :: ctot
    283291  REAL, DIMENSION(klon,klev) :: ctot_vol
    … …  
    307315  REAL, DIMENSION(klon) :: ztupnew
    308316  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
    309   REAL, DIMENSION(klon) :: zqice_sedim ! for sedimentation of ice crystals
     317  REAL, DIMENSION(klon) :: cldfra_above, icesed_flux ! for sedimentation of ice crystals
    310318  REAL, DIMENSION(klon) :: zrflclr, zrflcld
    311319  REAL, DIMENSION(klon) :: ziflclr, ziflcld
    … …  
    324332  REAL :: delta_z
    325333  ! for contrails
    326   REAL, DIMENSION(klon) :: contfra, qcont
     334  REAL, DIMENSION(klon) :: contfra, perscontfra, qcont
     335  LOGICAL, DIMENSION(klon) :: pt_pron_clds
    327336  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrail)
    328337  ! Constants used for calculating ratios that are advected (using a parent-child
    … …  
    424433dcf_con(:,:)    = 0.
    425434dcf_mix(:,:)    = 0.
     435dcfsed(:,:)     = 0.
    426436dqi_adj(:,:)    = 0.
    427437dqi_sub(:,:)    = 0.
    428438dqi_con(:,:)    = 0.
    429439dqi_mix(:,:)    = 0.
     440dqised(:,:)     = 0.
    430441dqvc_adj(:,:)   = 0.
    431442dqvc_sub(:,:)   = 0.
    432443dqvc_con(:,:)   = 0.
    433444dqvc_mix(:,:)   = 0.
     445dqvcsed(:,:)    = 0.
    434446qvc(:)          = 0.
    435447shear(:)        = 0.
    … …  
    467479zqupnew(:)    = 0.
    468480zqvapclr(:)   = 0.
    469 zqice_sedim(:)= 0.
     481cldfra_above(:)= 0.
     482icesed_flux(:)= 0.
    470483
    471484
    … …  
    506519        !c_iso init of iso
    507520    ENDDO
     521    IF ( ok_ice_supersat ) THEN
     522      cldfra_in(:) = cf_seri(:,k)
     523      qvc_in(:) = qvc_seri(:,k)
     524    ENDIF
    508525
    509526    ! --------------------------------------------------------------------
    … …  
    519536      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    520537                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
    521                         zqvapclr, zqupnew, zqice_sedim, &
    522                         cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, &
     538                        zqvapclr, zqupnew, icesed_flux, &
     539                        cldfra_in, qvc_in, qliq_in, qice_in, &
    523540                        zrfl, zrflclr, zrflcld, &
    524541                        zifl, ziflclr, ziflcld, &
    525                         dqreva(:,k), dqssub(:,k), &
    526                         dcfsed(:,k), dqvcsed(:,k) &
     542                        dqreva(:,k), dqssub(:,k) &
    527543                        )
    528544
    … …  
    531547
    532548      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    533                         zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
     549                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
    534550                        zrfl, zrflclr, zrflcld, &
    535551                        zifl, ziflclr, ziflcld, &
    … …  
    646662            lognormale(:) = .FALSE.
    647663       
     664        ENDIF
     665
     666
     667        IF ( ok_ice_supersat ) THEN
     668
     669          !--Initialisation
     670          IF ( ok_plane_contrail ) THEN
     671            contfra(:)     = 0.
     672            qcont(:)       = 0.
     673            perscontfra(:) = 0.
     674          ENDIF
     675
     676          DO i = 1, klon
     677
     678            pt_pron_clds(i) = ( ( ( ( zt(i) .LE. temp_nowater ) .OR. ok_weibull_warm_clouds ) &
     679                .AND. ( .NOT. ok_no_issr_strato .OR. ( stratomask(i,k) .EQ. 0. ) ) ) &
     680                .AND. ( cfcon(i,k) .LT. ( 1. - eps ) ) )
     681
     682            IF ( pt_pron_clds(i) ) THEN
     683
     684              !--If deep convection is activated, the condensation scheme activates
     685              !--only in the environment. NB. the clear sky fraction will the be
     686              !--maximised by 1. - cfcon(i,k)
     687              IF ( ptconv(i,k) ) zq(i) = zq(i) - ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
     688
     689              IF ( ( cfcon(i,k) * qccon(i,k) ) .LT. ( cfcon_old(i,k) * qccon_old(i,k) ) ) THEN
     690                !--If deep convection is weakening, we add the clouds that are not anymore
     691                !--'in' deep convection to the advected clouds
     692                cldfra_in(i) = cldfra_in(i) + MAX(0., cfcon_old(i,k) - cfcon(i,k))
     693                qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * MAX(0., cfcon_old(i,k) - cfcon(i,k))
     694                qice_in(i) = qice_in(i) + ( qccon_old(i,k) * cfcon_old(i,k) &
     695                                          - qccon(i,k) * cfcon(i,k) )
     696              ELSEIF ( cfcon(i,k) .GT. cfcon_old(i,k) ) THEN
     697                !--Else if deep convection is strengthening, it consumes the existing cloud
     698                !--fraction (which does not at this moment represent deep convection)
     699                !--NB. if deep convection is strengthening while the fraction decreases,
     700                !--clear sky water vapor will be transfered in priority
     701                cldfra_in(i) = cldfra_in(i) * ( 1. &
     702                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
     703                qvc_in(i)    = qvc_in(i)    * ( 1. &
     704                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
     705                qice_in(i)   = qice_in(i)   * ( 1. &
     706                             - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
     707              ENDIF
     708
     709              !--Barriers
     710              cldfra_in(i) = MAX(0., MIN(1. - cfcon(i,k), cldfra_in(i)))
     711              qvc_in(i)    = MAX(0., MIN(zq(i), qvc_in(i)))
     712              qice_in(i)   = MAX(0., MIN(zq(i) - qvc_in(i), qice_in(i)))
     713
     714              !--Calculate the shear value (input for condensation and ice supersat)
     715              !--Cell thickness [m]
     716              delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(i) * RD
     717              IF ( iftop ) THEN
     718                ! top
     719                shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
     720                               + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
     721              ELSEIF ( k .EQ. 1 ) THEN
     722                ! surface
     723                shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
     724                               + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
     725              ELSE
     726                ! other layers
     727                shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
     728                                   - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
     729                               + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
     730                                   - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
     731              ENDIF
     732            ENDIF
     733          ENDDO
    648734        ENDIF
    649735
    … …  
    717803                  IF (ok_ice_supersat) THEN
    718804
    719                     !--Calculate the shear value (input for condensation and ice supersat)
    720                     DO i = 1, klon
    721                       !--Cell thickness [m]
    722                       delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
    723                       IF ( iftop ) THEN
    724                         ! top
    725                         shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
    726                                        + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
    727                       ELSEIF ( k .EQ. 1 ) THEN
    728                         ! surface
    729                         shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
    730                                        + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
    731                       ELSE
    732                         ! other layers
    733                         shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
    734                                            - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
    735                                        + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
    736                                            - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
    737                       ENDIF
    738                     ENDDO
     805                    IF ( iftop ) THEN
     806                      cldfra_above(:) = 0.
     807                    ELSE
     808                      cldfra_above(:) = rneb(:,k+1)
     809                    ENDIF
    739810
    740811                    !---------------------------------------------
    … …  
    744815                    CALL condensation_ice_supersat( &
    745816                        klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), &
    746                         cldfracv(:,k), qcondcv(:,k), &
    747                         cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, &
    748                         shear, tke_dissip(:,k), cell_area, stratomask(:,k), &
    749                         Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
     817                        cfcon(:,k), cldfra_in, qvc_in, qliq_in, qice_in, &
     818                        shear, tke_dissip(:,k), cell_area, Tbef, zq, zqs, &
     819                        gammasat, ratqs(:,k), keepgoing, pt_pron_clds, &
     820                        cldfra_above, icesed_flux,&
    750821                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
    751                         dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
    752                         dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
    753                         dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
    754                         cfa_seri(:,k), qta_seri(:,k), flight_dist(:,k), flight_h2o(:,k), &
    755                         contfra, qcont, Tcritcont(:,k), qcritcont(:,k), &
    756                         potcontfraP(:,k), potcontfraNP(:,k), &
     822                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), dcfsed(:,k), &
     823                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), dqised(:,k), &
     824                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), dqvcsed(:,k), &
     825                        cfa_seri(:,k), pcf_seri(:,k), qva_seri(:,k), qia_seri(:,k), &
     826                        flight_dist(:,k), flight_h2o(:,k), contfra, perscontfra, qcont, &
     827                        Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
    757828                        dcfa_ini(:,k), dqia_ini(:,k), dqta_ini(:,k), &
    758829                        dcfa_sub(:,k), dqia_sub(:,k), dqta_sub(:,k), &
    … …  
    9511022
    9521023    IF (ok_plane_contrail) THEN
    953       !!--Contrails do not precipitate. We remove then from the variables temporarily
    954       !DO i = 1, klon
    955       !  rneb(i,k) = rneb(i,k) - contfra(i)
    956       !  zoliqi(i) = zoliqi(i) - ( qcont(i) - zqs(i) * contfra(i) )
    957       !ENDDO
    9581024      !--Contrails precipitate as natural clouds. We save the partition of ice
    9591025      !--between natural clouds and contrails
    … …  
    9751041                            ctot_vol(:,k), ptconv(:,k), &
    9761042                            zt, zq, zoliql, zoliqi, zfice, &
    977                             rneb(:,k), zqice_sedim, znebprecipclr, znebprecipcld, &
     1043                            rneb(:,k), icesed_flux, znebprecipclr, znebprecipcld, &
    9781044                            zrfl, zrflclr, zrflcld, &
    9791045                            zifl, ziflclr, ziflcld, &
    … …  
    9911057
    9921058      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    993                             ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
    994                             zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
     1059                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
     1060                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
    9951061                            rneb(:,k), znebprecipclr, znebprecipcld, &
    9961062                            zneb, tot_zneb, zrho_up, zvelo_up, &
    9971063                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
    998                             zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
     1064                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k), dqised(:,k) &
    9991065                            )
    10001066
    … …  
    10021068   
    10031069    IF (ok_plane_contrail) THEN
    1004       !!--Contrails are reintroduced in the variables
    1005       !DO i = 1, klon
    1006       !  rneb(i,k) = rneb(i,k) + contfra(i)
    1007       !  zoliqi(i) = zoliqi(i) + ( qcont(i) - zqs(i) * contfra(i) )
    1008       !ENDDO
    10091070      !--Contrails fraction is left unchanged, but contrails water has changed
    10101071      DO i = 1, klon
    … …  
    11061167    ! P6 > write diagnostics and outputs
    11071168    !------------------------------------------------------------
     1169
     1170    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
     1171    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
    11081172   
    11091173    !--AB Write diagnostics and tracers for ice supersaturation
    11101174    IF ( ok_plane_contrail ) THEN
    11111175      cfa_seri(:,k) = contfra(:)
    1112       qta_seri(:,k) = qcont(:)
    1113       qice_cont(:,k) = qcont(:) - zqs(:) * contfra(:)
     1176      pcf_seri(:,k) = perscontfra(:)
     1177      qva_seri(:,k) = zqs(:) * contfra(:)
     1178      qia_seri(:,k) = qcont(:) - zqs(:) * contfra(:)
     1179      DO i = 1, klon
     1180        IF ( ( rneb(i,k) - cfa_seri(i,k) ) .GT. eps ) THEN
     1181          qice_perscont(i,k) = ( radocond(i,k) - qia_seri(i,k) ) &
     1182              * perscontfra(i) / ( rneb(i,k) - cfa_seri(i,k) )
     1183        ELSE
     1184          qice_perscont(i,k) = 0.
     1185        ENDIF
     1186      ENDDO
    11141187    ENDIF
    11151188
    1116     IF ( ok_ice_supersat .AND. .NOT. ok_unadjusted_clouds) qvc(:) = zqs(:) * rneb(:,k)
    1117 
    11181189    IF ( ok_ice_supersat ) THEN
    1119       CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
    1120       CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
    11211190
    11221191      DO i = 1, klon
     1192
     1193        !--If prognostic clouds are activated, deep convection vapor is
     1194        !--re-added to the total water vapor
     1195        IF ( ptconv(i,k) .AND. pt_pron_clds(i) ) &
     1196            zq(i) = zq(i) + ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
    11231197
    11241198        cf_seri(i,k) = rneb(i,k)
    … …  
    12361310  ENDDO
    12371311
     1312  IF ( ok_ice_sedim ) THEN
     1313    DO i = 1, klon
     1314      snow(i) = snow(i) + icesed_flux(i)
     1315    ENDDO
     1316  ENDIF
     1317
    12381318  IF (ncoreczq>0) THEN
    12391319      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'

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

    r5602 r5609  
    9494!**********************************************************************************
    9595SUBROUTINE condensation_ice_supersat( &
    96       klon, dtime, pplay, paprsdn, paprsup, cldfracv, qcondcv, &
    97       cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, stratomask, &
    98       temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, &
    99       cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
    100       dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, &
    101       contfra_in, qcont_in, flight_dist, flight_h2o, contfra, qcont, &
    102       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
     96      klon, dtime, pplay, paprsdn, paprsup, cfcon, &
     97      cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, &
     98      temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, pt_pron_clds, &
     99      cldfra_above, icesed_flux, &
     100      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, dcf_sed, &
     101      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, &
     102      dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, dqvc_sed, &
     103      contfra_in, perscontfra_in, qva_in, qia_in, flight_dist, flight_h2o, &
     104      contfra, perscontfra, qcont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    103105      dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub, &
    104106      dcfa_cir, dqta_cir, dcfa_mix, dqia_mix, dqta_mix)
    … …  
    118120
    119121USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf, calc_gammasat, GAMMAINC
    120 USE lmdz_lscp_ini,   ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI, EPS_W
    121 USE lmdz_lscp_ini,   ONLY: eps, temp_nowater, ok_weibull_warm_clouds
    122 USE lmdz_lscp_ini,   ONLY: ok_unadjusted_clouds, ok_no_issr_strato
    123 
    124 USE lmdz_lscp_ini,   ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp
    125 USE lmdz_lscp_ini,   ONLY: beta_pdf_lscp, temp_thresh_pdf_lscp
     122USE lmdz_lscp_ini,   ONLY: RLSTT, RTT, RD, RG, RV, RPI, EPS_W
     123USE lmdz_lscp_ini,   ONLY: eps, temp_nowater, ok_unadjusted_clouds
     124
     125USE lmdz_lscp_ini,   ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
     126USE lmdz_lscp_ini,   ONLY: mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp
    126127USE lmdz_lscp_ini,   ONLY: std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp
    127128USE lmdz_lscp_ini,   ONLY: coef_mixing_lscp, coef_shear_lscp
    … …  
    143144REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsdn       ! pressure at the lower interface [Pa]
    144145REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsup       ! pressure at the upper interface [Pa]
    145 REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfracv      ! cloud fraction from deep convection [-]
    146 REAL,     INTENT(IN)   , DIMENSION(klon) :: qcondcv       ! in-cloud condensed specific humidity from deep convection [kg/kg]
     146REAL,     INTENT(IN)   , DIMENSION(klon) :: cfcon         ! cloud fraction from deep convection [-]
    147147REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfra_in     ! cloud fraction [-]
    148148REAL,     INTENT(IN)   , DIMENSION(klon) :: qvc_in        ! gridbox-mean water vapor in cloud [kg/kg]
    … …  
    152152REAL,     INTENT(IN)   , DIMENSION(klon) :: pbl_eps       ! TKE dissipation [m2/s3]
    153153REAL,     INTENT(IN)   , DIMENSION(klon) :: cell_area     ! cell area [m2]
    154 REAL,     INTENT(IN)   , DIMENSION(klon) :: stratomask    ! fraction of stratosphere in the mesh (1 or 0)
    155154REAL,     INTENT(IN)   , DIMENSION(klon) :: temp          ! temperature [K]
    156155REAL,     INTENT(IN)   , DIMENSION(klon) :: qtot_in       ! total specific humidity (without precip) [kg/kg]
    … …  
    159158REAL,     INTENT(IN)   , DIMENSION(klon) :: ratqs         ! ratio between the variance of the total water distribution and its average [-]
    160159LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
     160LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: pt_pron_clds  ! .TRUE. if clouds are prognostic in this mesh
     161REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfra_above  ! cloud fraction IN THE LAYER ABOVE [-]
     162REAL,     INTENT(IN)   , DIMENSION(klon) :: icesed_flux   ! sedimentated ice flux [kg/s/m2]
    161163!
    162164! Input for aviation
    163165!
    164166REAL,     INTENT(IN)   , DIMENSION(klon) :: contfra_in    ! input linear contrails fraction [-]
    165 REAL,     INTENT(IN)   , DIMENSION(klon) :: qcont_in      ! input linear contrails total specific humidity [kg/kg]
     167REAL,     INTENT(IN)   , DIMENSION(klon) :: perscontfra_in! input contrail induced cirrus fraction [-]
     168REAL,     INTENT(IN)   , DIMENSION(klon) :: qva_in        ! input linear contrails total specific humidity [kg/kg]
     169REAL,     INTENT(IN)   , DIMENSION(klon) :: qia_in        ! input linear contrails total specific humidity [kg/kg]
    166170REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_dist   ! aviation distance flown concentration [m/s/m3]
    167171REAL,     INTENT(IN)   , DIMENSION(klon) :: flight_h2o    ! aviation emitted H2O concentration [kgH2O/s/m3]
    … …  
    186190REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_con  ! cloud fraction tendency because of condensation [s-1]
    187191REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_mix  ! cloud fraction tendency because of cloud mixing [s-1]
     192REAL,     INTENT(INOUT), DIMENSION(klon) :: dcf_sed  ! cloud fraction tendency because of sedimentation [s-1]
    188193REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_adj  ! specific ice content tendency because of temperature adjustment [kg/kg/s]
    189194REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_sub  ! specific ice content tendency because of sublimation [kg/kg/s]
    190195REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_con  ! specific ice content tendency because of condensation [kg/kg/s]
    191196REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_mix  ! specific ice content tendency because of cloud mixing [kg/kg/s]
     197REAL,     INTENT(INOUT), DIMENSION(klon) :: dqi_sed  ! specific ice content tendency because of sedimentation [kg/kg/s]
    192198REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_adj ! specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
    193199REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sub ! specific cloud water vapor tendency because of sublimation [kg/kg/s]
    194200REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_con ! specific cloud water vapor tendency because of condensation [kg/kg/s]
    195201REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_mix ! specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
     202REAL,     INTENT(INOUT), DIMENSION(klon) :: dqvc_sed ! specific cloud water vapor tendency because of sedimentation [kg/kg/s]
    196203!
    197204!  Diagnostics for aviation
    … …  
    199206!
    200207REAL,     INTENT(INOUT), DIMENSION(klon) :: contfra      ! linear contrail fraction [-]
     208REAL,     INTENT(INOUT), DIMENSION(klon) :: perscontfra  ! linear contrail induced cirrus fraction [-]
    201209REAL,     INTENT(INOUT), DIMENSION(klon) :: qcont        ! linear contrail specific humidity [kg/kg]
    202210REAL,     INTENT(INOUT), DIMENSION(klon) :: Tcritcont    ! critical temperature for contrail formation [K]
    … …  
    220228INTEGER :: i
    221229LOGICAL :: ok_warm_cloud
    222 REAL, DIMENSION(klon) :: qtot, qcld, qzero
     230REAL, DIMENSION(klon) :: qcld, qzero
    223231REAL, DIMENSION(klon) :: clrfra, qclr
    224232!
    … …  
    228236!
    229237! for unadjusted clouds
    230 REAL :: qvapincld, qiceincld, qvapincld_new
    231 !
    232 ! for sublimation
    233 REAL :: dqt_sub
     238REAL :: qiceincld, qvapincld, qvapincld_new
     239!
     240! for deposition / sublimation
     241REAL :: pres_sat, kappa_depsub, tau_depsub
     242REAL :: air_thermal_conduct, water_vapor_diff
     243REAL :: N_ice
     244!
     245! for dissipation
     246REAL :: pdf_shape
    234247!
    235248! for condensation
    236249REAL, DIMENSION(klon) :: qsatl, dqsatl
    237 REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_loc
    238 REAL :: rhl_clr
    239 REAL :: pdf_ratqs, pdf_skew
    240 REAL :: pdf_x, pdf_y, pdf_T
    241 REAL :: pdf_e3, pdf_e4
    242 REAL :: cf_cond, qt_cond, dqt_con
     250REAL, DIMENSION(klon) :: pdf_alpha, pdf_scale, pdf_gamma
     251REAL :: rhl_clr, pdf_loc
     252REAL :: pdf_e3, pdf_x, pdf_y
     253REAL :: dqt_con
     254!
     255! for sedimentation
     256REAL, DIMENSION(klon) :: qice_sedim
     257REAL :: clrfra_sed
    243258!
    244259! for mixing
    245260REAL :: a_mix, bovera, Povera, N_cld_mix, L_mix
    246 REAL :: clrfra_mix, sigma_mix
     261REAL :: cldfra_mix, clrfra_mix, sigma_mix
    247262REAL :: L_shear, shear_fra
    248 REAL :: qvapinclr_lim
     263REAL :: qiceinmix, qvapinmix_lim, qvapinclr_lim
    249264REAL :: pdf_fra_above_nuc, pdf_q_above_nuc
    250265REAL :: pdf_fra_above_lim, pdf_q_above_lim
    251 REAL :: pdf_fra_below_lim
     266REAL :: pdf_fra_below_lim, pdf_q_below_lim
     267REAL :: mixed_fraction
    252268!
    253269! for cell properties
    254 REAL, DIMENSION(klon) :: dz
     270REAL :: rho, rhodz, dz
     271!
     272! for contrails
     273REAL :: perscontfra_ratio
     274REAL :: contrails_conversion_factor
    255275
    256276qzero(:) = 0.
    257 qtot = qtot_in
    258277
    259278!--Calculation of qsat w.r.t. liquid
    … …  
    269288  !--formed elsewhere)
    270289  IF (keepgoing(i)) THEN
    271 
    272     !--Initialisation
    273     issrfra(i)  = 0.
    274     qissr(i)    = 0.
    275     contfra(i)  = 0.
    276     qcont(i)    = 0.
    277290
    278291    !--If the temperature is higher than the threshold below which
    … …  
    281294    !--If ok_weibull_warm_clouds = .TRUE., the Weibull law is used for
    282295    !--all clouds, and the lognormal scheme is not activated
    283     IF ( ( ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds ) .OR. &
    284          ( ok_no_issr_strato .AND. ( stratomask(i) .EQ. 1. ) ) ) THEN
    285 
    286       pdf_std   = ratqs(i) * qtot(i)
    287       pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) )
    288       pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) )
     296    IF ( .NOT. pt_pron_clds(i) ) THEN
     297
     298      pdf_std   = ratqs(i) * qtot_in(i)
     299      pdf_k     = -SQRT( LOG( 1. + (pdf_std / qtot_in(i))**2 ) )
     300      pdf_delta = LOG( qtot_in(i) / ( gamma_cond(i) * qsat(i) ) )
    289301      pdf_a     = pdf_delta / ( pdf_k * SQRT(2.) )
    290302      pdf_b     = pdf_k / (2. * SQRT(2.))
    … …  
    303315      ELSE
    304316        cldfra(i) = 0.5 * pdf_e1
    305         qincld(i) = qtot(i) * pdf_e2 / pdf_e1
     317        qincld(i) = qtot_in(i) * pdf_e2 / pdf_e1
    306318        qvc(i) = qsat(i) * cldfra(i)
    307319      ENDIF
    … …  
    319331      ok_warm_cloud = ( temp(i) .GT. temp_nowater )
    320332
    321       !--We remove the convection water from the total available water
    322       qtot(i) = qtot(i) - ( qcondcv(i) + qsat(i) ) * cldfracv(i)
    323 
    324333      !--The following barriers ensure that the traced cloud properties
    325334      !--are consistent. In some rare cases, i.e. the cloud water vapor
    326335      !--can be greater than the total water in the gridbox
    327       cldfra(i) = MAX(0., MIN(1. - cldfracv(i), cldfra_in(i)))
    328       qcld(i)   = MAX(0., MIN(qtot(i), qliq_in(i) + qice_in(i) + qvc_in(i)))
     336      cldfra(i) = MAX(0., MIN(1. - cfcon(i), cldfra_in(i)))
     337      qcld(i)   = MAX(0., MIN(qtot_in(i), qliq_in(i) + qice_in(i) + qvc_in(i)))
    329338      qvc(i)    = MAX(0., MIN(qcld(i), qvc_in(i)))
    330339
    331340      !--Initialise clear fraction properties
    332       clrfra(i) = ( 1. - cldfracv(i) ) - cldfra(i)
    333       qclr(i) = qtot(i) - qcld(i)
     341      clrfra(i) = MAX(0., MIN(1., ( 1. - cfcon(i) ) - cldfra(i)))
     342      qclr(i) = qtot_in(i) - qcld(i)
    334343
    335344      dcf_sub(i)  = 0.
    … …  
    344353      dqi_mix(i)  = 0.
    345354      dqvc_mix(i) = 0.
     355      dcf_sed(i)  = 0.
     356      dqi_sed(i)  = 0.
     357      dqvc_sed(i) = 0.
     358
     359      IF ( icesed_flux(i) .GT. 0. ) THEN
     360        !--If ice sedimentation is activated, the quantity of sedimentated ice was added
     361        !--to the total water vapor in the precipitation routine. Here we remove it
     362        !--(it will be reincluded later)
     363        qice_sedim(i) = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
     364        qclr(i) = qclr(i) - qice_sedim(i)
     365      ENDIF
    346366
    347367      !--Initialisation of the cell properties
    348368      !--Dry density [kg/m3]
    349       !rho = pplay(i) / temp(i) / RD
     369      rho = pplay(i) / temp(i) / RD
    350370      !--Dry air mass [kg/m2]
    351       !rhodz = ( paprsdn(i) - paprsup(i) ) / RG
     371      rhodz = ( paprsdn(i) - paprsup(i) ) / RG
    352372      !--Cell thickness [m]
    353       !dz = rhodz / rho
    354       dz(i) = ( paprsdn(i) - paprsup(i) ) / pplay(i) * temp(i) * RD / RG
     373      dz = rhodz / rho
     374
     375      !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
     376      !--hypothesis is lost, and the vapor in the cloud is purely prognostic.
     377      !
     378      !--The deposition equation is
     379      !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
     380      !--from Lohmann et al. (2016), where
     381      !--alpha is the deposition coefficient [-]
     382      !--mi is the mass of one ice crystal [kg]
     383      !--C is the capacitance of an ice crystal [m]
     384      !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
     385      !--R_v is the specific gas constant for humid air [J/kg/K]
     386      !--T is the temperature [K]
     387      !--esi is the saturation pressure w.r.t. ice [Pa]
     388      !--Dv is the diffusivity of water vapor [m2/s]
     389      !--Ls is the specific latent heat of sublimation [J/kg/K]
     390      !--ka is the thermal conductivity of dry air [J/m/s/K]
     391      !
     392      !--alpha is a coefficient to take into account the fact that during deposition, a water
     393      !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
     394      !--coherent (with the same structure). It has no impact for sublimation.
     395      !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
     396      !--during deposition, and alpha = 1. during sublimation.
     397      !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
     398      !-- C = capa_cond_cirrus * rm_ice
     399      !
     400      !--We have qice = Nice * mi, where Nice is the ice crystal
     401      !--number concentration per kg of moist air
     402      !--HYPOTHESIS 1: the ice crystals are spherical, therefore
     403      !-- mi = 4/3 * pi * rm_ice**3 * rho_ice
     404      !--HYPOTHESIS 2: the ice crystals are monodisperse with the
     405      !--initial radius rm_ice_0.
     406      !--NB. this is notably different than the assumption
     407      !--of a distributed qice in the cloud made in the sublimation process;
     408      !--should it be consistent?
     409      !
     410      !--As the deposition process does not create new ice crystals,
     411      !--and because we assume a same rm_ice value for all crystals
     412      !--therefore the sublimation process does not destroy ice crystals
     413      !--(or, in a limit case, it destroys all ice crystals), then
     414      !--Nice is a constant during the sublimation/deposition process.
     415      !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
     416      !
     417      !--The deposition equation then reads:
     418      !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
     419      !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
     420      !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
     421      !--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
     422      !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
     423      !--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
     424      !--and we have
     425      !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
     426      !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
     427      !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
     428      !
     429      !--This system of equations can be resolved with an exact
     430      !--explicit numerical integration, having one variable resolved
     431      !--explicitly, the other exactly. The exactly resolved variable is
     432      !--the one decreasing, so it is qvc if the process is
     433      !--condensation, qi if it is sublimation.
     434      !
     435      !--kappa is computed as an initialisation constant, as it depends only
     436      !--on temperature and other pre-computed values
     437      pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i)
     438      !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
     439      air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184
     440      !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
     441      water_vapor_diff = 0.211 * ( temp(i) / RTT )**1.94 * ( 101325. / pplay(i) ) * 1.e-4
     442      !--Median ice crystal concentration [#/m3] in cirrus clouds from Kramer et al (2020)
     443      N_ice = 0.003 * 1e6
     444      !--NB. the greater kappa_depsub, the lower the efficiency of the
     445      !--deposition/sublimation process
     446      kappa_depsub = ( 4. / 3. * RPI * N_ice / rho * rho_ice )**(1./3.) &
     447                   * qsat(i) / ( 4. * RPI * capa_cond_cirrus * N_ice / rho ) &
     448                   * ( RV * temp(i) / water_vapor_diff / pres_sat &
     449                     + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )
    355450
    356451      !--If contrails are activated
    357452      IF ( ok_plane_contrail ) THEN
    358         contfra(i) = contfra_in(i)
    359         qcont(i) = qcont_in(i)
     453        contfra(i) = MAX(0., MIN(cldfra(i), contfra_in(i)))
     454        perscontfra(i) = MAX(0., MIN(cldfra(i), perscontfra_in(i)))
     455        qcont(i) = MAX(0., MIN(qcld(i), qva_in(i) + qia_in(i)))
    360456
    361457        dcfa_ini(i) = 0.
    … …  
    370466        dqia_mix(i) = 0.
    371467        dqta_mix(i) = 0.
     468      ELSE
     469        contfra(i) = 0.
     470        perscontfra(i) = 0.
     471        qcont(i) = 0.
    372472      ENDIF
    373473
    … …  
    379479      !--If there is a contrail
    380480      IF ( contfra(i) .GT. eps ) THEN
     481        !--We remove contrails from the main class
     482        cldfra(i) = cldfra(i) - contfra(i)
     483        qcld(i) = qcld(i) - qcont(i)
     484        qvc(i) = qvc(i) - qsat(i) * contfra(i)
     485
    381486        !--The contrail is always adjusted to saturation
    382487        qiceincld = ( qcont(i) / contfra(i) - qsat(i) )
    383488       
    384489        !--If the ice water content is too low, the cloud is purely sublimated
    385         IF ( qiceincld .LT. eps ) THEN
     490        IF ( qiceincld .LT. ( 0.005 * qsat(i) ) ) THEN
    386491          dcfa_sub(i) = - contfra(i)
    387492          dqia_sub(i) = - qiceincld * contfra(i)
    … …  
    389494          contfra(i) = 0.
    390495          qcont(i) = 0.
    391           clrfra(i) = clrfra(i) + dcfa_sub(i)
    392           qclr(i) = qclr(i) + dqta_sub(i)
    393         ELSE
    394           !--We remove contrails from the main class
    395           cldfra(i) = cldfra(i) - contfra(i)
    396           qcld(i) = qcld(i) - qcont(i)
    397           qvc(i) = qvc(i) - qsat(i) * contfra(i)
     496          clrfra(i) = MIN(1., clrfra(i) - dcfa_sub(i))
     497          qclr(i) = qclr(i) - dqta_sub(i)
    398498        ENDIF   ! qiceincld .LT. eps
    399499      ENDIF ! contfra(i) .GT. eps
     500
     501      !--If there is a contrail induced cirrus, we save the ratio
     502      perscontfra_ratio = perscontfra(i) / MAX(eps, cldfra(i))
    400503
    401504
    … …  
    421524          qcld(i) = 0.
    422525          qvc(i) = 0.
    423           clrfra(i) = clrfra(i) - dcf_sub(i)
     526          clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
    424527          qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
    425 
    426           !--If all the ice has been sublimated, we sublimate
    427           !--completely the cloud and do not activate the sublimation
    428           !--process
    429           qvapincld_new = 0.
    430528
    431529        !--Else, the cloud is adjusted and sublimated
    432530        ELSE
    433531
    434           !--The vapor in cloud cannot be higher than the
    435           !--condensation threshold
    436           qvapincld = MIN(qvapincld, gamma_cond(i) * qsat(i))
    437           qiceincld = ( qcld(i) / cldfra(i) - qvapincld )
     532          !IF ( qiceincld .LT. ( 0.005 * qsat(i) ) ) THEN
     533          !  dcf_sub(i) = ( qiceincld / ( 0.005 * qsat(i) ) - 1. ) * cldfra(i)
     534          !  dqvc_sub(i) = qvapincld * dcf_sub(i)
     535
     536          !  cldfra(i) = cldfra(i) + dcf_sub(i)
     537          !  qcld(i) = qcld(i) + dqvc_sub(i)
     538          !  qvc(i) = qvc(i) + dqvc_sub(i)
     539          !  clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
     540          !  qclr(i) = qclr(i) - dqvc_sub(i)
     541          !ENDIF
    438542
    439543          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    440             qvapincld_new = QVAPINCLD_DEPSUB( &
    441                 qvapincld, qiceincld, temp(i), qsat(i), pplay(i), dtime)
    442             IF ( qvapincld_new .EQ. 0. ) THEN
    443               !--If all the ice has been sublimated, we sublimate
    444               !--completely the cloud and do not activate the dissipation
    445               !--process
    446               !--Tendencies and diagnostics
    447               dcf_sub(i) = - cldfra(i)
    448               dqvc_sub(i) = - qvc(i)
    449               dqi_sub(i) = - ( qcld(i) - qvc(i) )
    450 
    451               cldfra(i) = 0.
    452               qcld(i) = 0.
    453               qvc(i) = 0.
    454               clrfra(i) = clrfra(i) - dcf_sub(i)
    455               qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
    456             ENDIF
     544            IF ( qvapincld .GE. qsat(i) ) THEN
     545              !--If the cloud is initially supersaturated
     546              !--Exact explicit integration (qvc exact, qice explicit)
     547              tau_depsub = 1. / ( depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub )
     548            ELSE
     549              !--If the cloud is initially subsaturated
     550              !--Exact explicit integration (qvc exact, qice explicit)
     551              !--Same but depo_coef_cirrus = 1
     552              tau_depsub = 1. / ( qiceincld**(1./3.) * kappa_depsub )
     553            ENDIF ! qvapincld .GT. qsat
     554            qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime / tau_depsub )
     555            !--If all the ice is sublimated
     556            IF ( qvapincld_new .GE. ( qvapincld + qiceincld ) ) qvapincld_new = 0.
    457557          ELSE
    458558            !--We keep the saturation adjustment hypothesis, and the vapor in the
    459559            !--cloud is set equal to the saturation vapor
    460             qvapincld_new = qsat(i)
     560            IF ( ( qvapincld + qiceincld ) .GT. qsat(i) ) THEN
     561              qvapincld_new = qsat(i)
     562            ELSE
     563              qvapincld_new = 0.
     564            ENDIF
    461565          ENDIF ! ok_unadjusted_clouds
    462          
    463           !--Adjustment of the IWC to the new vapor in cloud
    464           !--(this can be either positive or negative)
    465           dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
    466           dqi_adj(i) = - dqvc_adj(i)
    467 
    468           !--Add tendencies
    469           !--The vapor in the cloud is updated, but not qcld as it is constant
    470           !--through this process, as well as cldfra which is unmodified
    471           qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
    472566         
    473567
    … …  
    477571
    478572          !--If the dissipation process must be activated
    479           IF ( qvapincld_new .GT. 0. ) THEN
    480             !--Normal distribution around qvapincld + qiceincld with width
    481             !--constant in the RHi space
    482             pdf_y = ( qvapincld_new - qvapincld - qiceincld ) &
    483                   / ( std_subl_pdf_lscp / 100. * qsat(i))
    484             pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) )
    485             pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI )
    486 
    487             dcf_sub(i) = - cldfra(i) * pdf_e1
    488             dqt_sub    = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 &
    489                                        - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 )
     573          !IF ( cldfra(i) .GT. eps ) THEN
     574          IF ( qvapincld_new .GT. qvapincld ) THEN
     575            !--Gamma distribution starting at qvapincld
     576            pdf_shape = ( mu_subl_pdf_lscp + 1. ) / qiceincld
     577            pdf_y = pdf_shape * ( qvapincld_new - qvapincld )
     578            pdf_e1 = GAMMAINC ( mu_subl_pdf_lscp + 1. , pdf_y )
     579            pdf_e2 = GAMMAINC ( mu_subl_pdf_lscp + 2. , pdf_y )
    490580
    491581            !--Tendencies and diagnostics
    492             dqvc_sub(i) = dqt_sub
     582            dcf_sub(i)  = - cldfra(i) * pdf_e1
     583            dqi_sub(i)  = - cldfra(i) * pdf_e2 / pdf_shape
     584            dqvc_sub(i) = dcf_sub(i) * qvapincld
    493585
    494586            !--Add tendencies
    495             cldfra(i) = cldfra(i) + dcf_sub(i)
    496             qcld(i)   = qcld(i) + dqt_sub
     587            cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i))
     588            qcld(i)   = qcld(i) + dqvc_sub(i) + dqi_sub(i)
    497589            qvc(i)    = qvc(i) + dqvc_sub(i)
    498             clrfra(i) = clrfra(i) - dcf_sub(i)
    499             qclr(i)   = qclr(i) - dqt_sub
     590            clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
     591            qclr(i)   = qclr(i) - dqvc_sub(i) - dqi_sub(i)
    500592          ENDIF ! qvapincld_new .GT. 0.
     593
     594
     595          !------------------------------------
     596          !--       PHASE ADJUSTMENT        --
     597          !------------------------------------
     598
     599          IF ( qvapincld_new .EQ. 0. ) THEN
     600            !--If all the ice has been sublimated, we sublimate
     601            !--completely the cloud and do not activate the dissipation
     602            !--process
     603            !--Tendencies and diagnostics
     604            dcf_sub(i) = - cldfra(i)
     605            dqvc_sub(i) = - qvc(i)
     606            dqi_sub(i) = - ( qcld(i) - qvc(i) )
     607
     608            cldfra(i) = 0.
     609            qcld(i) = 0.
     610            qvc(i) = 0.
     611            clrfra(i) = MIN(1., clrfra(i) - dcf_sub(i))
     612            qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
     613          ELSE
     614            !--Adjustment of the IWC to the new vapor in cloud
     615            !--(this can be either positive or negative)
     616            dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )
     617            dqi_adj(i) = - dqvc_adj(i)
     618
     619            !--Add tendencies
     620            !--The vapor in the cloud is updated, but not qcld as it is constant
     621            !--through this process, as well as cldfra which is unmodified
     622            qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))
     623          ENDIF
    501624
    502625        ENDIF   ! qiceincld .LT. eps
    503626      ENDIF ! cldfra(i) .GT. eps
     627
     628      !--If there is a contrail induced cirrus, we restore it
     629      perscontfra(i) = perscontfra_ratio * cldfra(i)
    504630
    505631
    … …  
    519645        !--Calculation of the properties of the PDF
    520646        !--Parameterization from IAGOS observations
    521         !--pdf_e1 and pdf_e2 will be reused below
     647        !--pdf_alpha, pdf_scale and pdf_gamma will be reused below
    522648
    523649        !--Coefficient for standard deviation:
    … …  
    537663        pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3
    538664       
    539         IF ( ok_warm_cloud ) THEN
    540           !--If the statistical scheme is activated, the standard deviation is adapted
    541           !--to depend on the pressure level. It is multiplied by ratqs, so that near the
    542           !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
    543           pdf_std = pdf_std * ratqs(i)
    544         ENDIF
    545 
    546         pdf_e2 = GAMMA(1. + 1. / pdf_alpha(i))
    547         pdf_scale(i) = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha(i)) - pdf_e2**2 ))
    548         pdf_loc(i) = rhl_clr - pdf_scale(i) * pdf_e2
     665        !IF ( ok_warm_cloud ) THEN
     666        !  !--If the statistical scheme is activated, the standard deviation is adapted
     667        !  !--to depend on the pressure level. It is multiplied by ratqs, so that near the
     668        !  !--surface std is almost 0, and upper than about 450 hPa the std is left untouched
     669        !  pdf_std = pdf_std * ratqs(i)
     670        !ENDIF
     671
     672        pdf_gamma(i) = GAMMA(1. + 1. / pdf_alpha(i))
     673        pdf_scale(i) = MAX(eps, pdf_std / SQRT( &
     674                                    GAMMA(1. + 2. / pdf_alpha(i)) - pdf_gamma(i)**2 ))
     675        pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
    549676
    550677        !--Calculation of the newly condensed water and fraction (pronostic)
    … …  
    553680
    554681        pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
    555         pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     682        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    556683        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    557         pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
    558         cf_cond = EXP( - pdf_y )
    559         qt_cond = ( pdf_e3 + pdf_loc(i) * cf_cond ) * qsatl(i) / 100.
    560 
    561         IF ( cf_cond .GT. eps ) THEN
    562 
    563           dcf_con(i) = clrfra(i) * cf_cond
    564           dqt_con = clrfra(i) * qt_cond
     684        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
     685        pdf_fra_above_nuc = EXP( - pdf_y )
     686        pdf_q_above_nuc = ( pdf_e3 + pdf_loc * pdf_fra_above_nuc ) * qsatl(i) / 100.
     687
     688        IF ( pdf_fra_above_nuc .GT. eps ) THEN
     689
     690          dcf_con(i) = clrfra(i) * pdf_fra_above_nuc
     691          dqt_con = clrfra(i) * pdf_q_above_nuc
    565692
    566693          !--Barriers (should be useless
    … …  
    574701            qvapincld = gamma_cond(i) * qsat(i)
    575702            qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i)
    576             qvapincld_new = QVAPINCLD_DEPSUB( &
    577                 qvapincld, qiceincld, temp(i), qsat(i), pplay(i), dtime / 2.)
    578             qvapincld = qvapincld_new
     703            tau_depsub = 1. / ( depo_coef_cirrus * qiceincld**(1./3.) * kappa_depsub )
     704            qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime / 2. / tau_depsub )
    579705          ELSE
    580706            !--We keep the saturation adjustment hypothesis, and the vapor in the
    581707            !--newly formed cloud is set equal to the saturation vapor.
    582             qvapincld = qsat(i)
     708            qvapincld_new = qsat(i)
    583709          ENDIF
    584710
    585711          !--Tendency on cloud vapor and diagnostic
    586           dqvc_con(i) = qvapincld * dcf_con(i)
     712          dqvc_con(i) = qvapincld_new * dcf_con(i)
    587713          dqi_con(i) = dqt_con - dqvc_con(i)
    588714
    589715          !--Add tendencies
    590           cldfra(i)  = cldfra(i) + dcf_con(i)
     716          cldfra(i)  = MIN(1., cldfra(i) + dcf_con(i))
    591717          qcld(i)    = qcld(i) + dqt_con
    592718          qvc(i)     = qvc(i) + dqvc_con(i)
    593           clrfra(i)  = clrfra(i) - dcf_con(i)
     719          clrfra(i)  = MAX(0., clrfra(i) - dcf_con(i))
    594720          qclr(i)    = qclr(i) - dqt_con
    595721
    … …  
    597723      ENDIF ! ( 1. - cldfra(i) ) .GT. eps
    598724
     725
     726      !--If there is a contrail induced cirrus, we save the ratio
     727      perscontfra_ratio = perscontfra(i) / MAX(eps, cldfra(i))
    599728
    600729      !--------------------------------------
    … …  
    658787        clrfra_mix = N_cld_mix * RPI / cell_area(i) &
    659788                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
     789        !--The fraction of clear sky mixed is
     790        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
     791        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
     792                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
    660793
    661794
    … …  
    666799        !--semi-major axis (a_mix), on the entire cell heigh dz.
    667800        !--The increase in size is
    668         L_shear = coef_shear_lscp * shear(i) * dz(i) * dtime
     801        L_shear = coef_shear_lscp * shear(i) * dz * dtime
    669802        !--therefore, the fraction of clear sky mixed is
    670803        !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area
    … …  
    678811        !--mixed clouds are different.
    679812        clrfra_mix = clrfra_mix + shear_fra
    680 
     813        cldfra_mix = cldfra_mix + shear_fra
    681814
    682815        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
    683816
    684         clrfra_mix = MIN( clrfra_mix, clrfra(i) )
     817        clrfra_mix = MAX(eps, MIN(clrfra(i), clrfra_mix))
     818        cldfra_mix = MAX(eps, MIN(cldfra(i), cldfra_mix))
    685819
    686820        !--We compute the limit vapor in clear sky where the mixed cloud could not
    … …  
    690824        !--diagnostic, and if the cloud size is increased, we add the new vapor to the
    691825        !--cloud's vapor without condensing or sublimating ice crystals
    692         qvapinclr_lim = ( ( cldfra(i) + clrfra_mix ) * qsat(i) - qcld(i) ) / clrfra_mix
     826        IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     827          qiceinmix = ( qcld(i) - qvc(i) ) / cldfra(i) / ( 1. + clrfra_mix / cldfra_mix )
     828          tau_depsub = 1. / ( qiceinmix**(1./3.) * kappa_depsub )
     829          qvapinmix_lim = qsat(i) - qiceinmix / ( 1. - EXP( - dtime / tau_depsub ) )
     830          qvapinclr_lim = qvapinmix_lim * ( 1. + cldfra_mix / clrfra_mix ) &
     831                        - qvc(i) / cldfra(i) * cldfra_mix / clrfra_mix
     832        ELSE
     833          !--NB. if tau_depsub = 0, we get the same result as above
     834          qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
     835                        - qcld(i) / cldfra(i) * cldfra_mix / clrfra_mix
     836        ENDIF
    693837
    694838        IF ( qvapinclr_lim .LT. 0. ) THEN
    … …  
    696840          dcf_mix(i)  = clrfra_mix
    697841          dqvc_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
    698           dqi_mix(i)  = 0.
    699842        ELSE
    700843          !--We then calculate the clear sky part where the humidity is lower than
    … …  
    704847          !--because the clear sky fraction can only be reduced by condensation.
    705848          !--Thus the `pdf_xxx` variables are well defined.
    706           pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
    707           pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     849
     850          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     851          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
     852          pdf_x = qvapinclr_lim / qsatl(i) * 100.
     853          pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    708854          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    709           pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
    710           pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i)
    711           pdf_q_above_nuc = ( pdf_e3 * clrfra(i) &
    712                           + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100.
    713 
    714           pdf_x = qvapinclr_lim / qsatl(i) * 100.
    715           pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    716           pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    717           pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     855          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
    718856          pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
    719857          pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
    720                           + pdf_loc(i) * pdf_fra_above_lim ) * qsatl(i) / 100.
     858                          + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
    721859
    722860          pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
    723           pdf_fra_above_lim = pdf_fra_above_lim - pdf_fra_above_nuc
    724           pdf_q_above_lim = pdf_q_above_lim - pdf_q_above_nuc
     861          pdf_q_below_lim = qclr(i) - pdf_q_above_lim
    725862
    726863          !--sigma_mix is the ratio of the surroundings of the clouds where mixing
    … …  
    733870            dcf_mix(i)  = clrfra_mix * sigma_mix
    734871            dqvc_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
    735             dqi_mix(i)  = 0.
    736872          ENDIF
    737873
    738874          IF ( pdf_fra_below_lim .GT. eps ) THEN
    739             dcf_mix(i)  = dcf_mix(i)  - cldfra(i) * ( 1. - sigma_mix )
    740             dqvc_mix(i) = dqvc_mix(i) - cldfra(i) * ( 1. - sigma_mix ) &
    741                                                   * qvc(i) / cldfra(i)
    742             dqi_mix(i)  = dqi_mix(i)  - cldfra(i) * ( 1. - sigma_mix ) &
    743                                                   * ( qcld(i) - qvc(i) ) / cldfra(i)
     875            dcf_mix(i)  = dcf_mix(i)  - cldfra_mix * ( 1. - sigma_mix )
     876            dqvc_mix(i) = dqvc_mix(i) - cldfra_mix * ( 1. - sigma_mix ) &
     877                                      * qvc(i) / cldfra(i)
     878            dqi_mix(i)  = dqi_mix(i)  - cldfra_mix * ( 1. - sigma_mix ) &
     879                                      * ( qcld(i) - qvc(i) ) / cldfra(i)
    744880          ENDIF
     881
    745882        ENDIF
    746 
    747         !--Add tendencies
    748         cldfra(i)  = cldfra(i) + dcf_mix(i)
    749         qcld(i)    = qcld(i) + dqvc_mix(i) + dqi_mix(i)
    750         qvc(i)     = qvc(i) + dqvc_mix(i)
    751         clrfra(i)  = clrfra(i) - dcf_mix(i)
    752         qclr(i)    = qclr(i) - dqvc_mix(i) - dqi_mix(i)
    753      
    754883      ENDIF ! ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )
    755884
    … …  
    812941        clrfra_mix = N_cld_mix * RPI / cell_area(i) &
    813942                   * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 )
     943        !--The fraction of clear sky mixed is
     944        !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area
     945        cldfra_mix = N_cld_mix * RPI / cell_area(i) &
     946                   * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 )
    814947       
    815948       
    … …  
    823956        !--With this, the clouds increase in size along b only, by a factor
    824957        !--L_shear * SQRT(2.) / 2. (to account for the 45 degrees orientation of the wind)
    825         L_shear = coef_shear_contrails * shear(i) * dz(i) * dtime
     958        L_shear = coef_shear_contrails * shear(i) * dz * dtime
    826959        !--therefore, the fraction of clear sky mixed is
    827960        !-- N_cld_mix * ( a * (b + L_shear * SQRT(2.) / 2.) - a * b ) * RPI / 2. / cell_area
    … …  
    835968        !--mixed clouds are different.
    836969        clrfra_mix = clrfra_mix + shear_fra
     970        cldfra_mix = cldfra_mix + shear_fra
    837971       
    838972       
    839973        !-- PART 3 - CALCULATION OF THE MIXING PROPERTIES
    840974       
    841         clrfra_mix = MIN( clrfra_mix, clrfra(i) )
     975        clrfra_mix = MAX(eps, MIN(clrfra(i), clrfra_mix))
     976        cldfra_mix = MAX(eps, MIN(cldfra(i), cldfra_mix))
    842977       
    843978        !--We compute the limit vapor in clear sky where the mixed cloud could not
    … …  
    847982        !--diagnostic, and if the cloud size is increased, we add the new vapor to the
    848983        !--cloud's vapor without condensing or sublimating ice crystals
    849         qvapinclr_lim = ( ( contfra(i) + clrfra_mix ) * qsat(i) - qcont(i) ) / clrfra_mix
     984        qvapinclr_lim = qsat(i) * ( 1. + cldfra_mix / clrfra_mix ) &
     985                      - qcont(i) / contfra(i) * cldfra_mix / clrfra_mix
    850986
    851987        IF ( qvapinclr_lim .LT. 0. ) THEN
    … …  
    853989          dcfa_mix(i) = clrfra_mix
    854990          dqta_mix(i) = clrfra_mix * qclr(i) / clrfra(i)
    855           dqia_mix(i) = 0.
    856991        ELSE
    857992          !--We then calculate the clear sky part where the humidity is lower than
    … …  
    861996          !--because the clear sky fraction can only be reduced by condensation.
    862997          !--Thus the `pdf_xxx` variables are well defined.
    863           pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
    864           pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     998
     999          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     1000          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
     1001          pdf_x = qvapinclr_lim / qsatl(i) * 100.
     1002          pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    8651003          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    866           pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
    867           pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i)
    868           pdf_q_above_nuc = ( pdf_e3 * clrfra(i) &
    869                           + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100.
    870 
    871           pdf_x = qvapinclr_lim / qsatl(i) * 100.
    872           pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    873           pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    874           pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     1004          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
    8751005          pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
    8761006          pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
    877                           + pdf_loc(i) * pdf_fra_above_lim ) * qsatl(i) / 100.
     1007                          + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
    8781008
    8791009          pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
    880           pdf_fra_above_lim = pdf_fra_above_lim - pdf_fra_above_nuc
    881           pdf_q_above_lim = pdf_q_above_lim - pdf_q_above_nuc
    8821010
    8831011          !--sigma_mix is the ratio of the surroundings of the clouds where mixing
    … …  
    8931021            dcfa_mix(i) = clrfra_mix * sigma_mix
    8941022            dqta_mix(i) = clrfra_mix * sigma_mix * pdf_q_above_lim / pdf_fra_above_lim
    895             dqia_mix(i) = 0.
    8961023          ENDIF
    8971024
    8981025          IF ( pdf_fra_below_lim .GT. eps ) THEN
    899             dcfa_mix(i) = dcfa_mix(i) - contfra(i) * ( 1. - sigma_mix )
    900             dqta_mix(i) = dqta_mix(i) - contfra(i) * ( 1. - sigma_mix ) &
    901                                       * qcont(i) / contfra(i)
    902             dqia_mix(i) = dqia_mix(i) - contfra(i) * ( 1. - sigma_mix ) &
    903                                       * ( qcont(i) / contfra(i) - qsat(i) )
     1026            qvapincld = qcont(i) / contfra(i)
     1027            qiceincld = qvapincld - qsat(i)
     1028            dcfa_mix(i) = dcfa_mix(i) - cldfra_mix * ( 1. - sigma_mix )
     1029            dqta_mix(i) = dqta_mix(i) - cldfra_mix * ( 1. - sigma_mix ) * qvapincld
     1030            dqia_mix(i) = dqia_mix(i) - cldfra_mix * ( 1. - sigma_mix ) * qiceincld
    9041031          ENDIF
     1032
     1033        ENDIF
     1034      ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
     1035
     1036      IF ( contfra(i) .GT. eps ) THEN
     1037        !--We balance the mixing tendencies between the different cloud classes
     1038        mixed_fraction = dcf_mix(i) + dcfa_mix(i)
     1039        IF ( mixed_fraction .GT. clrfra(i) ) THEN
     1040          mixed_fraction = clrfra(i) / mixed_fraction
     1041          dcf_mix(i)  = dcf_mix(i)  * mixed_fraction
     1042          dqvc_mix(i) = dqvc_mix(i) * mixed_fraction
     1043          dqi_mix(i)  = dqi_mix(i)  * mixed_fraction
     1044          dcfa_mix(i) = dcfa_mix(i) * mixed_fraction
     1045          dqta_mix(i) = dqta_mix(i) * mixed_fraction
    9051046        ENDIF
    9061047
    … …  
    9101051        clrfra(i)  = clrfra(i) - dcfa_mix(i)
    9111052        qclr(i)    = qclr(i) - dqta_mix(i)
    912 
    913       ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
     1053      ENDIF
     1054
     1055      !--Add tendencies
     1056      cldfra(i)  = cldfra(i) + dcf_mix(i)
     1057      qcld(i)    = qcld(i) + dqvc_mix(i) + dqi_mix(i)
     1058      qvc(i)     = qvc(i) + dqvc_mix(i)
     1059      clrfra(i)  = clrfra(i) - dcf_mix(i)
     1060      qclr(i)    = qclr(i) - dqvc_mix(i) - dqi_mix(i)
     1061
     1062      !--If there is a contrail induced cirrus, we restore it
     1063      perscontfra(i) = perscontfra_ratio * cldfra(i)
     1064
     1065
     1066      !---------------------------------------
     1067      !--         ICE SEDIMENTATION         --
     1068      !---------------------------------------
     1069      !
     1070      !--If ice supersaturation is activated, the cloud properties are prognostic.
     1071      !--The falling ice is then considered a new cloud in the gridbox.
     1072      !--BEWARE with this parameterisation, we can create a new cloud with a much
     1073      !--different ice water content and water vapor content than the existing cloud
     1074      !--(if it exists). This implies than unphysical fluxes of ice and vapor
     1075      !--occur between the existing cloud and the newly formed cloud (from sedimentation).
     1076      !--Note also that currently, the precipitation scheme assume a maximum
     1077      !--random overlap, meaning that the new formed clouds will not be affected
     1078      !--by vertical wind shear.
     1079      !
     1080      IF ( icesed_flux(i) .GT. 0. ) THEN
     1081
     1082        clrfra_sed = MIN(clrfra(i), cldfra_above(i) - cldfra(i))
     1083
     1084        IF ( ( clrfra_sed .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN
     1085
     1086          qiceincld = qice_sedim(i) / cldfra_above(i)
     1087          IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     1088            tau_depsub = 1. / ( qiceincld**(1./3.) * kappa_depsub )
     1089            qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime / tau_depsub ) )
     1090          ELSE
     1091            qvapinclr_lim = qsat(i) - qiceincld
     1092          ENDIF
     1093
     1094          rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     1095          pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
     1096          pdf_x = qvapinclr_lim / qsatl(i) * 100.
     1097          pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     1098          pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     1099          pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
     1100          pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
     1101          pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
     1102                          + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
     1103
     1104          IF ( pdf_fra_above_lim .GT. eps ) THEN
     1105            dcf_sed(i)  = clrfra_sed * pdf_fra_above_lim / clrfra(i)
     1106            dqvc_sed(i) = dcf_sed(i) * pdf_q_above_lim / pdf_fra_above_lim
     1107            !--We include the sedimentated ice:
     1108            dqi_sed(i)  = qiceincld &    !--We include the sedimentated ice:
     1109                        * ( dcf_sed(i) & !--the part that falls in the newly formed cloud and
     1110                          + cldfra(i) )  !--the part that falls in the existing cloud
     1111                       
     1112          ENDIF
     1113
     1114        ELSE
     1115
     1116          dqi_sed(i) = qice_sedim(i)
     1117
     1118        ENDIF ! clrfra(i) .GT. eps
     1119
     1120        !--Add tendencies
     1121        cldfra(i)  = MIN(1., cldfra(i) + dcf_sed(i))
     1122        qcld(i)    = qcld(i) + dqvc_sed(i) + dqi_sed(i)
     1123        qvc(i)     = qvc(i) + dqvc_sed(i)
     1124        clrfra(i)  = MAX(0., clrfra(i) - dcf_sed(i))
     1125        !--We re-include sublimated sedimentated ice into clear sky water vapor
     1126        qclr(i)    = qclr(i) - dqvc_sed(i) + qice_sedim(i) - dqi_sed(i)
     1127
     1128      ENDIF ! qice_sedim(i) .GT. 0.
     1129
    9141130
    9151131      !--We put back contrails in the clouds class
    … …  
    9191135
    9201136
    921       !--Diagnostics
    922       dcf_sub(i)  = dcf_sub(i)  / dtime
    923       dcf_con(i)  = dcf_con(i)  / dtime
    924       dcf_mix(i)  = dcf_mix(i)  / dtime
    925       dqi_adj(i)  = dqi_adj(i)  / dtime
    926       dqi_sub(i)  = dqi_sub(i)  / dtime
    927       dqi_con(i)  = dqi_con(i)  / dtime
    928       dqi_mix(i)  = dqi_mix(i)  / dtime
    929       dqvc_adj(i) = dqvc_adj(i) / dtime
    930       dqvc_sub(i) = dqvc_sub(i) / dtime
    931       dqvc_con(i) = dqvc_con(i) / dtime
    932       dqvc_mix(i) = dqvc_mix(i) / dtime
    933        
    9341137      !--Diagnose ISSRs
    9351138      IF ( clrfra(i) .GT. eps ) THEN
    936         !--We then calculate the part that is greater than qsat
    937         !--and lower than gamma_cond * qsat, which is the ice supersaturated region
    938         pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100.
    939         pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
     1139        rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     1140        pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
     1141        pdf_x = qsat(i) / qsatl(i) * 100.
     1142        pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    9401143        pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    941         pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
    942         pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i)
    943         pdf_q_above_nuc = ( pdf_e3 * clrfra(i) &
    944                         + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100.
    945 
    946         pdf_x = qsat(i) / qsatl(i) * 100.
    947         pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)
    948         pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    949         pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2
     1144        pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
    9501145        issrfra(i) = EXP( - pdf_y ) * clrfra(i)
    951         qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc(i) * issrfra(i) ) * qsatl(i) / 100.
    952 
    953         issrfra(i) = issrfra(i) - pdf_fra_above_nuc
    954         qissr(i) = qissr(i) - pdf_q_above_nuc
     1146        qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc * issrfra(i) ) * qsatl(i) / 100.
     1147      ELSE
     1148        issrfra(i) = 0.
     1149        qissr(i) = 0.
    9551150      ENDIF
    9561151
    … …  
    9691164      ENDIF ! cldfra .LT. eps
    9701165
     1166      !--Diagnostics
     1167      dcf_sub(i)  = dcf_sub(i)  / dtime
     1168      dcf_con(i)  = dcf_con(i)  / dtime
     1169      dcf_mix(i)  = dcf_mix(i)  / dtime
     1170      dcf_sed(i)  = dcf_sed(i)  / dtime
     1171      dqi_adj(i)  = dqi_adj(i)  / dtime
     1172      dqi_sub(i)  = dqi_sub(i)  / dtime
     1173      dqi_con(i)  = dqi_con(i)  / dtime
     1174      dqi_mix(i)  = dqi_mix(i)  / dtime
     1175      dqi_sed(i)  = dqi_sed(i)  / dtime
     1176      dqvc_adj(i) = dqvc_adj(i) / dtime
     1177      dqvc_sub(i) = dqvc_sub(i) / dtime
     1178      dqvc_con(i) = dqvc_con(i) / dtime
     1179      dqvc_mix(i) = dqvc_mix(i) / dtime
     1180      dqvc_sed(i) = dqvc_sed(i) / dtime
     1181
    9711182    ENDIF ! ( temp(i) .GT. temp_nowater ) .AND. .NOT. ok_weibull_warm_clouds
    9721183
    … …  
    9831194  CALL contrails_formation( &
    9841195      klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, &
    985       flight_dist, clrfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &
     1196      flight_dist, clrfra, qclr, pdf_scale, pdf_alpha, pdf_gamma, &
     1197      keepgoing, pt_pron_clds, &
    9861198      Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    9871199      dcfa_ini, dqia_ini, dqta_ini)
    9881200
    9891201  DO i = 1, klon
    990     IF ( keepgoing(i) .AND. ( temp(i) .LE. temp_nowater ) ) THEN
     1202    IF ( keepgoing(i) .AND. pt_pron_clds(i) ) THEN
    9911203
    9921204      !--Convert existing contrail fraction into "natural" cirrus cloud fraction
    993       dcfa_cir(i) = contfra(i) * ( EXP( - dtime / linear_contrails_lifetime ) - 1. )
    994       dqta_cir(i) = qcont(i)   * ( EXP( - dtime / linear_contrails_lifetime ) - 1. )
     1205      IF ( cldfra(i) .GE. ( 1. - cfcon(i) - eps ) ) THEN
     1206        contrails_conversion_factor = 1.
     1207      ELSE
     1208        contrails_conversion_factor = ( 1. - &
     1209            !--First, the linear contrails are continuously degraded in induced cirrus
     1210            EXP( - dtime / linear_contrails_lifetime ) &
     1211            !--Second, if the cloud fills the entire gridbox, the linear contrails
     1212            !--cannot exist. The exponent is set so that this only happens for
     1213            !--very cloudy gridboxes
     1214            * ( 1. - cldfra(i) / ( 1. - cfcon(i) ) )**0.1 )
     1215      ENDIF
     1216      dcfa_cir(i) = - contrails_conversion_factor * contfra(i)
     1217      dqta_cir(i) = - contrails_conversion_factor * qcont(i)
    9951218
    9961219      !--Add tendencies
    9971220      issrfra(i) = MAX(0., issrfra(i) - dcfa_ini(i))
    9981221      qissr(i)   = MAX(0., qissr(i) - dqta_ini(i))
    999       cldfra(i)  = MAX(0., MIN(1. - cldfracv(i), cldfra(i) + dcfa_ini(i)))
    1000       qcld(i)    = MAX(0., MIN(qtot(i), qcld(i) + dqta_ini(i)))
     1222      clrfra(i)  = MAX(0., clrfra(i) - dcfa_ini(i))
     1223      qclr(i)    = MAX(0., qclr(i) - dqta_ini(i))
     1224
     1225      cldfra(i)  = MAX(0., MIN(1. - cfcon(i), cldfra(i) + dcfa_ini(i)))
     1226      qcld(i)    = MAX(0., MIN(qtot_in(i), qcld(i) + dqta_ini(i)))
    10011227      qvc(i)     = MAX(0., MIN(qcld(i), qvc(i) + dcfa_ini(i) * qsat(i)))
    10021228      contfra(i) = MAX(0., MIN(cldfra(i), contfra(i) + dcfa_cir(i) + dcfa_ini(i)))
    10031229      qcont(i)   = MAX(0., MIN(qcld(i), qcont(i) + dqta_cir(i) + dqta_ini(i)))
     1230      perscontfra(i) = perscontfra(i) - dcfa_cir(i)
    10041231
    10051232      !--Diagnostics
    … …  
    10241251        cldfra(i) = 0.
    10251252        contfra(i)= 0.
     1253        perscontfra(i) = 0.
    10261254        qcld(i)   = 0.
    10271255        qvc(i)    = 0.
     1256        qcont(i)  = 0.
    10281257        qincld(i) = qsat(i)
    10291258      ELSE
    … …  
    10361265      ENDIF
    10371266
     1267      IF ( perscontfra(i) .LT. eps ) perscontfra(i) = 0.
     1268
    10381269    ENDIF ! keepgoing
    10391270  ENDDO
    … …  
    10441275!**********************************************************************************
    10451276
    1046 FUNCTION QVAPINCLD_DEPSUB( &
    1047     qvapincld, qiceincld, temp, qsat, pplay, dtime)
    1048 
    1049 USE lmdz_lscp_ini, ONLY: RV, RLSTT, RTT, RD, EPS_W
    1050 USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice
    1051 USE lmdz_lscp_ini, ONLY: eps
    1052 
    1053 IMPLICIT NONE
    1054 
    1055 ! inputs
    1056 REAL :: qvapincld     !
    1057 REAL :: qiceincld     !
    1058 REAL :: temp          !
    1059 REAL :: qsat          !
    1060 REAL :: pplay         !
    1061 REAL :: dtime         ! time step [s]
    1062 ! outpu
    1063 REAL :: qvapincld_depsub
    1064 
    1065 !
    1066 ! local
    1067 REAL :: pres_sat, rho, kappa
    1068 REAL :: air_thermal_conduct, water_vapor_diff
    1069 REAL :: iwc
    1070 REAL :: iwc_log_inf100, iwc_log_sup100
    1071 REAL :: iwc_inf100, alpha_inf100, coef_inf100
    1072 REAL :: mu_sup100, sigma_sup100, coef_sup100
    1073 REAL :: Dm_ice, rm_ice
    1074 
    1075 !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment
    1076 !--hypothesis is lost, and the vapor in the cloud is purely prognostic.
    1077 !
    1078 !--The deposition equation is
    1079 !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )
    1080 !--from Lohmann et al. (2016), where
    1081 !--alpha is the deposition coefficient [-]
    1082 !--mi is the mass of one ice crystal [kg]
    1083 !--C is the capacitance of an ice crystal [m]
    1084 !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]
    1085 !--R_v is the specific gas constant for humid air [J/kg/K]
    1086 !--T is the temperature [K]
    1087 !--esi is the saturation pressure w.r.t. ice [Pa]
    1088 !--Dv is the diffusivity of water vapor [m2/s]
    1089 !--Ls is the specific latent heat of sublimation [J/kg/K]
    1090 !--ka is the thermal conductivity of dry air [J/m/s/K]
    1091 !
    1092 !--alpha is a coefficient to take into account the fact that during deposition, a water
    1093 !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays
    1094 !--coherent (with the same structure). It has no impact for sublimation.
    1095 !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))
    1096 !--during deposition, and alpha = 1. during sublimation.
    1097 !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus
    1098 !-- C = capa_cond_cirrus * rm_ice
    1099 !
    1100 !--We have qice = Nice * mi, where Nice is the ice crystal
    1101 !--number concentration per kg of moist air
    1102 !--HYPOTHESIS 1: the ice crystals are spherical, therefore
    1103 !-- mi = 4/3 * pi * rm_ice**3 * rho_ice
    1104 !--HYPOTHESIS 2: the ice crystals are monodisperse with the
    1105 !--initial radius rm_ice_0.
    1106 !--NB. this is notably different than the assumption
    1107 !--of a distributed qice in the cloud made in the sublimation process;
    1108 !--should it be consistent?
    1109 !
    1110 !--As the deposition process does not create new ice crystals,
    1111 !--and because we assume a same rm_ice value for all crystals
    1112 !--therefore the sublimation process does not destroy ice crystals
    1113 !--(or, in a limit case, it destroys all ice crystals), then
    1114 !--Nice is a constant during the sublimation/deposition process.
    1115 !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
    1116 !
    1117 !--The deposition equation then reads:
    1118 !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice
    1119 !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &
    1120 !--             / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &
    1121 !--                         * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )
    1122 !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &
    1123 !--  *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )
    1124 !--and we have
    1125 !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
    1126 !-- dqi/dt  =   qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2
    1127 !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))
    1128 !
    1129 !--This system of equations can be resolved with an exact
    1130 !--explicit numerical integration, having one variable resolved
    1131 !--explicitly, the other exactly. The exactly resolved variable is
    1132 !--the one decreasing, so it is qvc if the process is
    1133 !--condensation, qi if it is sublimation.
    1134 !
    1135 !--kappa is computed as an initialisation constant, as it depends only
    1136 !--on temperature and other pre-computed values
    1137 pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay
    1138 !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)
    1139 air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184
    1140 !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)
    1141 water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4
    1142 kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat &
    1143       * ( RV * temp / water_vapor_diff / pres_sat &
    1144         + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) )
    1145 !--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process
    1146 
    1147 !--Dry density [kg/m3]
    1148 rho = pplay / temp / RD
    1149 
    1150 !--Here, the initial vapor in the cloud is qvapincld, and we compute
    1151 !--the new vapor qvapincld_new
    1152 
    1153 !--rm_ice formula from McFarquhar and Heymsfield (1997)
    1154 iwc = qiceincld * rho * 1e3
    1155 iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837)
    1156 iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) )
    1157 iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) )
    1158 
    1159 alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100
    1160 coef_inf100 = iwc_inf100 * alpha_inf100**3 / 120.
    1161 
    1162 mu_sup100 = ( 5.2 + 0.0013 * ( temp - RTT ) ) &
    1163           + ( 0.026 - 1.2E-3 * ( temp - RTT ) ) * iwc_log_sup100
    1164 sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp - RTT ) ) &
    1165              + ( 0.018 - 2.1E-4 * ( temp - RTT ) ) * iwc_log_sup100
    1166 coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2 )
    1167 
    1168 Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2 ) &
    1169        * coef_sup100 ) / ( coef_inf100 + coef_sup100 )
    1170 rm_ice = Dm_ice / 2. * 1.E-6
    1171 
    1172 IF ( qvapincld .GE. qsat ) THEN
    1173   !--If the cloud is initially supersaturated
    1174   !--Exact explicit integration (qvc exact, qice explicit)
    1175   qvapincld_depsub = qsat + ( qvapincld - qsat ) &
    1176                    * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 )
    1177 ELSE
    1178   !--If the cloud is initially subsaturated
    1179   !--Exact explicit integration (qvc exact, qice explicit)
    1180   !--Same but depo_coef_cirrus = 1
    1181   qvapincld_depsub = qsat + ( qvapincld - qsat ) &
    1182                    * EXP( - dtime * qiceincld / kappa / rm_ice**2 )
    1183   IF ( qvapincld_depsub .GE. ( qvapincld + qiceincld ) ) &
    1184                    qvapincld_depsub = 0.
    1185 ENDIF ! qvapincld .GT. qsat
    1186 
    1187 END FUNCTION QVAPINCLD_DEPSUB
    1188 
    11891277END MODULE lmdz_lscp_condensation

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

    r5602 r5609  
    162162  !$OMP THREADPRIVATE(ok_weibull_warm_clouds)
    163163
     164  REAL, SAVE, PROTECTED :: ffallv_issr                       ! tuning coefficient crystal fall velocity, cirrus clouds (with ISSR)
     165  !$OMP THREADPRIVATE(ffallv_issr)
     166
    164167  REAL, SAVE, PROTECTED :: depo_coef_cirrus=.7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
    165168  !$OMP THREADPRIVATE(depo_coef_cirrus)
    … …  
    168171  !$OMP THREADPRIVATE(capa_cond_cirrus)
    169172
    170   REAL, SAVE, PROTECTED :: std_subl_pdf_lscp=2.              ! [%] standard deviation of the gaussian distribution of water inside ice clouds
    171   !$OMP THREADPRIVATE(std_subl_pdf_lscp)
    172  
    173   REAL, SAVE, PROTECTED :: beta_pdf_lscp=1.E-3               ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
     173  REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp=1./3.            ! [-] factor for the ice distribution inside cirrus clouds
     174  !$OMP THREADPRIVATE(mu_subl_pdf_lscp)
     175 
     176  REAL, SAVE, PROTECTED :: beta_pdf_lscp=8.75E-4             ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
    174177  !$OMP THREADPRIVATE(beta_pdf_lscp)
    175178 
    … …  
    192195  !$OMP THREADPRIVATE(b_homofreez)
    193196
    194   REAL, SAVE, PROTECTED :: delta_hetfreez=1.                 ! [-] value between 0 and 1 to simulate for heterogeneous freezing.
     197  REAL, SAVE, PROTECTED :: delta_hetfreez=0.85               ! [-] value between 0 and 1 to simulate for heterogeneous freezing.
    195198  !$OMP THREADPRIVATE(delta_hetfreez)
    196199 
    … …  
    344347  !$OMP THREADPRIVATE(ok_ice_sedim)
    345348
    346   REAL, SAVE, PROTECTED :: ice_fallspeed_sedim=0.1          ! Ice fallspeed velocity for sedimentation [m/s]
    347   !$OMP THREADPRIVATE(ice_fallspeed_sedim)
     349  REAL, SAVE, PROTECTED :: fallice_sedim=1.                 ! Tuning factor for ice fallspeed velocity for sedimentation [-]
     350  !$OMP THREADPRIVATE(fallice_sedim)
    348351  !--End of the parameters for poprecip
    349352
    … …  
    466469    CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld)
    467470    CALL getin_p('ok_ice_sedim',ok_ice_sedim)
    468     CALL getin_p('ice_fallspeed_sedim',ice_fallspeed_sedim)
     471    CALL getin_p('fallice_sedim',fallice_sedim)
    469472    ! for condensation and ice supersaturation
    470473    CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
    471474    CALL getin_p('ok_weibull_warm_clouds',ok_weibull_warm_clouds)
     475    ffallv_issr=ffallv_lsc
     476    CALL getin_p('ffallv_issr',ffallv_issr)
    472477    CALL getin_p('depo_coef_cirrus',depo_coef_cirrus)
    473478    CALL getin_p('capa_cond_cirrus',capa_cond_cirrus)
    474     CALL getin_p('std_subl_pdf_lscp',std_subl_pdf_lscp)
     479    CALL getin_p('mu_subl_pdf_lscp',mu_subl_pdf_lscp)
    475480    CALL getin_p('beta_pdf_lscp',beta_pdf_lscp)
    476481    CALL getin_p('temp_thresh_pdf_lscp',temp_thresh_pdf_lscp)
    … …  
    562567    WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
    563568    WRITE(lunout,*) 'lscp_ini, ok_ice_sedim:', ok_ice_sedim
    564     WRITE(lunout,*) 'lscp_ini, ice_fallspeed_sedim:', ice_fallspeed_sedim
     569    WRITE(lunout,*) 'lscp_ini, fallice_sedim:', fallice_sedim
    565570    ! for condensation and ice supersaturation
    566571    WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
    567572    WRITE(lunout,*) 'lscp_ini, ok_unadjusted_clouds:', ok_unadjusted_clouds
    568573    WRITE(lunout,*) 'lscp_ini, ok_weibull_warm_clouds:', ok_weibull_warm_clouds
     574    WRITE(lunout,*) 'lscp_ini, ffallv_issr', ffallv_issr
    569575    WRITE(lunout,*) 'lscp_ini, depo_coef_cirrus:', depo_coef_cirrus
    570576    WRITE(lunout,*) 'lscp_ini, capa_cond_cirrus:', capa_cond_cirrus
    571     WRITE(lunout,*) 'lscp_ini, std_subl_pdf_lscp:', std_subl_pdf_lscp
     577    WRITE(lunout,*) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
    572578    WRITE(lunout,*) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
    573579    WRITE(lunout,*) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp

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

    r5589 r5609  
    1818SUBROUTINE histprecip_precld( &
    1919           klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, &
    20            zmqc, zneb, znebprecip, znebprecipclr, &
     20           zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
    2121           zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
    2222           )
    2323
    24 USE lmdz_lscp_ini, ONLY : iflag_evap_prec
     24USE lmdz_lscp_ini, ONLY : iflag_evap_prec, ok_ice_sedim
    2525USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, ztfondue
    2626USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
    … …  
    4444REAL,    INTENT(INOUT), DIMENSION(klon) :: zq             !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg]
    4545REAL,    INTENT(INOUT), DIMENSION(klon) :: zmqc           !--specific humidity in the precipitation falling from the upper layer [kg/kg]
     46REAL,    INTENT(IN),    DIMENSION(klon) :: icesed_flux    !--sedimentated ice flux [kg/s/m2]
    4647
    4748REAL,    INTENT(IN),    DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
    … …  
    7475
    7576REAL :: zmelt
     77REAL :: qice_sedim
    7678INTEGER :: i
    7779
    … …  
    118120  ENDDO
    119121
     122ENDIF
     123
     124!--------------------
     125!-- ICE SEDIMENTATION
     126!--------------------
     127IF ( ok_ice_sedim ) THEN
     128  DO i =1, klon
     129    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
     130    !--added to the total water content of the gridbox
     131    IF ( icesed_flux(i) .GT. 0. ) THEN
     132      qice_sedim = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
     133      !--Vapor is updated after evaporation/sublimation (it is increased)
     134      zq(i) = zq(i) + qice_sedim
     135      !--Air and precip temperature (i.e., gridbox temperature)
     136      !--is updated due to latent heat cooling
     137      zt(i) = zt(i) &
     138              - qice_sedim * RLSTT / RCPD &
     139              / ( 1. + RVTMP2 * ( zq(i) + zmqc(i) ) )
     140    ENDIF ! ok_ice_sedim
     141  ENDDO
    120142ENDIF
    121143
    … …  
    307329
    308330SUBROUTINE histprecip_postcld( &
    309         klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, zdqsdT_raw, &
    310         zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
     331        klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, pt_pron_clds, &
     332        zdqsdT_raw, zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
    311333        rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
    312334        zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
    313         zradocond, zradoice, dqrauto, dqsauto &
     335        zradocond, zradoice, dqrauto, dqsauto, dqised &
    314336        )
    315337
    … …  
    321343                          iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, &
    322344                          niter_lscp
     345USE lmdz_lscp_ini, ONLY : ok_ice_sedim, fallice_sedim, eps
    323346USE lmdz_lscp_tools, ONLY : fallice_velocity
    324347
    … …  
    335358REAL,    INTENT(IN),    DIMENSION(klon) :: ctot_vol       !--volumic cloud fraction [-]
    336359LOGICAL, INTENT(IN),    DIMENSION(klon) :: ptconv         !--true if we are in a convective point
     360LOGICAL, INTENT(IN),    DIMENSION(klon) :: pt_pron_clds   !--true if we are in a prognostic clouds point
    337361REAL,    INTENT(IN),    DIMENSION(klon) :: zdqsdT_raw     !--derivative of qsat w.r.t. temperature times L/Cp [SI]
    338362
    … …  
    360384REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
    361385REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
     386REAL,    INTENT(OUT),   DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
    362387
    363388REAL,    INTENT(OUT),   DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
    … …  
    365390REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
    366391REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
     392REAL,    INTENT(INOUT), DIMENSION(klon) :: dqised         !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
    367393
    368394
    … …  
    390416REAL, DIMENSION(klon) :: ziflprev
    391417
     418! Local variable for sedimentation
     419REAL :: iwcg, tempc, icevelo
     420REAL :: qice_sedim
     421
    392422! Misc
    393423INTEGER :: i, n
    … …  
    398428zqpreci(:) = 0.
    399429ziflprev(:) = 0.
     430
     431icesed_flux(:) = 0.
    400432
    401433
    … …  
    474506  ENDDO
    475507
    476   CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, zvelo)
     508  CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, pt_pron_clds, zvelo)
    477509
    478510  DO i = 1, klon
    … …  
    656688ENDDO
    657689
     690!---------------------------------------------------------
     691!--                  ICE SEDIMENTATION
     692!---------------------------------------------------------
     693IF ( ok_ice_sedim ) THEN
     694
     695  DO i = 1, klon
     696    IF ( ( zoliqi(i) .GT. 0. ) .AND. ( rneb(i) .GT. eps ) .AND. pt_pron_clds(i) ) THEN
     697
     698      iwcg = MAX( zrho(i) * zoliqi(i) / zneb(i) * 1000., eps ) ! iwc in g/m3
     699      tempc = zt(i) - RTT ! celcius temp
     700      ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
     701      IF ( tempc .GE. -60. ) THEN
     702        icevelo = EXP( 1.9714 + 0.00216078 * tempc &
     703                - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) )
     704      ELSE
     705        icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15
     706      ENDIF
     707      icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
     708
     709      qice_sedim = zoliqi(i) * MIN(1., icevelo * dtime / zdz(i))
     710      !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two
     711      !--times in the same timestep)
     712      !qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)
     713      !--Add tendencies
     714      zoliqi(i) = zoliqi(i) - qice_sedim
     715      zoliq(i)  = zoliq(i)  - qice_sedim
     716      !--Convert to flux
     717      icesed_flux(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
     718
     719      !--Diagnostic tendencies
     720      dqised(i) = dqised(i) - qice_sedim / dtime
     721    ENDIF
     722  ENDDO
     723ENDIF
     724
    658725! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
    659726! if iflag_evap_prec>=4
    … …  
    713780SUBROUTINE poprecip_precld( &
    714781           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
    715            qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, qice_sedim, &
     782           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, icesed_flux, &
    716783           cldfra, qvc, qliq, qice, &
    717784           rain, rainclr, raincld, snow, snowclr, snowcld, &
    718            dqreva, dqssub, dcfsed, dqvcsed &
    719            )
     785           dqreva, dqssub)
    720786
    721787USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
    … …  
    749815REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
    750816REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
    751 REAL,    INTENT(IN),    DIMENSION(klon) :: qice_sedim     !--ice water content that sedimentated from the layer above [kg/kg]
     817REAL,    INTENT(IN),    DIMENSION(klon) :: icesed_flux    !--sedimentated ice water flux [kg/s/m2]
    752818
    753819REAL,    INTENT(INOUT), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
    … …  
    766832REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
    767833REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
    768 REAL,    INTENT(OUT),   DIMENSION(klon) :: dcfsed         !--cloud fraction tendency due to sedimentation of ice crystals [s-1]
    769 REAL,    INTENT(OUT),   DIMENSION(klon) :: dqvcsed        !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
    770834
    771835
    … …  
    791855!--Specific heat constant
    792856REAL :: cpair, cpw
     857!--Specific ice water content sedimentated
     858REAL :: qice_sedim
    793859
    794860!--Initialisation
    … …  
    796862dqreva(:) = 0.
    797863dqssub(:) = 0.
    798 IF ( ok_ice_sedim .AND. ok_ice_supersat ) THEN
    799   dcfsed(:) = 0.
    800   dqvcsed(:) = 0.
    801 ENDIF
    802864
    803865!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
    … …  
    9431005    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
    9441006    !--added to the total water content of the gridbox
    945     IF ( ok_ice_sedim .AND. ( qice_sedim(i) .GT. 0. ) ) THEN
     1007    IF ( icesed_flux(i) .GT. 0. ) THEN
     1008      qice_sedim = icesed_flux(i) / dhum_to_dflux(i)
    9461009      !--Vapor is updated after evaporation/sublimation (it is increased)
    947       qvap(i) = qvap(i) + qice_sedim(i)
     1010      qvap(i) = qvap(i) + qice_sedim
    9481011      !--Air and precip temperature (i.e., gridbox temperature)
    9491012      !--is updated due to latent heat cooling
    9501013      temp(i) = temp(i) &
    951               - qice_sedim(i) * RLSTT / RCPD &
     1014              - qice_sedim * RLSTT / RCPD &
    9521015              / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) )
    953    
    954       IF ( ok_ice_supersat ) THEN
    955         !--If ice supersaturation is activated, the cloud properties are prognostic.
    956         !--The falling ice is then considered a new cloud in the gridbox.
    957         !--BEWARE with this parameterisation, we can create a new cloud with a much
    958         !--different ice water content and water vapor content than the existing cloud
    959         !--(if it exists). This implies than unphysical fluxes of ice and vapor
    960         !--occur between the existing cloud and the newly formed cloud (from sedimentation).
    961         !--Note also that currently, the precipitation scheme assume a maximum
    962         !--random overlap, meaning that the new formed clouds will not be affected
    963         !--by vertical wind shear.
    964         ! Maybe we should reduce dcfsed?
    965         dcfsed(i) = precipfracclr_tmp(i)
    966         dqvcsed(i) = qvapclr * dcfsed(i)
    967         !--Add tendencies
    968         cldfra(i) = cldfra(i) + dcfsed(i)
    969         qvc(i) = qvc(i) + dqvcsed(i)
    970         qice(i) = qice(i) + qice_sedim(i)
    971         !--Normalisation
    972         dcfsed(i) = dcfsed(i) / dtime
    973         dqvcsed(i) = dqvcsed(i) / dtime
    974       ENDIF
    9751016    ENDIF ! ok_ice_sedim
    9761017
    … …  
    13261367SUBROUTINE poprecip_postcld( &
    13271368           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
    1328            temp, qvap, qliq, qice, icefrac, cldfra, qice_sedim, &
     1369           temp, qvap, qliq, qice, icefrac, cldfra, icesed_flux, &
    13291370           precipfracclr, precipfraccld, &
    13301371           rain, rainclr, raincld, snow, snowclr, snowcld, &
    … …  
    13471388                          rain_fallspeed_clr, rain_fallspeed_cld,              &
    13481389                          snow_fallspeed_clr, snow_fallspeed_cld,              &
    1349                           ok_ice_sedim, ice_fallspeed_sedim
     1390                          ok_ice_sedim, fallice_sedim
    13501391
    13511392
    … …  
    13681409REAL,    INTENT(IN),    DIMENSION(klon) :: icefrac        !--ice fraction [-]
    13691410REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction [-]
    1370 REAL,    INTENT(INOUT), DIMENSION(klon) :: qice_sedim     !--quantity of sedimentated ice [kg/kg]
    13711411
    13721412REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
    … …  
    13821422REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
    13831423
     1424REAL,    INTENT(OUT),   DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
    13841425REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
    13851426REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
    … …  
    13931434REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
    13941435REAL,    INTENT(OUT),   DIMENSION(klon) :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
    1395 REAL,    INTENT(OUT),  DIMENSION(klon) :: dqised         !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
     1436REAL,    INTENT(INOUT), DIMENSION(klon) :: dqised         !--ice water content tendency due to sedimentation of ice crystals [kg/kg/s]
    13961437
    13971438
    … …  
    14361477
    14371478!--Ice sedimentation
    1438 REAL :: dz, qice_sedim_tmp
     1479REAL :: iwcg, tempc, icevelo
     1480REAL :: dz, qice_sedim
    14391481
    14401482
    … …  
    14431485
    14441486qzero(:)    = 0.
     1487icesed_flux(:) = 0.
    14451488
    14461489dqrcol(:)   = 0.
    … …  
    14531496dqrfreez(:) = 0.
    14541497dqsfreez(:) = 0.
    1455 IF ( ok_ice_sedim ) dqised(:) = 0.
    14561498
    14571499
    … …  
    19171959  !---------------------------------------------------------
    19181960  IF ( ok_ice_sedim ) THEN
    1919     qice_sedim_tmp = 0.
    1920 
    1921     IF ( temp(i) .LE. temp_nowater ) THEN
     1961
     1962    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. eps ) ) THEN
     1963
    19221964      rho = pplay(i) / temp(i) / RD
     1965      iwcg = MAX( rho * qice(i) / cldfra(i) * 1000., eps ) ! iwc in g/m3
     1966      tempc = temp(i) - RTT ! celcius temp
     1967      ! so-called 'empirical parameterization' in Stubenrauch et al. 2019
     1968      IF ( tempc .GE. -60. ) THEN
     1969        icevelo = EXP( 1.9714 + 0.00216078 * tempc &
     1970                - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) )
     1971      ELSE
     1972        icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15
     1973      ENDIF
     1974      icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
     1975
    19231976      dz = ( paprsdn(i) - paprsup(i) ) / RG / rho
    1924       qice_sedim_tmp = qice(i) * MIN(1., ice_fallspeed_sedim * dtime / dz)
     1977      qice_sedim = qice(i) * MIN(1., icevelo * dtime / dz)
    19251978      !--Add tendencies
    1926       qice(i) = qice(i) - qice_sedim_tmp
    1927     ENDIF
    1928 
    1929     !--Diagnostic tendencies
    1930     dqised(i) = ( qice_sedim(i) - qice_sedim_tmp ) / dtime
    1931     !--Save quantity that will be sedimentated in the layer below
    1932     qice_sedim(i) = qice_sedim_tmp
     1979      qice(i) = qice(i) - qice_sedim
     1980      !--Convert to flux
     1981      icesed_flux(i) = qice_sedim * dhum_to_dflux(i)
     1982
     1983      !--Diagnostic tendencies
     1984      dqised(i) = dqised(i) - qice_sedim / dtime
     1985    ENDIF
    19331986  ENDIF
    19341987

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

    r5577 r5609  
    66
    77!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    8 SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo)
     8SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,ptpronclds,velo)
    99
    1010    ! Ref:
    … …  
    1616    ! 3212–3234. https://doi.org/10.1029/2019MS001642
    1717   
    18     use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc
     18    use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc, ffallv_issr
    1919    use lmdz_lscp_ini, only: cice_velo, dice_velo
    20     use lmdz_lscp_ini, only: ok_bug_ice_fallspeed
     20    use lmdz_lscp_ini, only: ok_bug_ice_fallspeed, eps
    2121
    2222    IMPLICIT NONE
    … …  
    2828    REAL, INTENT(IN), DIMENSION(klon) :: pres      ! air pressure [Pa]
    2929    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv    ! convective point  [-]
     30    LOGICAL, INTENT(IN), DIMENSION(klon) :: ptpronclds! prognostic clouds point  [-]
    3031
    3132    REAL, INTENT(OUT), DIMENSION(klon) :: velo    ! fallspeed velocity of crystals [m/s]
    … …  
    4243        IF (ptconv(i)) THEN
    4344            fallv_tun=ffallv_con
     45        ELSEIF (ptpronclds(i)) THEN
     46            fallv_tun=ffallv_issr
    4447        ELSE
    4548            fallv_tun=ffallv_lsc
    … …  
    5154        ELSE
    5255            ! AB the threshold is way too high, we reduce it
    53             iwcg=MAX(iwc(i)*1000.,1E-10) ! iwc in g/m3. We set a minimum value to prevent from division by 0
     56            iwcg=MAX(iwc(i)*1000.,eps) ! iwc in g/m3. We set a minimum value to prevent from division by 0
    5457        ENDIF
    5558        phpa=pres(i)/100.    ! pressure in hPa

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

    r5601 r5609  
    1818  USE surface_data,     ONLY : type_ocean, version_ocean
    1919  USE phyetat0_get_mod, ONLY : phyetat0_get, phyetat0_srf
    20   USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, phys_tstep, &
     20  USE phys_state_var_mod, ONLY : ancien_ok, cwcon, clwcon, detr_therm, phys_tstep, &
    2121       qsol, fevap, z0m, z0h, agesno, &
    2222       du_gwd_rando, du_gwd_front, entr_therm, f0, fm_therm, &
    2323       falb_dir, falb_dif, prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien, &
    2424       ftsol, pbl_tke, pctsrf, q_ancien, ql_ancien, qs_ancien, qbs_ancien, &
    25        cf_ancien, qvc_ancien, cfa_ancien, qta_ancien, radpas, radsol, rain_fall, ratqs, &
    26        rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
     25       cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien, qva_ancien, qia_ancien, &
     26       radpas, radsol, rain_fall, &
     27       ratqs, rnebcon, rugoro, sig1, snow_fall, bs_fall, solaire_etat0, sollw, sollwdown, &
    2728       solsw, solswfdiff, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
    2829       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, beta_aridity, wake_fip, wake_pe, &
    … …  
    415416    ancien_ok=ancien_ok.AND.phyetat0_get(cf_ancien,"CFANCIEN","CFANCIEN",0.)
    416417    ancien_ok=ancien_ok.AND.phyetat0_get(qvc_ancien,"QVCANCIEN","QVCANCIEN",0.)
     418    found=phyetat0_get(cwcon,"CWCON","CWCON",0.)
    417419  ELSE
    418420    cf_ancien(:,:)=0.
    … …  
    423425  IF ( ok_plane_contrail ) THEN
    424426    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.)
     427    ancien_ok=ancien_ok.AND.phyetat0_get(pcf_ancien,"PCFANCIEN","PCFANCIEN",0.)
     428    ancien_ok=ancien_ok.AND.phyetat0_get(qva_ancien,"QVAANCIEN","QVAANCIEN",0.)
     429    ancien_ok=ancien_ok.AND.phyetat0_get(qia_ancien,"QIAANCIEN","QIAANCIEN",0.)
    426430  ELSE
    427431    cfa_ancien(:,:)=0.
    428     qta_ancien(:,:)=0.
     432    pcf_ancien(:,:)=0.
     433    qva_ancien(:,:)=0.
     434    qia_ancien(:,:)=0.
    429435  ENDIF
    430436
    … …  
    458464  IF ( ok_plane_contrail ) THEN
    459465    IF ( ( maxval(cfa_ancien).EQ.minval(cfa_ancien) ) .OR. &
    460          ( maxval(qta_ancien).EQ.minval(qta_ancien) ) ) THEN
     466         ( maxval(pcf_ancien).EQ.minval(pcf_ancien) ) .OR. &
     467         ( maxval(qva_ancien).EQ.minval(qva_ancien) ) .OR. &
     468         ( maxval(qia_ancien).EQ.minval(qia_ancien) ) ) THEN
    461469       ancien_ok=.false.
    462470     ENDIF

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

    r5601 r5609  
    2323                                prw_ancien, prlw_ancien, prsw_ancien, prbsw_ancien,      &
    2424                                ql_ancien, qs_ancien, qbs_ancien, cf_ancien, &
    25                                 qvc_ancien, cfa_ancien, qta_ancien,          &
    26                                 u_ancien, v_ancien,                          &
    27                                 clwcon, rnebcon, ratqs, pbl_tke,             &
     25                                qvc_ancien, cfa_ancien, pcf_ancien,          &
     26                                qva_ancien, qia_ancien, u_ancien, v_ancien,  &
     27                                cwcon, clwcon, rnebcon, ratqs, pbl_tke,      &
    2828                                wake_delta_pbl_tke, zmax0, f0, sig1, w01,    &
    2929                                wake_deltat, wake_deltaq, wake_s, awake_s,   &
    … …  
    255255      CALL put_field(pass,"CFANCIEN", "CFANCIEN", cf_ancien)
    256256      CALL put_field(pass,"QVCANCIEN", "QVCANCIEN", qvc_ancien)
     257      CALL put_field(pass,"CWCON", "CWCON", cwcon)
    257258    ENDIF
    258259
    259260    IF ( ok_plane_contrail ) THEN
    260261      CALL put_field(pass,"CFAANCIEN", "CFAANCIEN", cfa_ancien)
    261       CALL put_field(pass,"QTAANCIEN", "QTAANCIEN", qta_ancien)
     262      CALL put_field(pass,"PCFANCIEN", "PCFANCIEN", pcf_ancien)
     263      CALL put_field(pass,"QVAANCIEN", "QVAANCIEN", qva_ancien)
     264      CALL put_field(pass,"QIAANCIEN", "QIAANCIEN", qia_ancien)
    262265    ENDIF
    263266

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

    r5601 r5609  
    669669      REAL, SAVE, ALLOCATABLE :: cfa_seri(:,:), d_cfa_dyn(:,:)
    670670      !$OMP THREADPRIVATE(cfa_seri, d_cfa_dyn)
    671       REAL, SAVE, ALLOCATABLE :: qta_seri(:,:), d_qta_dyn(:,:)
    672       !$OMP THREADPRIVATE(qta_seri, d_qta_dyn)
     671      REAL, SAVE, ALLOCATABLE :: pcf_seri(:,:), d_pcf_dyn(:,:)
     672      !$OMP THREADPRIVATE(pcf_seri, d_pcf_dyn)
     673      REAL, SAVE, ALLOCATABLE :: qva_seri(:,:), d_qva_dyn(:,:)
     674      !$OMP THREADPRIVATE(qva_seri, d_qva_dyn)
     675      REAL, SAVE, ALLOCATABLE :: qia_seri(:,:), d_qia_dyn(:,:)
     676      !$OMP THREADPRIVATE(qia_seri, d_qia_dyn)
    673677      REAL, SAVE, ALLOCATABLE :: flight_dist(:,:), flight_h2o(:,:)
    674678      !$OMP THREADPRIVATE(flight_dist, flight_h2o)
    … …  
    677681      REAL, SAVE, ALLOCATABLE :: potcontfraP(:,:), potcontfraNP(:,:)
    678682      !$OMP THREADPRIVATE(potcontfraP, potcontfraNP)
    679       REAL, SAVE, ALLOCATABLE :: qice_cont(:,:), dcfa_cir(:,:), dqta_cir(:,:)
    680       !$OMP THREADPRIVATE(qice_cont, dcfa_cir, dqta_cir)
     683      REAL, SAVE, ALLOCATABLE :: qice_perscont(:,:)
     684      !$OMP THREADPRIVATE(qice_perscont)
     685      REAL, SAVE, ALLOCATABLE :: dcfa_cir(:,:), dqta_cir(:,:)
     686      !$OMP THREADPRIVATE(dcfa_cir, dqta_cir)
    681687      REAL, SAVE, ALLOCATABLE :: dcfa_ini(:,:), dqia_ini(:,:), dqta_ini(:,:)
    682688      !$OMP THREADPRIVATE(dcfa_ini, dqia_ini, dqta_ini)
    … …  
    12441250      ALLOCATE(d_q_avi(klon,klev))
    12451251      ALLOCATE(cfa_seri(klon,klev), d_cfa_dyn(klon,klev))
    1246       ALLOCATE(qta_seri(klon,klev), d_qta_dyn(klon,klev))
     1252      ALLOCATE(pcf_seri(klon,klev), d_pcf_dyn(klon,klev))
     1253      ALLOCATE(qva_seri(klon,klev), d_qva_dyn(klon,klev))
     1254      ALLOCATE(qia_seri(klon,klev), d_qia_dyn(klon,klev))
    12471255      ALLOCATE(flight_dist(klon,klev), flight_h2o(klon,klev))
    12481256      ALLOCATE(Tcritcont(klon,klev), qcritcont(klon,klev))
    12491257      ALLOCATE(potcontfraP(klon,klev), potcontfraNP(klon,klev))
    1250       ALLOCATE(qice_cont(klon,klev), dcfa_cir(klon,klev), dqta_cir(klon,klev))
     1258      ALLOCATE(qice_perscont(klon,klev))
     1259      ALLOCATE(dcfa_cir(klon,klev), dqta_cir(klon,klev))
    12511260      ALLOCATE(dcfa_ini(klon,klev), dqia_ini(klon,klev), dqta_ini(klon,klev))
    12521261      ALLOCATE(dcfa_sub(klon,klev), dqia_sub(klon,klev), dqta_sub(klon,klev))
    … …  
    16611670
    16621671!-- LSCP - aviation and contrails variables
    1663       DEALLOCATE(d_q_avi, cfa_seri, d_cfa_dyn, qta_seri, d_qta_dyn, flight_dist, flight_h2o)
     1672      DEALLOCATE(cfa_seri, d_cfa_dyn, pcf_seri, d_pcf_dyn)
     1673      DEALLOCATE(qva_seri, d_qva_dyn, qia_seri, d_qia_dyn)
     1674      DEALLOCATE(d_q_avi, flight_dist, flight_h2o)
    16641675      DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
    1665       DEALLOCATE(qice_cont, dcfa_cir, dqta_cir, dcfa_ini, dqia_ini, dqta_ini)
     1676      DEALLOCATE(qice_perscont, dcfa_cir, dqta_cir, dcfa_ini, dqia_ini, dqta_ini)
    16661677      DEALLOCATE(dcfa_sub, dqia_sub, dqta_sub, dcfa_mix, dqia_mix, dqta_mix)
    16671678      DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont, conttau, contemi)

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

    r5601 r5609  
    21922192  TYPE(ctrl_out), SAVE :: o_dcfadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    21932193    '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) /))
     2194  TYPE(ctrl_out), SAVE :: o_pcfseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2195    'pcfseri', 'Contrail induced cirrus fraction', '-', (/ ('',i=1,10) /))
     2196  TYPE(ctrl_out), SAVE :: o_dpcfdyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2197    'dpcfdyn', 'Dynamics contrail induced cirrus fraction tendency', 's-1', (/ ('',i=1,10) /))
     2198  TYPE(ctrl_out), SAVE :: o_qvaseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2199    'qvaseri', 'Linear contrail total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
     2200  TYPE(ctrl_out), SAVE :: o_dqvadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2201    'dqvadyn', 'Dynamics linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
     2202  TYPE(ctrl_out), SAVE :: o_qiaseri = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2203    'qiaseri', 'Linear contrail total specific humidity', 'kg/kg', (/ ('',i=1,10) /))
     2204  TYPE(ctrl_out), SAVE :: o_dqiadyn = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
     2205    'dqiadyn', 'Dynamics linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('',i=1,10) /))
    21982206  TYPE(ctrl_out), SAVE :: o_Tcritcont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    21992207    'Tcritcont', 'Temperature threshold for contrail formation', 'K', (/ ('',i=1,10) /))
    … …  
    22042212  TYPE(ctrl_out), SAVE :: o_potcontfraNP = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/),&
    22052213    'potcontfraNP', 'Potential non-persistent contrail fraction', '-', (/ ('', i=1,10)/))
    2206   TYPE(ctrl_out), SAVE :: o_qice_cont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    2207     'qice_cont', 'Linear contrails ice specific humidity', 'kg/kg', (/ ('', i=1, 10) /))
     2214  TYPE(ctrl_out), SAVE :: o_qice_perscont = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2215    'qice_perscont', 'Contrail induced cirrus ice specific humidity', 'kg/kg', (/ ('', i=1, 10) /))
    22082216  TYPE(ctrl_out), SAVE :: o_dcfaini = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    22092217    'dcfaini', 'Initial formation linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))

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

    r5601 r5609  
    229229         o_issrfra250to300, o_issrfra300to400, o_issrfra400to500, &
    230230!-- LSCP - aviation variables
    231          o_cfaseri, o_dcfadyn, o_qtaseri, o_dqtadyn, o_dqavi, &
     231         o_cfaseri, o_dcfadyn, o_pcfseri, o_dpcfdyn, &
     232         o_qvaseri, o_dqvadyn, o_qiaseri, o_dqiadyn, o_dqavi, &
    232233         o_Tcritcont, o_qcritcont, o_potcontfraP, o_potcontfraNP, &
    233          o_flight_dist, o_flight_h2o, o_qice_cont, o_dcfacir, o_dqtacir, &
     234         o_flight_dist, o_flight_h2o, o_qice_perscont, o_dcfacir, o_dqtacir, &
    234235         o_dcfaini, o_dqiaini, o_dqtaini, o_dcfasub, o_dqiasub, o_dqtasub, &
    235236         o_dcfamix, o_dqiamix, o_dqtamix, &
    … …  
    360361         issrfra100to150, issrfra150to200, issrfra200to250, &
    361362         issrfra250to300, issrfra300to400, issrfra400to500, &
    362          cfa_seri, d_cfa_dyn, qta_seri, d_qta_dyn, d_q_avi, &
     363         cfa_seri, d_cfa_dyn, pcf_seri, d_pcf_dyn, &
     364         qva_seri, d_qva_dyn, qia_seri, d_qia_dyn, d_q_avi, &
    363365         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    364          flight_dist, flight_h2o, qice_cont, dcfa_cir, dqta_cir, &
     366         flight_dist, flight_h2o, qice_perscont, dcfa_cir, dqta_cir, &
    365367         dcfa_ini, dqia_ini, dqta_ini, dcfa_sub, dqia_sub, dqta_sub, &
    366368         dcfa_mix, dqia_mix, dqta_mix, &
    … …  
    21052107           CALL histwrite_phy(o_dqsfreez, dqsfreez)
    21062108           CALL histwrite_phy(o_dqsrim, dqsrim)
    2107            IF ( ok_ice_sedim ) THEN
    2108             CALL histwrite_phy(o_dqised, dqised)
    2109             CALL histwrite_phy(o_dcfsed, dcfsed)
    2110             CALL histwrite_phy(o_dqvcsed, dqvcsed)
    2111            ENDIF
    21122109           ELSE
    21132110            CALL histwrite_phy(o_dqreva, dqreva)
    … …  
    21162113            CALL histwrite_phy(o_dqsauto, dqsauto)
    21172114           ENDIF
     2115           CALL histwrite_phy(o_qsatl, qsatliq)
     2116           CALL histwrite_phy(o_qsati, qsatice)
    21182117       ENDIF
    21192118
    … …  
    21412140         CALL histwrite_phy(o_dqvccon, dqvc_con)
    21422141         CALL histwrite_phy(o_dqvcmix, dqvc_mix)
    2143          CALL histwrite_phy(o_qsatl, qsatliq)
    2144          CALL histwrite_phy(o_qsati, qsatice)
    21452142         CALL histwrite_phy(o_issrfra100to150, issrfra100to150)
    21462143         CALL histwrite_phy(o_issrfra150to200, issrfra150to200)
    … …  
    21492146         CALL histwrite_phy(o_issrfra300to400, issrfra300to400)
    21502147         CALL histwrite_phy(o_issrfra400to500, issrfra400to500)
     2148         IF ( ok_ice_sedim ) THEN
     2149          CALL histwrite_phy(o_dqised, dqised)
     2150          CALL histwrite_phy(o_dcfsed, dcfsed)
     2151          CALL histwrite_phy(o_dqvcsed, dqvcsed)
     2152         ENDIF
    21512153       ENDIF
    21522154!-- LSCP - aviation variables
    … …  
    21542156         CALL histwrite_phy(o_cfaseri, cfa_seri)
    21552157         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)
     2158         CALL histwrite_phy(o_pcfseri, pcf_seri)
     2159         CALL histwrite_phy(o_dpcfdyn, d_pcf_dyn)
     2160         CALL histwrite_phy(o_qvaseri, qva_seri)
     2161         CALL histwrite_phy(o_dqvadyn, d_qva_dyn)
     2162         CALL histwrite_phy(o_qiaseri, qia_seri)
     2163         CALL histwrite_phy(o_dqiadyn, d_qia_dyn)
    21582164         CALL histwrite_phy(o_flight_dist, flight_dist)
    21592165         CALL histwrite_phy(o_Tcritcont, Tcritcont)
    … …  
    21612167         CALL histwrite_phy(o_potcontfraP, potcontfraP)
    21622168         CALL histwrite_phy(o_potcontfraNP, potcontfraNP)
    2163          CALL histwrite_phy(o_qice_cont, qice_cont)
     2169         CALL histwrite_phy(o_qice_perscont, qice_perscont)
    21642170         CALL histwrite_phy(o_dcfacir, dcfa_cir)
    21652171         CALL histwrite_phy(o_dqtacir, dqta_cir)
    … …  
    22192225       CALL histwrite_phy(o_dqlphy2d,  zx_tmp_fi2d)
    22202226
    2221        IF (nqo.EQ.3) THEN
     2227       IF (nqo .GE. 3) THEN
    22222228       CALL histwrite_phy(o_dqsphy,  d_qx(:,:,isol))
    22232229       IF (vars_defined) CALL water_int(klon,klev,d_qx(:,:,isol),zmasse,zx_tmp_fi2d)

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

    r5601 r5609  
    9494      REAL, ALLOCATABLE, SAVE :: cf_ancien(:,:), qvc_ancien(:,:)
    9595!$OMP THREADPRIVATE(cf_ancien, qvc_ancien)
    96       REAL, ALLOCATABLE, SAVE :: cfa_ancien(:,:), qta_ancien(:,:)
    97 !$OMP THREADPRIVATE(cfa_ancien, qta_ancien)
     96      REAL, ALLOCATABLE, SAVE :: cfa_ancien(:,:), pcf_ancien(:,:)
     97!$OMP THREADPRIVATE(cfa_ancien, pcf_ancien)
     98      REAL, ALLOCATABLE, SAVE :: qva_ancien(:,:), qia_ancien(:,:)
     99!$OMP THREADPRIVATE(qva_ancien, qia_ancien)
     100      REAL, ALLOCATABLE, SAVE :: cwcon(:,:), cwcon0(:,:)
     101!$OMP THREADPRIVATE(cwcon, cwcon0)
    98102!!! RomP >>>
    99103      REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
    … …  
    590594      ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
    591595      ALLOCATE(cf_ancien(klon,klev), qvc_ancien(klon,klev))
    592       ALLOCATE(cfa_ancien(klon,klev), qta_ancien(klon,klev))
     596      ALLOCATE(cfa_ancien(klon,klev), pcf_ancien(klon,klev))
     597      ALLOCATE(qva_ancien(klon,klev), qia_ancien(klon,klev))
     598      ALLOCATE(cwcon(klon,klev), cwcon0(klon,klev))
    593599!!! Rom P >>>
    594600      ALLOCATE(tr_ancien(klon,klev,nbtr))
    … …  
    817823      DEALLOCATE(qtc_cv,sigt_cv,detrain_cv,fm_cv)
    818824      DEALLOCATE(u_ancien, v_ancien)
    819       DEALLOCATE(cf_ancien, qvc_ancien, cfa_ancien, qta_ancien)
     825      DEALLOCATE(cf_ancien, qvc_ancien, cfa_ancien, pcf_ancien)
     826      DEALLOCATE(qva_ancien, qia_ancien)
     827      DEALLOCATE(cwcon, cwcon0)
    820828      DEALLOCATE(tr_ancien)                           !RomP
    821829      DEALLOCATE(ratqs, pbl_tke,coefh,coefm)

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

    r5602 r5609  
    328328       dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    329329       !-- LSCP - aviation and contrails variables
    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, &
     330       cfa_seri, pcf_seri, qva_seri, qia_seri, d_cfa_dyn, d_pcf_dyn, d_qva_dyn, d_qia_dyn, &
     331       d_q_avi, flight_dist, flight_h2o, qice_perscont, &
     332       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    333333       cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
    334334       conttau, contemi, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont, &
    … …  
    518518    !
    519519    ! indices de traceurs eau vapeur, liquide, glace, fraction nuageuse LS (optional), blowing snow (optional)
    520     INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, iqvc, icfa, iqta
    521 !$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, icfa, iqta)
     520    INTEGER,SAVE :: ivap, iliq, isol, ibs, icf, iqvc, icfa, ipcf, iqia, iqva
     521!$OMP THREADPRIVATE(ivap, iliq, isol, ibs, icf, iqvc, icfa, ipcf, iqia, iqva)
    522522    !
    523523    !
    … …  
    13631363       isol = strIdx(tracers(:)%name, addPhase('H2O', 's'))
    13641364       ibs  = strIdx(tracers(:)%name, addPhase('H2O', 'b'))
    1365        icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
    1366        iqvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
    1367        icfa = strIdx(tracers(:)%name, addPhase('H2O', 'a'))
    1368        iqta = strIdx(tracers(:)%name, addPhase('H2O', 'q'))
     1365       icf  = strIdx(tracers(:)%name, 'CLDFRA')
     1366       iqvc = strIdx(tracers(:)%name, 'CLDVAP_g')
     1367       icfa = strIdx(tracers(:)%name, 'CONTFRA')
     1368       ipcf = strIdx(tracers(:)%name, 'PERSCONTFRA')
     1369       iqva = strIdx(tracers(:)%name, 'CONTWATER_g')
     1370       iqia = strIdx(tracers(:)%name, 'CONTWATER_s')
    13691371!       CALL init_etat0_limit_unstruct
    13701372!       IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
    … …  
    14141416       ENDIF
    14151417
    1416        IF (ok_ice_supersat.AND.(nqo.LT.5)) THEN
    1417           WRITE (lunout, *) ' ok_ice_supersat=y requires 5 H2O tracers ', &
    1418                '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c) but nqo=', nqo, '. Might as well stop here.'
     1418       IF (ok_ice_supersat.AND.((icf.EQ.0).OR.(iqvc.EQ.0))) THEN
     1419          WRITE (lunout, *) ' ok_ice_supersat=y requires 5 tracers ', &
     1420               '(H2O_g, H2O_l, H2O_s, CLDFRA, CLDVAP_g) but not all are ', &
     1421               'provided. Might as well stop here.'
    14191422          abort_message='see above'
    14201423          CALL abort_physic(modname,abort_message,1)
    … …  
    14331436       ENDIF
    14341437
    1435        IF (ok_plane_contrail.AND.(nqo.LT.6)) THEN
    1436           WRITE (lunout, *) ' ok_plane_contrail=y requires 6 H2O tracers ', &
    1437               '(H2O_g, H2O_l, H2O_s, H2O_f, H2O_c, H2O_a) but nqo=', nqo, '. Might as well stop here.'
     1438       IF (ok_plane_contrail.AND.((icfa.EQ.0).OR.(iqva.EQ.0).OR.(iqia.EQ.0).OR.(ipcf.EQ.0))) THEN
     1439          WRITE (lunout, *) ' ok_plane_contrail=y requires 8 tracers ', &
     1440              '(H2O_g, H2O_l, H2O_s, CLDFRA, CLDVAP_g, CONTFRA, PERSCONTFRA, CONTWATER_s) ', &
     1441              'but not all are provided. Might as well stop here.'
    14381442          abort_message='see above'
    14391443          CALL abort_physic(modname,abort_message,1)
    … …  
    14461450       ENDIF
    14471451
    1448         IF (ok_bs) THEN
    1449          IF ((ok_ice_supersat.AND.nqo .LT.6).OR.(.NOT.ok_ice_supersat.AND.nqo.LT.4)) THEN
    1450              WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
    1451                                'but nqo=', nqo
    1452              abort_message='see above'
    1453              CALL abort_physic(modname,abort_message, 1)
    1454          ENDIF
    1455         ENDIF
     1452       IF (ok_bs.AND.(nqo.LT.4)) THEN
     1453          WRITE (lunout, *) 'activation of blowing snow needs a specific H2O tracer', &
     1454                            'but nqo=', nqo
     1455          abort_message='see above'
     1456          CALL abort_physic(modname,abort_message, 1)
     1457       ENDIF
    14561458
    14571459       Ncvpaseq1 = 0
    … …  
    16311633       rnebcon(:,:) = 0.0
    16321634       clwcon(:,:) = 0.0
     1635       !--AB for prognostic clouds
     1636       cwcon0(:,:) = 0.0
     1637       cwcon(:,:) = 0.0
    16331638
    16341639       !
    … …  
    24542459          q_seri(i,k)  = qx(i,k,ivap)
    24552460          ql_seri(i,k) = qx(i,k,iliq)
     2461          qs_seri(i,k) = 0.
    24562462          qbs_seri(i,k)= 0.
    24572463          cf_seri(i,k) = 0.
    24582464          qvc_seri(i,k)= 0.
    24592465          cfa_seri(i,k)= 0.
    2460           qta_seri(i,k)= 0.
     2466          pcf_seri(i,k)= 0.
     2467          qva_seri(i,k)= 0.
     2468          qia_seri(i,k)= 0.
    24612469          !CR: ATTENTION, on rajoute la variable glace
    2462           IF (nqo.EQ.2) THEN             !--vapour and liquid only
    2463              qs_seri(i,k) = 0.
    2464           ELSE IF (nqo.EQ.3) THEN        !--vapour, liquid and ice
     2470          IF (nqo .GE. 3) THEN        !--vapour, liquid and ice
    24652471             qs_seri(i,k) = qx(i,k,isol)
    2466           ELSE IF (nqo.GE.4) THEN        !--vapour, liquid, ice, blowing snow, cloud fraction and cloudy water vapor to total water vapor ratio
    2467              qs_seri(i,k) = qx(i,k,isol)
    2468              IF (ok_ice_supersat) THEN
    2469                cf_seri(i,k) = qx(i,k,icf)
    2470                qvc_seri(i,k) = qx(i,k,iqvc)
    2471              ENDIF
    2472              IF (ok_plane_contrail) THEN
    2473                cfa_seri(i,k) = qx(i,k,icfa)
    2474                qta_seri(i,k) = qx(i,k,iqta)
    2475              ENDIF
    2476              IF (ok_bs) THEN
    2477                qbs_seri(i,k)= qx(i,k,ibs)
    2478              ENDIF
     2472          ENDIF
     2473          IF (ok_bs) THEN
     2474             qbs_seri(i,k) = qx(i,k,ibs)
     2475          ENDIF
     2476          IF (ok_ice_supersat) THEN
     2477             cf_seri(i,k)  = qx(i,k,icf)
     2478             qvc_seri(i,k) = qx(i,k,iqvc)
     2479          ENDIF
     2480          IF (ok_plane_contrail) THEN
     2481             cfa_seri(i,k) = qx(i,k,icfa)
     2482             pcf_seri(i,k) = qx(i,k,ipcf)
     2483             qva_seri(i,k) = qx(i,k,iqva)
     2484             qia_seri(i,k) = qx(i,k,iqia)
    24792485          ENDIF
    24802486       ENDDO
    … …  
    25522558       d_qvc_dyn(:,:)= (qvc_seri(:,:)-qvc_ancien(:,:))/phys_tstep
    25532559       d_cfa_dyn(:,:)= (cfa_seri(:,:)-cfa_ancien(:,:))/phys_tstep
    2554        d_qta_dyn(:,:)= (qta_seri(:,:)-qta_ancien(:,:))/phys_tstep
     2560       d_pcf_dyn(:,:)= (pcf_seri(:,:)-pcf_ancien(:,:))/phys_tstep
     2561       d_qva_dyn(:,:)= (qva_seri(:,:)-qva_ancien(:,:))/phys_tstep
     2562       d_qia_dyn(:,:)= (qia_seri(:,:)-qia_ancien(:,:))/phys_tstep
    25552563       CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
    25562564       d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
    … …  
    25752583       d_qvc_dyn(:,:)= 0.0
    25762584       d_cfa_dyn(:,:)= 0.0
    2577        d_qta_dyn(:,:)= 0.0
     2585       d_pcf_dyn(:,:)= 0.0
     2586       d_qva_dyn(:,:)= 0.0
     2587       d_qia_dyn(:,:)= 0.0
    25782588       d_q_dyn2d(:)  = 0.0
    25792589       d_ql_dyn2d(:) = 0.0
    … …  
    27002710            qi_seri_lscp(i,k) = qs_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27012711            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) )
     2712            qva_seri(i,k) = qva_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
     2713            qia_seri(i,k) = qia_seri(i,k) / ( q_seri(i,k) + ql_seri(i,k) + qs_seri(i,k) )
    27032714          ELSE
    27042715            ql_seri_lscp(i,k) = 0.
    27052716            qi_seri_lscp(i,k) = 0.
    27062717            qvc_seri(i,k) = 0.
    2707             qta_seri(i,k) = 0.
     2718            qva_seri(i,k) = 0.
     2719            qia_seri(i,k) = 0.
    27082720          ENDIF
    27092721        ENDDO
    … …  
    39273939        DO i = 1, klon
    39283940          qvc_seri(i,k) = qvc_seri(i,k) * q_seri(i,k)
     3941          cwcon0(i,k) = zqsat(i,k)
    39293942        ENDDO
    39303943      ENDDO
    … …  
    39333946      DO k = 1, klev
    39343947        DO i = 1, klon
    3935           qta_seri(i,k) = qta_seri(i,k) * q_seri(i,k)
     3948          qva_seri(i,k) = qva_seri(i,k) * q_seri(i,k)
     3949          qia_seri(i,k) = qia_seri(i,k) * q_seri(i,k)
    39363950        ENDDO
    39373951      ENDDO
    … …  
    39453959
    39463960    CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, &
    3947          t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, &
    3948          ptconv, rnebcon0, clwcon0, ratqs, sigma_qtherm, &
     3961         t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, ratqs, sigma_qtherm, &
     3962         ptconv, rnebcon, cwcon, clwcon, rnebcon0, cwcon0, clwcon0, &
    39493963         d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, &
    39503964         pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb,  &
    … …  
    39603974         dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
    39613975         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    3962          cfa_seri, qta_seri, flight_dist, flight_h2o, &
    3963          qice_cont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
     3976         cfa_seri, pcf_seri, qva_seri, qia_seri, flight_dist, flight_h2o, &
     3977         qice_perscont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    39643978         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
    39653979         qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
    … …  
    39673981         dqised, dcfsed, dqvcsed)
    39683982
     3983    IF (ok_ice_supersat) cwcon(:,:) = cwcon0(:,:)
    39693984
    39703985    ELSE
    … …  
    45134528               zfice, dNovrN, ptconv, rnebcon, clwcon, &
    45144529               !--AB contrails
    4515                cfa_seri, qice_cont, cldfra_nocont, cldtau_nocont, cldemi_nocont, &
    4516                conttau, contemi, cldh_nocont, contcov, &
     4530               cfa_seri, pcf_seri, qia_seri, qice_perscont, cldfra_nocont, &
     4531               cldtau_nocont, cldemi_nocont, conttau, contemi, cldh_nocont, contcov, &
    45174532               fiwp_nocont, fiwc_nocont, ref_ice_nocont)
    45184533
    … …  
    57095724             d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
    57105725          ENDIF
    5711           !--ice_supersat: nqo=5, we add cloud fraction and cloudy water vapor to total water vapor ratio
    5712           IF (nqo.ge.5 .and. ok_ice_supersat) THEN
     5726          IF (ok_bs) THEN
     5727             d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
     5728          ENDIF
     5729          IF (ok_ice_supersat) THEN
    57135730             d_qx(i,k,icf) = ( cf_seri(i,k) - qx(i,k,icf) ) / phys_tstep
    57145731             d_qx(i,k,iqvc) = ( qvc_seri(i,k) - qx(i,k,iqvc) ) / phys_tstep
    5715              IF (nqo.ge.6 .and. ok_plane_contrail) THEN
    5716                d_qx(i,k,icfa) = ( cfa_seri(i,k) - qx(i,k,icfa) ) / phys_tstep
    5717                d_qx(i,k,iqta) = ( qta_seri(i,k) - qx(i,k,iqta) ) / phys_tstep
    5718              ENDIF
    57195732          ENDIF
    5720 
    5721            IF (nqo.ge.4 .and. ok_bs) THEN
    5722              d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
     5733          IF (ok_plane_contrail) THEN
     5734             d_qx(i,k,icfa) = ( cfa_seri(i,k) - qx(i,k,icfa) ) / phys_tstep
     5735             d_qx(i,k,ipcf) = ( pcf_seri(i,k) - qx(i,k,ipcf) ) / phys_tstep
     5736             d_qx(i,k,iqva) = ( qva_seri(i,k) - qx(i,k,iqva) ) / phys_tstep
     5737             d_qx(i,k,iqia) = ( qia_seri(i,k) - qx(i,k,iqia) ) / phys_tstep
    57235738          ENDIF
    5724 
    57255739       ENDDO
    57265740    ENDDO
    … …  
    57545768    qvc_ancien(:,:)= qvc_seri(:,:)
    57555769    cfa_ancien(:,:)= cfa_seri(:,:)
    5756     qta_ancien(:,:)= qta_seri(:,:)
     5770    pcf_ancien(:,:)= pcf_seri(:,:)
     5771    qva_ancien(:,:)= qva_seri(:,:)
     5772    qia_ancien(:,:)= qia_seri(:,:)
    57575773    CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
    57585774    CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
Note: See TracChangeset for help on using the changeset viewer.