Changeset 5728 for LMDZ6


Ignore:
Timestamp:
Jun 27, 2025, 7:30:06 PM (5 weeks ago)
Author:
aborella
Message:

Correction to rare floating point errors + added support for the sedimentation of contrails

Location:
LMDZ6/branches/contrails
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r5717 r5728  
    739739        <field id="pfraclr" long_name="LS precip fraction clear-sky part"    unit="-" />               
    740740        <field id="pfracld" long_name="LS precip fraction cloudy part"    unit="-" />
     741        <field id="distcltop"    long_name="distance to cloud top"    unit="m" />
     742        <field id="tempcltop"    long_name="temperature of cloud top" unit="K" />
    741743        <field id="qrainlsc" long_name="LS specific rain water content"    unit="kg/kg" />
    742744        <field id="qsnowlsc" long_name="LS specific snow water content"    unit="kg/kg" />
    … …  
    940942        <field id="dqilmix"    long_name="Mixing linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
    941943        <field id="dqtlmix"    long_name="Mixing linear contrail total specific humidity tendency"    unit="kg/kg/s" />
     944        <field id="dcflsed"    long_name="Ice sedimentation linear contrail fraction tendency"    unit="s-1" />
     945        <field id="dqilsed"    long_name="Ice sedimentation linear contrail ice specific humidity tendency"    unit="kg/kg/s" />
     946        <field id="dqtlsed"    long_name="Ice sedimentation linear contrail total specific humidity tendency"    unit="kg/kg/s" />
    942947        <field id="dcfcsub"    long_name="Sublimation contrail cirrus fraction tendency"    unit="s-1" />
    943948        <field id="dqicsub"    long_name="Sublimation contrail cirrus ice specific humidity tendency"    unit="kg/kg/s" />
    … …  
    946951        <field id="dqicmix"    long_name="Mixing contrail cirrus ice specific humidity tendency"    unit="kg/kg/s" />
    947952        <field id="dqtcmix"    long_name="Mixing contrail cirrus total specific humidity tendency"    unit="kg/kg/s" />
     953        <field id="dcfcsed"    long_name="Ice sedimentation contrail cirrus fraction tendency"    unit="s-1" />
     954        <field id="dqicsed"    long_name="Ice sedimentation contrail cirrus ice specific humidity tendency"    unit="kg/kg/s" />
     955        <field id="dqtcsed"    long_name="Ice sedimentation contrail cirrus total specific humidity tendency"    unit="kg/kg/s" />
    948956        <field id="flightdist" long_name="Aviation flown distance concentration"    unit="m/s/m3" />
    949957        <field id="flighth2o"  long_name="Aviation emitted H2O concentration"    unit="kg H2O/s/m3" />

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

    r5717 r5728  
    9797      cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, &
    9898      temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, pt_pron_clds, &
    99       dzsed_abv, flsed_abv, cfsed_abv_in, &
     99      dzsed_abv, flsed_abv, cfsed_abv, &
    100100      dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv, &
    101101      dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv, &
    … …  
    110110      dcfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, &
    111111      dcfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix, &
     112      dcfl_sed, dqtl_sed, dqil_sed, dcfc_sed, dqtc_sed, dqic_sed, &
    112113      dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix)
    113114
    … …  
    172173REAL,     INTENT(IN)   , DIMENSION(klon) :: dzsed_abv     ! sedimentated cloud height IN THE LAYER ABOVE [m]
    173174REAL,     INTENT(IN)   , DIMENSION(klon) :: flsed_abv     ! sedimentated ice flux FROM THE LAYER ABOVE [kg/s/m2]
    174 REAL,     INTENT(IN)   , DIMENSION(klon) :: cfsed_abv_in  ! cloud fraction IN THE LAYER ABOVE [-]
     175REAL,     INTENT(IN)   , DIMENSION(klon) :: cfsed_abv     ! cloud fraction IN THE LAYER ABOVE [-]
    175176REAL,     INTENT(IN)   , DIMENSION(klon) :: dzsed_lincont_abv   ! sedimentated linear contrails height IN THE LAYER ABOVE [m]
    176177REAL,     INTENT(IN)   , DIMENSION(klon) :: flsed_lincont_abv   ! sedimentated ice flux in linear contrails FROM THE LAYER ABOVE [kg/s/m2]
    … …  
    251252REAL,     INTENT(INOUT), DIMENSION(klon) :: dqil_mix     ! linear contrails ice specific humidity tendency because of mixing [kg/kg/s]
    252253REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtl_mix     ! linear contrails total specific humidity tendency because of mixing [kg/kg/s]
     254REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfl_sed     ! linear contrails cloud fraction tendency because of ice sedimentation [s-1]
     255REAL,     INTENT(INOUT), DIMENSION(klon) :: dqil_sed     ! linear contrails ice specific humidity tendency because of ice sedimentation [kg/kg/s]
     256REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtl_sed     ! linear contrails total specific humidity tendency because of ice sedimentation [kg/kg/s]
    253257REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_sub     ! contrail cirrus cloud fraction tendency because of sublimation [s-1]
    254258REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_sub     ! contrail cirrus ice specific humidity tendency because of sublimation [kg/kg/s]
    … …  
    257261REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_mix     ! contrail cirrus ice specific humidity tendency because of mixing [kg/kg/s]
    258262REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_mix     ! contrail cirrus total specific humidity tendency because of mixing [kg/kg/s]
     263REAL,     INTENT(INOUT), DIMENSION(klon) :: dcfc_sed     ! contrail cirrus cloud fraction tendency because of ice sedimentation [s-1]
     264REAL,     INTENT(INOUT), DIMENSION(klon) :: dqic_sed     ! contrail cirrus ice specific humidity tendency because of ice sedimentation [kg/kg/s]
     265REAL,     INTENT(INOUT), DIMENSION(klon) :: dqtc_sed     ! contrail cirrus total specific humidity tendency because of ice sedimentation [kg/kg/s]
    259266!
    260267! Local
    … …  
    289296!
    290297! for sedimentation
    291 REAL, DIMENSION(klon) :: cfsed_abv, qised_abv
     298REAL, DIMENSION(klon) :: qised_abv
    292299REAL :: iwc, icefall_velo, qice_sedim
    293 REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3
     300REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3, sedfra_tmp
     301REAL :: dcf_sed1, dcf_sed2, dqvc_sed1, dqvc_sed2, dqi_sed1, dqi_sed2
    294302!
    295303! for mixing
    … …  
    309317REAL :: contrails_conversion_factor
    310318REAL, DIMENSION(klon) :: qised_lincont_abv, qised_circont_abv
     319REAL :: dcfl_sed1, dcfl_sed2, dqtl_sed1, dqtl_sed2, dqil_sed1, dqil_sed2
     320REAL :: dcfc_sed1, dcfc_sed2, dqtc_sed1, dqtc_sed2, dqic_sed1, dqic_sed2
    311321
    312322qzero(:) = 0.
    … …  
    401411
    402412      qised_abv(i) = 0.
    403       cfsed_abv(i) = cfsed_abv_in(i)
    404413      IF ( dzsed_abv(i) .GT. eps ) THEN
    405414        !--If ice sedimentation is activated, the quantity of sedimentated ice was added
    … …  
    514523        ENDIF
    515524
     525        qised_lincont_abv(i) = 0.
    516526        IF ( dzsed_lincont_abv(i) .GT. eps ) THEN
     527          !--If ice sedimentation is activated, the quantity of sedimentated ice was added
     528          !--to the total water vapor in the precipitation routine. Here we remove it
     529          !--(it will be reincluded later)
    517530          qised_lincont_abv(i) = flsed_lincont_abv(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
    518           qised_abv(i) = qised_abv(i) - qised_lincont_abv(i)
    519           cfsed_abv(i) = cfsed_abv(i) - cfsed_lincont_abv(i)
     531          qclr(i) = qclr(i) - qised_lincont_abv(i)
    520532        ENDIF
    521533
     534        qised_circont_abv(i) = 0.
    522535        IF ( dzsed_circont_abv(i) .GT. eps ) THEN
     536          !--If ice sedimentation is activated, the quantity of sedimentated ice was added
     537          !--to the total water vapor in the precipitation routine. Here we remove it
     538          !--(it will be reincluded later)
    523539          qised_circont_abv(i) = flsed_circont_abv(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
    524           qised_abv(i) = qised_abv(i) - qised_circont_abv(i)
    525           cfsed_abv(i) = cfsed_abv(i) - cfsed_circont_abv(i)
     540          qclr(i) = qclr(i) - qised_circont_abv(i)
    526541        ENDIF
    527542
    … …  
    537552        dqil_mix(i) = 0.
    538553        dqtl_mix(i) = 0.
     554        dcfl_sed(i) = 0.
     555        dqil_sed(i) = 0.
     556        dqtl_sed(i) = 0.
    539557        dcfc_sub(i) = 0.
    540558        dqic_sub(i) = 0.
    … …  
    543561        dqic_mix(i) = 0.
    544562        dqtc_mix(i) = 0.
     563        dcfc_sed(i) = 0.
     564        dqic_sed(i) = 0.
     565        dqtc_sed(i) = 0.
    545566      ELSE
    546567        lincontfra(i) = 0.
    … …  
    572593        cldfra(i) = MAX(0., cldfra(i) - lincontfra(i))
    573594        qcld(i) = MAX(0., qcld(i) - qlincont(i))
    574         qvc(i) = MAX(0., qvc(i) - qsat(i) * lincontfra(i))
     595        qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qsat(i) * lincontfra(i)))
    575596      ENDIF ! lincontfra(i) .GT. eps
    576597
    … …  
    592613        cldfra(i) = MAX(0., cldfra(i) - circontfra(i))
    593614        qcld(i) = MAX(0., qcld(i) - qcircont(i))
    594         qvc(i) = MAX(0., qvc(i) - qsat(i) * circontfra(i))
     615        qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qsat(i) * circontfra(i)))
    595616      ENDIF ! circontfra(i) .GT. eps
    596617
    … …  
    13771398
    13781399      IF ( ok_ice_sedim ) THEN
     1400        !--Initialisation
     1401        dcf_sed1  = 0.
     1402        dqvc_sed1 = 0.
     1403        dqi_sed1  = 0.
     1404        dcf_sed2  = 0.
     1405        dqvc_sed2 = 0.
     1406        dqi_sed2  = 0.
     1407        dcfl_sed1 = 0.
     1408        dqtl_sed1 = 0.
     1409        dqil_sed1 = 0.
     1410        dcfl_sed2 = 0.
     1411        dqtl_sed2 = 0.
     1412        dqil_sed2 = 0.
     1413        dcfc_sed1 = 0.
     1414        dqtc_sed1 = 0.
     1415        dqic_sed1 = 0.
     1416        dcfc_sed2 = 0.
     1417        dqtc_sed2 = 0.
     1418        dqic_sed2 = 0.
    13791419        !---------------------------------------
    13801420        !--         ICE SEDIMENTATION         --
    … …  
    13931433          qice_sedim = ( qcld(i) - qvc(i) ) * dzsed(i) / dz
    13941434          !--Tendencies
    1395           dcf_sed(i) = - cldfra(i) * dzsed(i) / dz
    1396           dqi_sed(i) = - qice_sedim
    1397           dqvc_sed(i) = - qvc(i) * dzsed(i) / dz
     1435          dcf_sed1 = - cldfra(i) * dzsed(i) / dz
     1436          dqi_sed1 = - qice_sedim
     1437          dqvc_sed1 = - qvc(i) * dzsed(i) / dz
    13981438          !--Convert to flux
    13991439          flsed(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
     1440          !--Save the vapor in the cloud (will be needed later)
     1441          qvapincld = qvc(i) / cldfra(i)
     1442          !--Add tendencies
     1443          cldfra(i)  = cldfra(i) + dcf_sed1
     1444          qcld(i)    = qcld(i) + dqvc_sed1 + dqi_sed1
     1445          qvc(i)     = qvc(i) + dqvc_sed1
     1446          clrfra(i)  = clrfra(i) - dcf_sed1
     1447          qclr(i)    = qclr(i) - dqvc_sed1 - dqi_sed1
     1448        ENDIF
     1449        !
     1450        IF ( lincontfra(i) .GT. eps ) THEN
     1451          icefall_velo = fallice_linear_contrails
     1452          dzsed_lincont(i) = icefall_velo * dtime * MAX(0., MIN(1., (2. - icefall_velo * dtime / dz)))
     1453          cfsed_lincont(i) = lincontfra(i)
     1454          qice_sedim = ( qlincont(i) - qsat(i) ) * dzsed_lincont(i) / dz
     1455          !--Tendencies
     1456          dcfl_sed1 = - lincontfra(i) * dzsed_lincont(i) / dz
     1457          dqil_sed1 = - qice_sedim
     1458          dqtl_sed1 = - qlincont(i) * dzsed_lincont(i) / dz
     1459          !--Convert to flux
     1460          flsed_lincont(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
     1461          !--Add tendencies
     1462          lincontfra(i)  = lincontfra(i) + dcfl_sed1
     1463          qlincont(i)    = qlincont(i) + dqtl_sed1
     1464          clrfra(i)  = clrfra(i) - dcfl_sed1
     1465          qclr(i)    = qclr(i) - dqtl_sed1
     1466        ENDIF
     1467        !
     1468        IF ( circontfra(i) .GT. eps ) THEN
     1469          icefall_velo = fallice_cirrus_contrails
     1470          dzsed_circont(i) = icefall_velo * dtime * MAX(0., MIN(1., (2. - icefall_velo * dtime / dz)))
     1471          cfsed_circont(i) = circontfra(i)
     1472          qice_sedim = ( qcircont(i) - qsat(i) ) * dzsed_circont(i) / dz
     1473          !--Tendencies
     1474          dcfc_sed1 = - circontfra(i) * dzsed_circont(i) / dz
     1475          dqic_sed1 = - qice_sedim
     1476          dqtc_sed1 = - qcircont(i) * dzsed_circont(i) / dz
     1477          !--Convert to flux
     1478          flsed_circont(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
     1479          !--Add tendencies
     1480          circontfra(i)  = circontfra(i) + dcfc_sed1
     1481          qcircont(i)    = qcircont(i) + dqtc_sed1
     1482          clrfra(i)  = clrfra(i) - dcfc_sed1
     1483          qclr(i)    = qclr(i) - dqtc_sed1
    14001484        ENDIF
    14011485        !
    … …  
    14171501          !--ice sedimentation
    14181502          !--(1) we use the clear-sky PDF to determine the humidity and increase the
    1419           !--cloud fraction / sublimate falling ice
     1503          !--cloud fraction / sublimate falling ice. NB it can also fall into
     1504          !--another cloud class
    14201505          !--(2) we do not increase cloud fraction and add the falling ice to the cloud
    14211506          !--(3) we increase the cloud fraction but use in-cloud water vapor as
    … …  
    14301515
    14311516          !--For region (1)
    1432           IF ( ( sedfra1 .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN
    1433 
    1434             IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    1435               tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
    1436               qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) )
    1437             ELSE
    1438               qvapinclr_lim = qsat(i) - qiceincld
     1517          IF ( sedfra1 .GT. eps ) THEN
     1518            IF ( clrfra(i) .GT. eps ) THEN
     1519              IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     1520                tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
     1521                qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) )
     1522              ELSE
     1523                qvapinclr_lim = qsat(i) - qiceincld
     1524              ENDIF
     1525
     1526              rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     1527              pdf_x = qvapinclr_lim / qsatl(i) * 100.
     1528              pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
     1529              pdf_x = qsat(i) / qsatl(i) * 100.
     1530              pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
     1531              IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
     1532                pdf_fra_above_lim = 0.
     1533                pdf_q_above_lim = 0.
     1534              ELSEIF ( pdf_y .LT. -10. ) THEN
     1535                pdf_fra_above_lim = clrfra(i)
     1536                pdf_q_above_lim = qclr(i)
     1537              ELSE
     1538                pdf_y = EXP( pdf_y )
     1539                pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     1540                pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
     1541                pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
     1542                pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
     1543                                  + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
     1544              ENDIF
     1545             
     1546              pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
     1547
     1548              !--The first three lines allow to favor ISSR rather than subsaturated sky for
     1549              !--the growth. The fourth line indicates that the ice is falling only in
     1550              !--the clear-sky area. The rest falls into other cloud classes
     1551              sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, &
     1552                      sedfra1 * chi_sedim * pdf_fra_above_lim &
     1553                      / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) &
     1554                      * clrfra(i) / (totfra_in(i) - cldfra(i))
     1555
     1556              IF ( pdf_fra_above_lim .GT. eps ) THEN
     1557                IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     1558                  dcf_sed2  = dcf_sed2  + sedfra_tmp
     1559                  dqvc_sed2 = dqvc_sed2 + sedfra_tmp * pdf_q_above_lim / pdf_fra_above_lim
     1560                  dqi_sed2  = dqi_sed2  + sedfra_tmp * qiceincld
     1561                ELSE
     1562                  dcf_sed2  = dcf_sed2  + sedfra_tmp
     1563                  dqvc_sed2 = dqvc_sed2 + sedfra_tmp * qsat(i)
     1564                  dqi_sed2  = dqi_sed2  + sedfra_tmp &
     1565                                   * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
     1566                ENDIF
     1567              ENDIF
     1568            ENDIF ! clrfra(i) .GT. eps
     1569
     1570            !--If the sedimented ice falls into a linear contrail, it increases its content
     1571            IF ( lincontfra(i) .GT. eps ) THEN
     1572                sedfra_tmp = sedfra1 * lincontfra(i) / (totfra_in(i) - cldfra(i))
     1573                dqil_sed2  = dqil_sed2 + sedfra_tmp * qiceincld
    14391574            ENDIF
    14401575
    1441             rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
    1442             pdf_x = qvapinclr_lim / qsatl(i) * 100.
    1443             pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
    1444             pdf_x = qsat(i) / qsatl(i) * 100.
    1445             pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
    1446             IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
    1447               pdf_fra_above_lim = 0.
    1448               pdf_q_above_lim = 0.
    1449             ELSEIF ( pdf_y .LT. -10. ) THEN
    1450               pdf_fra_above_lim = clrfra(i)
    1451               pdf_q_above_lim = qclr(i)
    1452             ELSE
    1453               pdf_y = EXP( pdf_y )
    1454               pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    1455               pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
    1456               pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
    1457               pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
    1458                                 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
     1576            !--If the sedimented ice falls into a contrail cirrus, it increases its content
     1577            IF ( circontfra(i) .GT. eps ) THEN
     1578                sedfra_tmp = sedfra1 * circontfra(i) / (totfra_in(i) - cldfra(i))
     1579                dqic_sed2  = dqic_sed2 + sedfra_tmp * qiceincld
    14591580            ENDIF
    1460            
    1461             sedfra1 = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, &
    1462                     sedfra1 * chi_sedim * pdf_fra_above_lim &
    1463                     / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim )))
    1464 
    1465             IF ( pdf_fra_above_lim .GT. eps ) THEN
    1466               IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    1467                 dcf_sed(i)  = dcf_sed(i)  + sedfra1
    1468                 dqvc_sed(i) = dqvc_sed(i) + sedfra1 * pdf_q_above_lim / pdf_fra_above_lim
    1469                 dqi_sed(i)  = dqi_sed(i)  + sedfra1 * qiceincld
    1470               ELSE
    1471                 dcf_sed(i)  = dcf_sed(i)  + sedfra1
    1472                 dqvc_sed(i) = dqvc_sed(i) + sedfra1 * qsat(i)
    1473                 dqi_sed(i)  = dqi_sed(i)  + sedfra1 &
    1474                                  * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
    1475               ENDIF
    1476             ENDIF
    1477           ENDIF ! ( sedfra1 .GT. eps .AND. ( clrfra(i) .GT. eps )
     1581          ENDIF ! sedfra1 .GT. eps
    14781582
    14791583          !--For region (2)
    1480           dqi_sed(i) = dqi_sed(i) + sedfra2 * qiceincld
     1584          dqi_sed2 = dqi_sed2 + sedfra2 * qiceincld
    14811585
    14821586          !--For region (3)
    14831587          IF ( cldfra(i) .GT. eps ) THEN
    1484             dcf_sed(i)  = dcf_sed(i)  + sedfra3
    1485             dqvc_sed(i) = dqvc_sed(i) + sedfra3 * qvc(i) / cldfra(i)
    1486             dqi_sed(i)  = dqi_sed(i)  + sedfra3 * qiceincld
     1588            dcf_sed2  = dcf_sed2  + sedfra3
     1589            dqvc_sed2 = dqvc_sed2 + sedfra3 * qvapincld
     1590            dqi_sed2  = dqi_sed2  + sedfra3 * qiceincld
    14871591          ENDIF
    1488         ENDIF ! qised_abv(i) .GT. 0.
     1592        ENDIF ! sedfra_abv .GT. eps
     1593
     1594        IF ( ok_plane_contrail ) THEN
     1595          sedfra_abv = MIN(dzsed_lincont_abv(i), dz) / dz * cfsed_lincont_abv(i)
     1596          IF ( sedfra_abv .GT. eps ) THEN
     1597
     1598            !--Three regions to be defined: (1) clear-sky, (2) cloudy, and
     1599            !--(3) region where it was previously cloudy but now clear-sky because of
     1600            !--ice sedimentation
     1601            !--(1) we use the clear-sky PDF to determine the humidity and increase the
     1602            !--cloud fraction / sublimate falling ice. NB it can also fall into
     1603            !--another cloud class
     1604            !--(2) we do not increase cloud fraction and add the falling ice to the cloud
     1605            !--(3) we increase the cloud fraction but use in-cloud water vapor as
     1606            !--background vapor
     1607            sedfra2 = MIN(cfsed_lincont(i), cfsed_lincont_abv(i)) &
     1608                    * MAX(0., MIN(dz, dzsed_lincont_abv(i)) - dzsed_lincont(i)) / dz
     1609            sedfra3 = MIN(cfsed_lincont(i), cfsed_lincont_abv(i)) &
     1610                    * MIN(MIN(dz, dzsed_lincont_abv(i)), dzsed_lincont(i)) / dz
     1611            sedfra1 = sedfra_abv - sedfra3 - sedfra2
     1612
     1613            qiceincld = qised_lincont_abv(i) / sedfra_abv
     1614
     1615            !--For region (1)
     1616            IF ( sedfra1 .GT. eps ) THEN
     1617              IF ( clrfra(i) .GT. eps ) THEN
     1618                qvapinclr_lim = qsat(i) - qiceincld
     1619
     1620                rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     1621                pdf_x = qvapinclr_lim / qsatl(i) * 100.
     1622                pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
     1623                pdf_x = qsat(i) / qsatl(i) * 100.
     1624                pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
     1625                IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
     1626                  pdf_fra_above_lim = 0.
     1627                  pdf_q_above_lim = 0.
     1628                ELSEIF ( pdf_y .LT. -10. ) THEN
     1629                  pdf_fra_above_lim = clrfra(i)
     1630                  pdf_q_above_lim = qclr(i)
     1631                ELSE
     1632                  pdf_y = EXP( pdf_y )
     1633                  pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     1634                  pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
     1635                  pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
     1636                  pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
     1637                                    + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
     1638                ENDIF
     1639               
     1640                pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
     1641
     1642                !--The first three lines allow to favor ISSR rather than subsaturated sky for
     1643                !--the growth. The fourth line indicates that the ice is falling only in
     1644                !--the clear-sky area. The rest falls into other cloud classes
     1645                sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, &
     1646                        sedfra1 * chi_sedim * pdf_fra_above_lim &
     1647                        / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) &
     1648                        * clrfra(i) / (totfra_in(i) - lincontfra(i))
     1649
     1650                IF ( pdf_fra_above_lim .GT. eps ) THEN
     1651                  dcfl_sed2 = dcfl_sed2 + sedfra_tmp
     1652                  dqtl_sed2 = dqtl_sed2 + sedfra_tmp * qsat(i)
     1653                  dqil_sed2 = dqil_sed2 + sedfra_tmp &
     1654                                   * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
     1655                ENDIF
     1656              ENDIF ! clrfra(i) .GT. eps
     1657
     1658              !--If the sedimented ice falls into a natural cirrus, it increases its content
     1659              IF ( cldfra(i) .GT. eps ) THEN
     1660                  sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - lincontfra(i))
     1661                  dqi_sed2  = dqi_sed2 + sedfra_tmp * qiceincld
     1662              ENDIF
     1663
     1664              !--If the sedimented ice falls into a contrail cirrus, it increases its content
     1665              IF ( circontfra(i) .GT. eps ) THEN
     1666                  sedfra_tmp = sedfra1 * circontfra(i) / (totfra_in(i) - lincontfra(i))
     1667                  dqic_sed2  = dqic_sed2 + sedfra_tmp * qiceincld
     1668              ENDIF
     1669            ENDIF ! sedfra1 .GT. eps
     1670
     1671            !--For region (2)
     1672            dqil_sed2 = dqil_sed2 + sedfra2 * qiceincld
     1673
     1674            !--For region (3)
     1675            IF ( lincontfra(i) .GT. eps ) THEN
     1676              dcfl_sed2 = dcfl_sed2 + sedfra3
     1677              dqtl_sed2 = dqtl_sed2 + sedfra3 * (qsat(i) + qiceincld)
     1678              dqil_sed2 = dqil_sed2 + sedfra3 * qiceincld
     1679            ENDIF
     1680          ENDIF ! sedfra_abv .GT. eps
     1681
     1682          !
     1683          sedfra_abv = MIN(dzsed_circont_abv(i), dz) / dz * cfsed_circont_abv(i)
     1684          IF ( sedfra_abv .GT. eps ) THEN
     1685
     1686            !--Three regions to be defined: (1) clear-sky, (2) cloudy, and
     1687            !--(3) region where it was previously cloudy but now clear-sky because of
     1688            !--ice sedimentation
     1689            !--(1) we use the clear-sky PDF to determine the humidity and increase the
     1690            !--cloud fraction / sublimate falling ice. NB it can also fall into
     1691            !--another cloud class
     1692            !--(2) we do not increase cloud fraction and add the falling ice to the cloud
     1693            !--(3) we increase the cloud fraction but use in-cloud water vapor as
     1694            !--background vapor
     1695            sedfra2 = MIN(cfsed_circont(i), cfsed_circont_abv(i)) &
     1696                    * MAX(0., MIN(dz, dzsed_circont_abv(i)) - dzsed_circont(i)) / dz
     1697            sedfra3 = MIN(cfsed_circont(i), cfsed_circont_abv(i)) &
     1698                    * MIN(MIN(dz, dzsed_circont_abv(i)), dzsed_circont(i)) / dz
     1699            sedfra1 = sedfra_abv - sedfra3 - sedfra2
     1700
     1701            qiceincld = qised_circont_abv(i) / sedfra_abv
     1702
     1703            !--For region (1)
     1704            IF ( sedfra1 .GT. eps ) THEN
     1705              IF ( clrfra(i) .GT. eps ) THEN
     1706                qvapinclr_lim = qsat(i) - qiceincld
     1707
     1708                rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     1709                pdf_x = qvapinclr_lim / qsatl(i) * 100.
     1710                pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
     1711                pdf_x = qsat(i) / qsatl(i) * 100.
     1712                pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
     1713                IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
     1714                  pdf_fra_above_lim = 0.
     1715                  pdf_q_above_lim = 0.
     1716                ELSEIF ( pdf_y .LT. -10. ) THEN
     1717                  pdf_fra_above_lim = clrfra(i)
     1718                  pdf_q_above_lim = qclr(i)
     1719                ELSE
     1720                  pdf_y = EXP( pdf_y )
     1721                  pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     1722                  pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
     1723                  pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
     1724                  pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
     1725                                    + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
     1726                ENDIF
     1727               
     1728                pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim
     1729
     1730                !--The first three lines allow to favor ISSR rather than subsaturated sky for
     1731                !--the growth. The fourth line indicates that the ice is falling only in
     1732                !--the clear-sky area. The rest falls into other cloud classes
     1733                sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, &
     1734                        sedfra1 * chi_sedim * pdf_fra_above_lim &
     1735                        / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) &
     1736                        * clrfra(i) / (totfra_in(i) - circontfra(i))
     1737
     1738                IF ( pdf_fra_above_lim .GT. eps ) THEN
     1739                  dcfc_sed2 = dcfc_sed2 + sedfra_tmp
     1740                  dqtc_sed2 = dqtc_sed2 + sedfra_tmp * qsat(i)
     1741                  dqic_sed2 = dqic_sed2 + sedfra_tmp &
     1742                                   * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i))
     1743                ENDIF
     1744              ENDIF ! clrfra(i) .GT. eps
     1745
     1746              !--If the sedimented ice falls into a natural cirrus, it increases its content
     1747              IF ( cldfra(i) .GT. eps ) THEN
     1748                  sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - circontfra(i))
     1749                  dqi_sed2  = dqi_sed2 + sedfra_tmp * qiceincld
     1750              ENDIF
     1751
     1752              !--If the sedimented ice falls into a contrail cirrus, it increases its content
     1753              IF ( lincontfra(i) .GT. eps ) THEN
     1754                  sedfra_tmp = sedfra1 * lincontfra(i) / (totfra_in(i) - circontfra(i))
     1755                  dqil_sed2  = dqil_sed2 + sedfra_tmp * qiceincld
     1756              ENDIF
     1757            ENDIF ! sedfra1 .GT. eps
     1758
     1759            !--For region (2)
     1760            dqic_sed2 = dqic_sed2 + sedfra2 * qiceincld
     1761
     1762            !--For region (3)
     1763            IF ( circontfra(i) .GT. eps ) THEN
     1764              dcfc_sed2 = dcfc_sed2 + sedfra3
     1765              dqtc_sed2 = dqtc_sed2 + sedfra3 * (qsat(i) + qiceincld)
     1766              dqic_sed2 = dqic_sed2 + sedfra3 * qiceincld
     1767            ENDIF
     1768          ENDIF ! sedfra_abv .GT. eps
     1769        ENDIF ! ok_plane_contrails
    14891770
    14901771        !--Add tendencies
    1491         cldfra(i)  = MIN(totfra_in(i), cldfra(i) + dcf_sed(i))
    1492         qcld(i)    = qcld(i) + dqvc_sed(i) + dqi_sed(i)
    1493         qvc(i)     = qvc(i) + dqvc_sed(i)
    1494         clrfra(i)  = MAX(0., clrfra(i) - dcf_sed(i))
     1772        cldfra(i)   = cldfra(i) + dcf_sed2
     1773        qcld(i)     = qcld(i)   + dqvc_sed2 + dqi_sed2
     1774        qvc(i)      = qvc(i)    + dqvc_sed2
     1775        clrfra(i)   = clrfra(i) - dcf_sed2
     1776        qclr(i)     = qclr(i)   - dqvc_sed2 - dqi_sed2
    14951777        !--We re-include sublimated sedimentated ice into clear sky water vapor
    1496         qclr(i)    = qclr(i) - dqvc_sed(i) + qised_abv(i) - dqi_sed(i)
     1778        qclr(i)     = qclr(i)   + qised_abv(i)
     1779        !--Diagnostics
     1780        dcf_sed(i)  = dcf_sed1  + dcf_sed2
     1781        dqvc_sed(i) = dqvc_sed1 + dqvc_sed2
     1782        dqi_sed(i)  = dqi_sed1  + dqi_sed2
     1783
     1784        if ( ok_plane_contrail ) THEN
     1785          !--Add tendencies
     1786          lincontfra(i) = lincontfra(i) + dcfl_sed2
     1787          qlincont(i)   = qlincont(i)   + dqtl_sed2
     1788          clrfra(i)     = clrfra(i) - dcfl_sed2
     1789          qclr(i)       = qclr(i)   - dqtl_sed2
     1790          !--We re-include sublimated sedimentated ice into clear sky water vapor
     1791          qclr(i)       = qclr(i) + qised_lincont_abv(i)
     1792          !--Diagnostics
     1793          dcfl_sed(i) = dcfl_sed1 + dcfl_sed2
     1794          dqtl_sed(i) = dqtl_sed1 + dqtl_sed2
     1795          dqil_sed(i) = dqil_sed1 + dqil_sed2
     1796
     1797          !--Add tendencies
     1798          circontfra(i) = circontfra(i) + dcfc_sed2
     1799          qcircont(i)   = qcircont(i)   + dqtc_sed2
     1800          clrfra(i)     = clrfra(i) - dcfc_sed2
     1801          qclr(i)       = qclr(i)   - dqtc_sed2
     1802          !--We re-include sublimated sedimentated ice into clear sky water vapor
     1803          qclr(i)       = qclr(i) + qised_circont_abv(i)
     1804          !--Diagnostics
     1805          dcfc_sed(i) = dcfc_sed1 + dcfc_sed2
     1806          dqtc_sed(i) = dqtc_sed1 + dqtc_sed2
     1807          dqic_sed(i) = dqic_sed1 + dqic_sed2
     1808        ENDIF
    14971809
    14981810      ENDIF ! ok_ice_sedim

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

    r5717 r5728  
    137137USE phys_local_var_mod, ONLY : dcfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix
    138138USE phys_local_var_mod, ONLY : dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix
     139USE phys_local_var_mod, ONLY : dcfl_sed, dqil_sed, dqtl_sed, dcfc_sed, dqic_sed, dqtc_sed
    139140USE geometry_mod, ONLY: longitude_deg, latitude_deg
    140141
    … …  
    464465shear(:)        = 0.
    465466flsed(:)        = 0.
     467flsed_lincont(:)= 0.
     468flsed_circont(:)= 0.
    466469pt_pron_clds(:) = .FALSE.
    467470
    … …  
    570573      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
    571574                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
    572                         zqvapclr, zqupnew, flsed, &
     575                        zqvapclr, zqupnew, flsed, flsed_lincont, flsed_circont, &
    573576                        cldfra_in, qvc_in, qliq_in, qice_in, &
    574577                        zrfl, zrflclr, zrflcld, &
    … …  
    581584
    582585      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
    583                         zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, flsed, &
     586                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
     587                        flsed, flsed_lincont, flsed_circont, &
    584588                        zrfl, zrflclr, zrflcld, &
    585589                        zifl, ziflclr, ziflcld, &
    … …  
    939943                        dcfl_cir(:,k), dqtl_cir(:,k), &
    940944                        dcfl_mix(:,k), dqil_mix(:,k), dqtl_mix(:,k), &
     945                        dcfl_sed(:,k), dqil_sed(:,k), dqtl_sed(:,k), &
     946                        dcfc_sed(:,k), dqic_sed(:,k), dqtc_sed(:,k), &
    941947                        dcfc_sub(:,k), dqic_sub(:,k), dqtc_sub(:,k), &
    942948                        dcfc_mix(:,k), dqic_mix(:,k), dqtc_mix(:,k))
    … …  
    10991105
    11001106                        IF ( ok_ice_sedim ) THEN
    1101                           qice_sedim = flsed(i) / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
     1107                          qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) &
     1108                              / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
    11021109                          ! Add the ice that was sedimented, as it is not included in zqn
    11031110                          qlibef(i) = qlibef(i) + qice_sedim
    … …  
    11791186
    11801187                IF ( ok_ice_sedim ) THEN
    1181                   qice_sedim = flsed(i) / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
     1188                  qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) &
     1189                      / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
    11821190                  ! Remove the ice that was sedimented. As it is not included in zqn,
    11831191                  ! we only remove it from the total water
    … …  
    12481256                            ctot_vol(:,k), ptconv(:,k), &
    12491257                            zt, zq, zoliql, zoliqi, zfice, &
    1250                             rneb(:,k), flsed, znebprecipclr, znebprecipcld, &
     1258                            rneb(:,k), znebprecipclr, znebprecipcld, &
    12511259                            zrfl, zrflclr, zrflcld, &
    12521260                            zifl, ziflclr, ziflcld, &
    … …  
    12651273      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
    12661274                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
    1267                             zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, flsed, &
     1275                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
    12681276                            rneb(:,k), znebprecipclr, znebprecipcld, &
    12691277                            zneb, tot_zneb, zrho_up, zvelo_up, &
    … …  
    14291437        !--the sink of condensed water from precipitation
    14301438        IF ( ptconv(i,k) ) THEN
    1431           IF ( zoliq(i) .GT. 0. ) THEN
     1439          IF ( zcond(i) .GT. 0. ) THEN
    14321440            qvcon_old(i,k) = qvcon(i,k)
    14331441            qccon_old(i,k) = qccon(i,k) * zoliq(i) / zcond(i)
    … …  
    15721580  IF ( ok_ice_sedim ) THEN
    15731581    DO i = 1, klon
    1574       snow(i) = snow(i) + flsed(i)
     1582      snow(i) = snow(i) + flsed(i) + flsed_lincont(i) + flsed_circont(i)
    15751583    ENDDO
    15761584  ENDIF

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

    r5717 r5728  
    1818SUBROUTINE histprecip_precld( &
    1919           klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, &
    20            zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
     20           zmqc, zneb, znebprecip, znebprecipclr, flsed, flsed_lincont, flsed_circont, &
    2121           zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub &
    2222           )
    … …  
    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]
    46 REAL,    INTENT(IN),    DIMENSION(klon) :: icesed_flux    !--sedimentated ice flux [kg/s/m2]
     46REAL,    INTENT(IN),    DIMENSION(klon) :: flsed          !--sedimentated ice flux [kg/s/m2]
     47REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_lincont  !--linear contrails sedimentated ice flux [kg/s/m2]
     48REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_circont  !--contrail cirrus sedimentated ice flux [kg/s/m2]
    4749
    4850REAL,    INTENT(IN),    DIMENSION(klon) :: zneb           !--cloud fraction IN THE LAYER ABOVE [-]
    … …  
    132134    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
    133135    !--added to the total water content of the gridbox
    134     IF ( icesed_flux(i) .GT. 0. ) THEN
    135       qice_sedim = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
     136    IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
     137      qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) &
     138          / ( paprsdn(i) - paprsup(i) ) * RG * dtime
    136139
    137140      !--Thermalisation of sedimentated ice
    … …  
    339342SUBROUTINE histprecip_postcld( &
    340343        klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, pt_pron_clds, &
    341         zdqsdT_raw, zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
     344        zdqsdT_raw, zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
    342345        rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, &
    343346        zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
    … …  
    394397REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
    395398REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
    396 REAL,    INTENT(INOUT), DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
    397399
    398400REAL,    INTENT(OUT),   DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
    … …  
    438440zqpreci(:) = 0.
    439441ziflprev(:) = 0.
    440 
    441 ! AB ICE SEDIM
    442 !icesed_flux(:) = 0.
    443442
    444443
    … …  
    700699
    701700ENDDO
    702 
    703 !---------------------------------------------------------
    704 !--                  ICE SEDIMENTATION
    705 !---------------------------------------------------------
    706 !IF ( ok_ice_sedim ) THEN
    707 IF ( .FALSE. ) THEN ! AB ICE SEDIM
    708 
    709   DO i = 1, klon
    710     IF ( ( zoliqi(i) .GT. 0. ) .AND. ( rneb(i) .GT. eps ) .AND. pt_pron_clds(i) ) THEN
    711 
    712       iwcg = MAX( zrho(i) * zoliqi(i) / zneb(i) * 1000., eps ) ! iwc in g/m3
    713       !tempc = zt(i) - RTT ! celcius temp
    714       !! so-called 'empirical parameterization' in Stubenrauch et al. 2019
    715       !IF ( tempc .GE. -60. ) THEN
    716       !  icevelo = EXP( 1.9714 + 0.00216078 * tempc &
    717       !          - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) )
    718       !ELSE
    719       !  icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15
    720       !ENDIF
    721       !icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
    722 
    723       icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16
    724 
    725       qice_sedim = zoliqi(i) * MIN(1., icevelo * dtime / zdz(i))
    726       !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two
    727       !--times in the same timestep)
    728       qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)
    729       !--Add tendencies
    730       zoliqi(i) = zoliqi(i) - qice_sedim
    731       zoliq(i)  = zoliq(i)  - qice_sedim
    732       !--Convert to flux
    733       icesed_flux(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
    734 
    735       !--Diagnostic tendencies
    736       dqised(i) = dqised(i) - qice_sedim / dtime
    737     ENDIF
    738   ENDDO
    739 ENDIF
    740701
    741702! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
    … …  
    796757SUBROUTINE poprecip_precld( &
    797758           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
    798            qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, icesed_flux, &
     759           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
     760           flsed, flsed_lincont, flsed_circont, &
    799761           cldfra, qvc, qliq, qice, &
    800762           rain, rainclr, raincld, snow, snowclr, snowcld, &
    … …  
    831793REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
    832794REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
    833 REAL,    INTENT(IN),    DIMENSION(klon) :: icesed_flux    !--sedimentated ice water flux [kg/s/m2]
     795REAL,    INTENT(IN),    DIMENSION(klon) :: flsed          !--sedimentated ice water flux [kg/s/m2]
     796REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_lincont  !--sedimentated ice water flux [kg/s/m2]
     797REAL,    INTENT(IN),    DIMENSION(klon) :: flsed_circont  !--sedimentated ice water flux [kg/s/m2]
    834798
    835799REAL,    INTENT(INOUT), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
    … …  
    933897  cpw = RCPD * RVTMP2
    934898  DO i = 1, klon
    935     IF ( icesed_flux(i) .GT. 0. ) THEN
    936       qice_sedim = icesed_flux(i) / dhum_to_dflux(i)
     899    IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
     900      qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) / dhum_to_dflux(i)
    937901      !--No condensed water so cp=cp(vapor+dry air)
    938902      !-- RVTMP2=rcpv/rcpd-1
    … …  
    10421006    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
    10431007    !--added to the total water content of the gridbox
    1044     IF ( icesed_flux(i) .GT. 0. ) THEN
    1045       qice_sedim = icesed_flux(i) / dhum_to_dflux(i)
     1008    IF ( (flsed(i) + flsed_lincont(i) + flsed_circont(i)) .GT. 0. ) THEN
     1009      qice_sedim = (flsed(i) + flsed_lincont(i) + flsed_circont(i)) / dhum_to_dflux(i)
    10461010      !--Vapor is updated after evaporation/sublimation (it is increased)
    10471011      qvap(i) = qvap(i) + qice_sedim
    … …  
    14081372SUBROUTINE poprecip_postcld( &
    14091373           klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, &
    1410            temp, qvap, qliq, qice, icefrac, cldfra, icesed_flux, &
     1374           temp, qvap, qliq, qice, icefrac, cldfra, &
    14111375           precipfracclr, precipfraccld, &
    14121376           rain, rainclr, raincld, snow, snowclr, snowcld, &
    … …  
    14641428REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
    14651429
    1466 REAL,    INTENT(INOUT), DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
    14671430REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
    14681431REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
    … …  
    15281491
    15291492qzero(:)    = 0.
    1530 ! AB ICE SEDIM
    1531 !icesed_flux(:) = 0.
    15321493
    15331494dqrcol(:)   = 0.
    … …  
    20001961
    20011962
    2002   !---------------------------------------------------------
    2003   !--                  ICE SEDIMENTATION
    2004   !---------------------------------------------------------
    2005   !IF ( ok_ice_sedim ) THEN
    2006   IF ( .FALSE. ) THEN ! AB ICE SEDIM
    2007 
    2008     IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. eps ) ) THEN
    2009 
    2010       rho = pplay(i) / temp(i) / RD
    2011       iwcg = MAX( rho * qice(i) / cldfra(i) * 1000., eps ) ! iwc in g/m3
    2012       !tempc = temp(i) - RTT ! celcius temp
    2013       !! so-called 'empirical parameterization' in Stubenrauch et al. 2019
    2014       !IF ( tempc .GE. -60. ) THEN
    2015       !  icevelo = EXP( 1.9714 + 0.00216078 * tempc &
    2016       !          - ( 0.0000414122 * tempc**2 + 0.00538922 * tempc + 0.0516344 ) * LOG(iwcg) )
    2017       !ELSE
    2018       !  icevelo = 65. * iwcg**0.2 * ( 15000. / pplay(i) )**0.15
    2019       !ENDIF
    2020       !icevelo = fallice_sedim * icevelo / 100.0 ! from cm/s to m/s
    2021 
    2022       icevelo = fallice_sedim * 1.645 * ( iwcg / 1000. )**0.16
    2023 
    2024       dz = ( paprsdn(i) - paprsup(i) ) / RG / rho
    2025       qice_sedim = qice(i) * MIN(1., icevelo * dtime / dz)
    2026       !--We remove the ice that was sedimentated in the gridbox (it cannot sedimentate two
    2027       !--times in the same timestep)
    2028       qice_sedim = MAX(0., qice_sedim - dqised(i) * dtime)
    2029       !--Add tendencies
    2030       qice(i) = qice(i) - qice_sedim
    2031       !--Convert to flux
    2032       icesed_flux(i) = qice_sedim * dhum_to_dflux(i)
    2033 
    2034       !--Diagnostic tendencies
    2035       dqised(i) = dqised(i) - qice_sedim / dtime
    2036     ENDIF
    2037   ENDIF
    2038 
    2039 
    20401963  !--If the local flux of rain+snow in clear air is lower than min_precip,
    20411964  !--we reduce the precipiration fraction in the clear air so that the new

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

    r5717 r5728  
    707707      REAL, SAVE, ALLOCATABLE :: dcfl_mix(:,:), dqil_mix(:,:), dqtl_mix(:,:)
    708708      !$OMP THREADPRIVATE(dcfl_mix, dqil_mix, dqtl_mix)
     709      REAL, SAVE, ALLOCATABLE :: dcfl_sed(:,:), dqil_sed(:,:), dqtl_sed(:,:)
     710      !$OMP THREADPRIVATE(dcfl_sed, dqil_sed, dqtl_sed)
     711      REAL, SAVE, ALLOCATABLE :: dcfc_sed(:,:), dqic_sed(:,:), dqtc_sed(:,:)
     712      !$OMP THREADPRIVATE(dcfc_sed, dqic_sed, dqtc_sed)
    709713      REAL, SAVE, ALLOCATABLE :: dcfc_sub(:,:), dqic_sub(:,:), dqtc_sub(:,:)
    710714      !$OMP THREADPRIVATE(dcfc_sub, dqic_sub, dqtc_sub)
    … …  
    13041308      ALLOCATE(dcfl_sub(klon,klev), dqil_sub(klon,klev), dqtl_sub(klon,klev))
    13051309      ALLOCATE(dcfl_mix(klon,klev), dqil_mix(klon,klev), dqtl_mix(klon,klev))
     1310      ALLOCATE(dcfl_sed(klon,klev), dqil_sed(klon,klev), dqtl_sed(klon,klev))
    13061311      ALLOCATE(dcfc_sub(klon,klev), dqic_sub(klon,klev), dqtc_sub(klon,klev))
    13071312      ALLOCATE(dcfc_mix(klon,klev), dqic_mix(klon,klev), dqtc_mix(klon,klev))
     1313      ALLOCATE(dcfc_sed(klon,klev), dqic_sed(klon,klev), dqtc_sed(klon,klev))
    13081314      ALLOCATE(cldfra_nocont(klon,klev), cldtau_nocont(klon,klev), cldemi_nocont(klon,klev))
    13091315      ALLOCATE(cldh_nocont(klon), contcov(klon), conttau(klon,klev), contemi(klon,klev))
    … …  
    17361742      DEALLOCATE(dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, dcfl_mix, dqil_mix, dqtl_mix)
    17371743      DEALLOCATE(dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix)
     1744      DEALLOCATE(dcfl_sed, dqil_sed, dqtl_sed, dcfc_sed, dqic_sed, dqtc_sed)
    17381745      DEALLOCATE(cldfra_nocont, cldtau_nocont, cldemi_nocont, conttau, contemi)
    17391746      DEALLOCATE(cldh_nocont, contcov, fiwp_nocont, fiwc_nocont, ref_ice_nocont)

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

    r5717 r5728  
    25592559  TYPE(ctrl_out), SAVE :: o_dqtlmix = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    25602560    'dqtlmix', 'Mixing linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2561  TYPE(ctrl_out), SAVE :: o_dcflsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2562    'dcflsed', 'Ice sedimentation linear contrail fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2563  TYPE(ctrl_out), SAVE :: o_dqilsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2564    'dqilsed', 'Ice sedimentation linear contrail ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2565  TYPE(ctrl_out), SAVE :: o_dqtlsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2566    'dqtlsed', 'Ice sedimentation linear contrail total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2567  TYPE(ctrl_out), SAVE :: o_dcfcsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2568    'dcfcsed', 'Ice sedimentation contrail cirrus fraction tendency', 's-1', (/ ('', i=1, 10) /))
     2569  TYPE(ctrl_out), SAVE :: o_dqicsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2570    'dqicsed', 'Ice sedimentation contrail cirrus ice specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
     2571  TYPE(ctrl_out), SAVE :: o_dqtcsed = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     2572    'dqtcsed', 'Ice sedimentation contrail cirrus total specific humidity tendency', 'kg/kg/s', (/ ('', i=1, 10) /))
    25612573  TYPE(ctrl_out), SAVE :: o_dcfcsub = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    25622574    'dcfcsub', 'Sublimation contrail cirrus fraction tendency', 's-1', (/ ('', i=1, 10) /))

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

    r5717 r5728  
    245245         o_dcflini, o_dqilini, o_dqtlini, o_dcflsub, o_dqilsub, o_dqtlsub, &
    246246         o_dcflmix, o_dqilmix, o_dqtlmix, &
     247         o_dcflsed, o_dqilsed, o_dqtlsed, o_dcfcsed, o_dqicsed, o_dqtcsed, &
    247248         o_dcfcsub, o_dqicsub, o_dqtcsub, o_dcfcmix, o_dqicmix, o_dqtcmix, &
    248249         o_cldfra_nocont, o_cldtau_nocont, o_cldemi_nocont, o_cldh_nocont, &
    … …  
    382383         dcfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, &
    383384         dcfl_mix, dqil_mix, dqtl_mix, &
     385         dcfl_sed, dqil_sed, dqtl_sed, dcfc_sed, dqic_sed, dqtc_sed, &
    384386         dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix, &
    385387         cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
    … …  
    22562258         CALL histwrite_phy(o_dqilmix, dqil_mix)
    22572259         CALL histwrite_phy(o_dqtlmix, dqtl_mix)
     2260         CALL histwrite_phy(o_dcflsed, dcfl_sed)
     2261         CALL histwrite_phy(o_dqilsed, dqil_sed)
     2262         CALL histwrite_phy(o_dqtlsed, dqtl_sed)
    22582263         CALL histwrite_phy(o_dcfcsub, dcfc_sub)
    22592264         CALL histwrite_phy(o_dqicsub, dqic_sub)
    … …  
    22622267         CALL histwrite_phy(o_dqicmix, dqic_mix)
    22632268         CALL histwrite_phy(o_dqtcmix, dqtc_mix)
     2269         CALL histwrite_phy(o_dcfcsed, dcfc_sed)
     2270         CALL histwrite_phy(o_dqicsed, dqic_sed)
     2271         CALL histwrite_phy(o_dqtcsed, dqtc_sed)
    22642272         CALL histwrite_phy(o_cldfra_nocont, cldfra_nocont)
    22652273         CALL histwrite_phy(o_cldtau_nocont, cldtau_nocont)
Note: See TracChangeset for help on using the changeset viewer.