Changeset 5691


Ignore:
Timestamp:
Jun 6, 2025, 11:55:21 AM (3 weeks ago)
Author:
aborella
Message:

Adaptations made to contrails radiative transfer and to ice sedimentation, following discussion with Bernd Karcher

Location:
LMDZ6/branches/contrails/libf/phylmd
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r5643 r5691  
    3333  USE lmdz_cloud_optics_prop_ini , ONLY : rei_coef, rei_min_temp
    3434  USE lmdz_cloud_optics_prop_ini , ONLY : zepsec, novlp, iflag_ice_thermo, ok_new_lscp
    35   USE lmdz_cloud_optics_prop_ini , ONLY : ok_plane_contrail, re_ice_crystals_contrails
     35  USE lmdz_cloud_optics_prop_ini , ONLY : ok_plane_contrail
     36  USE lmdz_cloud_optics_prop_ini , ONLY : rei_linear_contrails, rei_cirrus_contrails
    3637 
    3738
    … …  
    691692        IF ( lincontfra(i,k) .GT. zepsec ) THEN
    692693          pclc_lincont(i,k) = lincontfra(i,k)
    693           rei_lincont = re_ice_crystals_contrails * 1.E6
     694          rei_lincont = rei_linear_contrails * 1.E6
    694695        ELSE
    695696          pclc_lincont(i,k) = 0.
    … …  
    702703          pclc_circont(i,k) = circontfra(i,k)
    703704
    704           tc = temp(i, k) - 273.15
    705           rei_circont = d_rei_dt*tc + rei_max
    706           IF (tc<=-81.4) rei_circont = rei_min
     705          !tc = temp(i, k) - 273.15
     706          !rei_circont = d_rei_dt*tc + rei_max
     707          !IF (tc<=-81.4) rei_circont = rei_min
     708          rei_circont = rei_cirrus_contrails * 1.E6
    707709
    708710        ELSE

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

    r5639 r5691  
    2424  REAL, PROTECTED :: rei_coef, rei_min_temp
    2525  REAL, PROTECTED :: zepsec
    26   REAL, PROTECTED :: re_ice_crystals_contrails=7.5E-6 ! [m] effective radius of ice crystals in contrails
     26  REAL, PROTECTED :: rei_linear_contrails=7.5E-6 ! [m] effective radius of ice crystals in linear contrails
     27  REAL, PROTECTED :: rei_cirrus_contrails=10.E-6 ! [m] effective radius of ice crystals in contrails cirrus
    2728  REAL, PARAMETER :: thres_tau=0.3, thres_neb=0.001
    2829  REAL, PARAMETER :: prmhc=440.*100. ! Pressure between medium and high level cloud in Pa
    … …  
    4344!$OMP THREADPRIVATE(rei_coef, rei_min_temp)
    4445!$OMP THREADPRIVATE(zepsec)
    45 !$OMP THREADPRIVATE(re_ice_crystals_contrails, seuil_neb)
     46!$OMP THREADPRIVATE(rei_linear_contrails, rei_cirrus_contrails, seuil_neb)
    4647
    4748 
    … …  
    110111    rei_min_temp = 175.
    111112    CALL getin_p('rei_min_temp',rei_min_temp)
    112     CALL getin_p('re_ice_crystals_contrails', re_ice_crystals_contrails)
     113    CALL getin_p('rei_linear_contrails', rei_linear_contrails)
     114    write(lunout,*)'rei_linear_contrails=',rei_linear_contrails
     115    CALL getin_p('rei_cirrus_contrails', rei_cirrus_contrails)
     116    write(lunout,*)'rei_cirrus_contrails=',rei_cirrus_contrails
    113117    CALL getin_p('seuil_neb_rad', seuil_neb)
    114118    write(lunout,*)'seuil_neb_rad=',seuil_neb

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

    r5684 r5691  
    124124USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
    125125USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato
    126 USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_precip_lincontrails, ok_ice_sedim
     126USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_precip_contrails, ok_ice_sedim
    127127USE lmdz_lscp_ini, ONLY : ok_nodeep_lscp, ok_nodeep_lscp_rad
    128128
    … …  
    322322  REAL, DIMENSION(klon) :: ztupnew
    323323  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
    324   REAL, DIMENSION(klon) :: cldfra_above, icesed_flux ! for sedimentation of ice crystals
    325324  REAL, DIMENSION(klon) :: zrflclr, zrflcld
    326325  REAL, DIMENSION(klon) :: ziflclr, ziflcld
    … …  
    330329  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
    331330  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
     331  ! for ice sedimentation
     332  REAL, DIMENSION(klon) :: dzsed, flsed, cfsed
     333  REAL, DIMENSION(klon) :: dzsed_abv, flsed_abv, cfsed_abv
     334  REAL :: qice_sedim
    332335
    333336  ! for quantity of condensates seen by radiation
    … …  
    342345  REAL, DIMENSION(klon) :: totfra_in, qtot_in
    343346  LOGICAL, DIMENSION(klon) :: pt_pron_clds
     347  REAL, DIMENSION(klon) :: dzsed_lincont, flsed_lincont, cfsed_lincont
     348  REAL, DIMENSION(klon) :: dzsed_circont, flsed_circont, cfsed_circont
     349  REAL, DIMENSION(klon) :: dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv
     350  REAL, DIMENSION(klon) :: dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv
    344351  REAL :: qice_cont
    345352  !--for Lamquin et al 2012 diagnostics
    … …  
    448455qvc(:)          = 0.
    449456shear(:)        = 0.
     457flsed(:)        = 0.
    450458pt_pron_clds(:) = .FALSE.
    451459
    … …  
    480488zqupnew(:)    = 0.
    481489zqvapclr(:)   = 0.
    482 cldfra_above(:)= 0.
    483 icesed_flux(:)= 0.
    484490
    485491
    … …  
    537543      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    538544                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
    539                         zqvapclr, zqupnew, icesed_flux, &
     545                        zqvapclr, zqupnew, flsed, &
    540546                        cldfra_in, qvc_in, qliq_in, qice_in, &
    541547                        zrfl, zrflclr, zrflcld, &
    … …  
    548554
    549555      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    550                         zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, icesed_flux, &
     556                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, flsed, &
    551557                        zrfl, zrflclr, zrflcld, &
    552558                        zifl, ziflclr, ziflcld, &
    … …  
    670676          !--Initialisation
    671677          IF ( ok_plane_contrail ) THEN
    672             lincontfra(:) = 0.
    673             circontfra(:) = 0.
    674             qlincont(:)   = 0.
    675             qcircont(:)   = 0.
     678            IF ( iftop ) THEN
     679              dzsed_lincont_abv(:) = 0.
     680              flsed_lincont_abv(:) = 0.
     681              cfsed_lincont_abv(:) = 0.
     682              dzsed_circont_abv(:) = 0.
     683              flsed_circont_abv(:) = 0.
     684              cfsed_circont_abv(:) = 0.
     685            ELSE
     686              dzsed_lincont_abv(:) = dzsed_lincont(:)
     687              flsed_lincont_abv(:) = flsed_lincont(:)
     688              cfsed_lincont_abv(:) = cfsed_lincont(:)
     689              dzsed_circont_abv(:) = dzsed_circont(:)
     690              flsed_circont_abv(:) = flsed_circont(:)
     691              cfsed_circont_abv(:) = cfsed_circont(:)
     692            ENDIF
     693            dzsed_lincont(:) = 0.
     694            flsed_lincont(:) = 0.
     695            cfsed_lincont(:) = 0.
     696            dzsed_circont(:) = 0.
     697            flsed_circont(:) = 0.
     698            cfsed_circont(:) = 0.
     699            lincontfra(:)    = 0.
     700            circontfra(:)    = 0.
     701            qlincont(:)      = 0.
     702            qcircont(:)      = 0.
    676703          ENDIF
     704
     705          IF ( iftop ) THEN
     706            dzsed_abv(:) = 0.
     707            flsed_abv(:) = 0.
     708            cfsed_abv(:) = 0.
     709          ELSE
     710            dzsed_abv(:) = dzsed(:)
     711            flsed_abv(:) = flsed(:)
     712            cfsed_abv(:) = cfsed(:)
     713          ENDIF
     714          dzsed(:) = 0.
     715          flsed(:) = 0.
     716          cfsed(:) = 0.
    677717
    678718          DO i = 1, klon
    … …  
    821861                  IF (ok_ice_supersat) THEN
    822862
    823                     IF ( iftop ) THEN
    824                       cldfra_above(:) = 0.
    825                     ELSE
    826                       cldfra_above(:) = rneb(:,k+1)
    827                     ENDIF
    828 
    829863                    !---------------------------------------------
    830864                    !--   CONDENSATION AND ICE SUPERSATURATION  --
    … …  
    836870                        shear, tke_dissip(:,k), cell_area, Tbef, qtot_in, zqs, &
    837871                        gammasat, ratqs(:,k), keepgoing, pt_pron_clds, &
    838                         cldfra_above, icesed_flux,&
     872                        dzsed_abv, flsed_abv, cfsed_abv, &
     873                        dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv, &
     874                        dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv, &
     875                        dzsed, flsed, cfsed, dzsed_lincont, flsed_lincont, cfsed_lincont, &
     876                        dzsed_circont, flsed_circont, cfsed_circont, &
    839877                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
    840878                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), dcfsed(:,k), &
    … …  
    10471085
    10481086            ENDIF
     1087
     1088            ! Remove the ice that was sedimentated. As it is not included in zqn,
     1089            ! we only remove it from the total water
     1090            IF ( ok_ice_sedim ) THEN
     1091              zq(i) = zq(i) - flsed(i) / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
     1092            ENDIF
    10491093           
    10501094           
    … …  
    10791123
    10801124    IF (ok_plane_contrail) THEN
     1125
     1126      !--Ice water content of contrails
     1127      qice_lincont(:,k) = qlincont(:) - zqs(:) * lincontfra(:)
     1128      qice_circont(:,k) = qcircont(:) - zqs(:) * circontfra(:)
     1129
    10811130      !--Contrails precipitate as natural clouds. We save the partition of ice
    10821131      !--between natural clouds and contrails
    10831132      !--NB. we use qlincont / qcircont as a temporary variable to save this partition
    1084       IF ( ok_precip_lincontrails ) THEN
     1133      IF ( ok_precip_contrails ) THEN
    10851134        DO i = 1, klon
    10861135          IF ( zoliqi(i) .GT. 0. ) THEN
    1087             qlincont(i) = ( qlincont(i) - zqs(i) * lincontfra(i) ) / zoliqi(i)
     1136            qlincont(i) = qice_lincont(i,k) / zoliqi(i)
     1137            qcircont(i) = qice_circont(i,k) / zoliqi(i)
    10881138          ELSE
    10891139            qlincont(i) = 0.
     1140            qcircont(i) = 0.
    10901141          ENDIF
    10911142        ENDDO
    … …  
    10941145        !--the cloud variables
    10951146        DO i = 1, klon
    1096           qice_cont = qlincont(i) - zqs(i) * lincontfra(i)
    1097           rneb(i,k) = rneb(i,k) - lincontfra(i)
     1147          qice_cont = qice_lincont(i,k) + qice_circont(i,k)
     1148          rneb(i,k) = rneb(i,k) - ( lincontfra(i) + circontfra(i) )
    10981149          zoliq(i) = zoliq(i) - qice_cont
    10991150          zoliqi(i) = zoliqi(i) - qice_cont
    11001151        ENDDO
    11011152      ENDIF
    1102 
    1103       DO i = 1, klon
    1104         IF ( zoliqi(i) .GT. 0. ) THEN
    1105           qcircont(i) = ( qcircont(i) - zqs(i) * circontfra(i) ) / zoliqi(i)
    1106         ELSE
    1107           qcircont(i) = 0.
    1108         ENDIF
    1109       ENDDO
    11101153    ENDIF
    11111154
    … …  
    11171160                            ctot_vol(:,k), ptconv(:,k), &
    11181161                            zt, zq, zoliql, zoliqi, zfice, &
    1119                             rneb(:,k), icesed_flux, znebprecipclr, znebprecipcld, &
     1162                            rneb(:,k), flsed, znebprecipclr, znebprecipcld, &
    11201163                            zrfl, zrflclr, zrflcld, &
    11211164                            zifl, ziflclr, ziflcld, &
    … …  
    11341177      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    11351178                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
    1136                             zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, icesed_flux, &
     1179                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, flsed, &
    11371180                            rneb(:,k), znebprecipclr, znebprecipcld, &
    11381181                            zneb, tot_zneb, zrho_up, zvelo_up, &
    … …  
    11451188    IF ( ok_plane_contrail ) THEN
    11461189      !--Contrails fraction is left unchanged, but contrails water has changed
    1147       !--We alse compute the ice content that will be seen by radiation (qradice_lincont/circont)
    1148       IF ( ok_precip_lincontrails ) THEN
     1190      !--We alse compute the ice content that will be seen by radiation
     1191      !--(qradice_lincont/circont)
     1192      IF ( ok_precip_contrails ) THEN
    11491193        DO i = 1, klon
    11501194          IF ( zoliqi(i) .GT. 0. ) THEN
    11511195            qradice_lincont(i,k) = zradocond(i) * qlincont(i)
    11521196            qlincont(i) = zqs(i) * lincontfra(i) + zoliqi(i) * qlincont(i)
     1197            qradice_circont(i,k) = zradocond(i) * qcircont(i)
     1198            qcircont(i) = zqs(i) * circontfra(i) + zoliqi(i) * qcircont(i)
    11531199          ELSE
    11541200            qradice_lincont(i,k) = 0.
    11551201            lincontfra(i) = 0.
    11561202            qlincont(i) = 0.
     1203            qradice_circont(i,k) = 0.
     1204            circontfra(i) = 0.
     1205            qcircont(i) = 0.
    11571206          ENDIF
    11581207        ENDDO
    11591208      ELSE
    1160         !--If linear contrails do not precipitate, they are put back into
     1209        !--If contrails do not precipitate, they are put back into
    11611210        !--the cloud variables
    11621211        DO i = 1, klon
    1163           qice_cont = qlincont(i) - zqs(i) * lincontfra(i)
    1164           rneb(i,k) = rneb(i,k) + lincontfra(i)
     1212          rneb(i,k) = rneb(i,k) + ( lincontfra(i) + circontfra(i) )
     1213          qice_cont = qice_lincont(i,k) + qice_circont(i,k)
    11651214          zoliq(i) = zoliq(i) + qice_cont
    11661215          zoliqi(i) = zoliqi(i) + qice_cont
    11671216          zradocond(i) = zradocond(i) + qice_cont
    11681217          zradoice(i) = zradoice(i) + qice_cont
    1169           qradice_lincont(i,k) = qice_cont
     1218          qradice_lincont(i,k) = qice_lincont(i,k)
     1219          qradice_circont(i,k) = qice_circont(i,k)
    11701220        ENDDO
    11711221      ENDIF
    1172 
    1173       DO i = 1, klon
    1174         IF ( zoliqi(i) .GT. 0. ) THEN
    1175           qradice_circont(i,k) = zradocond(i) * qcircont(i)
    1176           qcircont(i) = zqs(i) * circontfra(i) + zoliqi(i) * qcircont(i)
    1177         ELSE
    1178           qradice_circont(i,k) = 0.
    1179           circontfra(i) = 0.
    1180           qcircont(i) = 0.
    1181         ENDIF
    1182       ENDDO
    11831222    ENDIF
    11841223
    … …  
    12891328      qtl_seri(:,k) = qlincont(:)
    12901329      qtc_seri(:,k) = qcircont(:)
    1291       qice_lincont(:,k) = qlincont(:) - zqs(:) * lincontfra(:)
    1292       qice_circont(:,k) = qcircont(:) - zqs(:) * circontfra(:)
    12931330    ENDIF
    12941331
    … …  
    14451482  IF ( ok_ice_sedim ) THEN
    14461483    DO i = 1, klon
    1447       snow(i) = snow(i) + icesed_flux(i)
     1484      snow(i) = snow(i) + flsed(i)
    14481485    ENDDO
    14491486  ENDIF

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

    r5684 r5691  
    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       cldfra_above, icesed_flux, &
     99      dzsed_abv, flsed_abv, cfsed_abv_in, &
     100      dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv, &
     101      dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv, &
     102      dzsed, flsed, cfsed, dzsed_lincont, flsed_lincont, cfsed_lincont, &
     103      dzsed_circont, flsed_circont, cfsed_circont, &
    100104      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, dcf_sed, &
    101105      dqi_adj, dqi_sub, dqi_con, dqi_mix, dqi_sed, &
    … …  
    131135USE lmdz_lscp_ini,   ONLY: coef_mixing_lscp, coef_shear_lscp
    132136USE lmdz_lscp_ini,   ONLY: aspect_ratio_cirrus, cooling_rate_ice_thresh
     137USE lmdz_lscp_ini,   ONLY: ok_ice_sedim, fallice_sedim
    133138
    134139USE lmdz_lscp_ini,   ONLY: ok_plane_contrail, aspect_ratio_lincontrails
    135140USE lmdz_lscp_ini,   ONLY: coef_mixing_lincontrails, coef_shear_lincontrails
    136141USE lmdz_lscp_ini,   ONLY: chi_mixing_lincontrails, linear_contrails_lifetime
     142USE lmdz_lscp_ini,   ONLY: fallice_linear_contrails, fallice_cirrus_contrails
    137143USE lmdz_aviation,   ONLY: contrails_formation
    138144
    … …  
    163169LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: keepgoing     ! .TRUE. if a new condensation loop should be computed
    164170LOGICAL,  INTENT(IN)   , DIMENSION(klon) :: pt_pron_clds  ! .TRUE. if clouds are prognostic in this mesh
    165 REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfra_above  ! cloud fraction IN THE LAYER ABOVE [-]
    166 REAL,     INTENT(IN)   , DIMENSION(klon) :: icesed_flux   ! sedimentated ice flux [kg/s/m2]
     171REAL,     INTENT(IN)   , DIMENSION(klon) :: dzsed_abv     ! sedimentated cloud height IN THE LAYER ABOVE [m]
     172REAL,     INTENT(IN)   , DIMENSION(klon) :: flsed_abv     ! sedimentated ice flux FROM THE LAYER ABOVE [kg/s/m2]
     173REAL,     INTENT(IN)   , DIMENSION(klon) :: cfsed_abv_in  ! cloud fraction IN THE LAYER ABOVE [-]
     174REAL,     INTENT(IN)   , DIMENSION(klon) :: dzsed_lincont_abv   ! sedimentated linear contrails height IN THE LAYER ABOVE [m]
     175REAL,     INTENT(IN)   , DIMENSION(klon) :: flsed_lincont_abv   ! sedimentated ice flux in linear contrails FROM THE LAYER ABOVE [kg/s/m2]
     176REAL,     INTENT(IN)   , DIMENSION(klon) :: cfsed_lincont_abv   ! linear contrails fraction IN THE LAYER ABOVE [-]
     177REAL,     INTENT(IN)   , DIMENSION(klon) :: dzsed_circont_abv   ! sedimentated contrails cirrus height IN THE LAYER ABOVE [m]
     178REAL,     INTENT(IN)   , DIMENSION(klon) :: flsed_circont_abv   ! sedimentated ice flux in contrails cirrus FROM THE LAYER ABOVE [kg/s/m2]
     179REAL,     INTENT(IN)   , DIMENSION(klon) :: cfsed_circont_abv   ! contrails cirrus fraction IN THE LAYER ABOVE [-]
     180REAL,     INTENT(INOUT), DIMENSION(klon) :: dzsed         ! sedimentated cloud height [m]
     181REAL,     INTENT(INOUT), DIMENSION(klon) :: flsed         ! sedimentated ice flux [kg/s/m2]
     182REAL,     INTENT(INOUT), DIMENSION(klon) :: cfsed         ! sedimentated cloud fraction [-]
     183REAL,     INTENT(INOUT), DIMENSION(klon) :: dzsed_lincont ! sedimentated linear contrails height [m]
     184REAL,     INTENT(INOUT), DIMENSION(klon) :: flsed_lincont ! sedimentated ice flux in linear contrails [kg/s/m2]
     185REAL,     INTENT(INOUT), DIMENSION(klon) :: cfsed_lincont ! sedimentated linear contrails fraction [-]
     186REAL,     INTENT(INOUT), DIMENSION(klon) :: dzsed_circont ! sedimentated contrails cirrus height [m]
     187REAL,     INTENT(INOUT), DIMENSION(klon) :: flsed_circont ! sedimentated ice flux in contrails cirrus [kg/s/m2]
     188REAL,     INTENT(INOUT), DIMENSION(klon) :: cfsed_circont ! sedimentated contrails cirrus fraction [-]
    167189!
    168190! Input for aviation
    … …  
    266288!
    267289! for sedimentation
    268 REAL, DIMENSION(klon) :: qice_sedim
    269 REAL :: clrfra_sed
     290REAL, DIMENSION(klon) :: cfsed_abv, qised_abv
     291REAL :: qice_sedim
     292REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3
    270293!
    271294! for mixing
    … …  
    284307! for contrails
    285308REAL :: contrails_conversion_factor
     309REAL, DIMENSION(klon) :: qised_lincont_abv, qised_circont_abv
    286310
    287311qzero(:) = 0.
    … …  
    375399      dqvc_sed(i) = 0.
    376400
    377       IF ( icesed_flux(i) .GT. 0. ) THEN
     401      qised_abv(i) = 0.
     402      cfsed_abv(i) = cfsed_abv_in(i)
     403      IF ( dzsed_abv(i) .GT. eps ) THEN
    378404        !--If ice sedimentation is activated, the quantity of sedimentated ice was added
    379405        !--to the total water vapor in the precipitation routine. Here we remove it
    380406        !--(it will be reincluded later)
    381         qice_sedim(i) = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
    382         qclr(i) = qclr(i) - qice_sedim(i)
     407        qised_abv(i) = flsed_abv(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
     408        qclr(i) = qclr(i) - qised_abv(i)
    383409      ENDIF
    384410
    … …  
    487513        ENDIF
    488514
     515        IF ( dzsed_lincont_abv(i) .GT. eps ) THEN
     516          qised_lincont_abv(i) = flsed_lincont_abv(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
     517          qised_abv(i) = qised_abv(i) - qised_lincont_abv(i)
     518          cfsed_abv(i) = cfsed_abv(i) - cfsed_lincont_abv(i)
     519        ENDIF
     520
     521        IF ( dzsed_circont_abv(i) .GT. eps ) THEN
     522          qised_circont_abv(i) = flsed_circont_abv(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
     523          qised_abv(i) = qised_abv(i) - qised_circont_abv(i)
     524          cfsed_abv(i) = cfsed_abv(i) - cfsed_circont_abv(i)
     525        ENDIF
     526
    489527        dcfl_ini(i) = 0.
    490528        dqil_ini(i) = 0.
    … …  
    13141352
    13151353
    1316       !---------------------------------------
    1317       !--         ICE SEDIMENTATION         --
    1318       !---------------------------------------
    1319       !
    1320       !--If ice supersaturation is activated, the cloud properties are prognostic.
    1321       !--The falling ice is then partly considered a new cloud in the gridbox.
    1322       !--BEWARE with this parameterisation, we can create a new cloud with a much
    1323       !--different ice water content and water vapor content than the existing cloud
    1324       !--(if it exists). This implies than unphysical fluxes of ice and vapor
    1325       !--occur between the existing cloud and the newly formed cloud (from sedimentation).
    1326       !--Note also that currently, the precipitation scheme assume a maximum
    1327       !--random overlap, meaning that the new formed clouds will not be affected
    1328       !--by vertical wind shear.
    1329       !
    1330       IF ( icesed_flux(i) .GT. 0. ) THEN
    1331 
    1332         clrfra_sed = MIN(clrfra(i), cldfra_above(i) - cldfra(i))
    1333 
    1334         IF ( ( clrfra_sed .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN
    1335 
    1336           qiceincld = qice_sedim(i) / cldfra_above(i)
    1337           IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
    1338             tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
    1339             qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) )
    1340           ELSE
    1341             qvapinclr_lim = qsat(i) - qiceincld
     1354      IF ( ok_ice_sedim ) THEN
     1355        !---------------------------------------
     1356        !--         ICE SEDIMENTATION         --
     1357        !---------------------------------------
     1358        !
     1359        !--First, the current ice is sedimentated into the layer below. The ice fallspeed
     1360        !--velocity is calculated and the quantity of sedimentated ice is computed
     1361        !--accordingly. The decrease in cloud fraction is calculated such that
     1362        !--in-cloud ice water content is constant.
     1363        !
     1364        qice_sedim = 0.
     1365        IF ( cldfra(i) .GT. eps ) THEN
     1366          dzsed(i) = MIN(fallice_sedim * dtime, dz)
     1367          cfsed(i) = cldfra(i)
     1368          qice_sedim = ( qcld(i) - qvc(i) ) * dzsed(i) / dz
     1369          !--Tendencies
     1370          dcf_sed(i) = - cldfra(i) * dzsed(i) / dz
     1371          dqi_sed(i) = - qice_sedim
     1372          dqvc_sed(i) = - qvc(i) * dzsed(i) / dz
     1373          !--Convert to flux
     1374          flsed(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime
     1375        ENDIF
     1376        !
     1377        !--Then, the sedimentated ice from above is added to the gridbox
     1378        !--The falling ice is partly considered a new cloud in the gridbox.
     1379        !--BEWARE with this parameterisation, we can create a new cloud with a much
     1380        !--different ice water content and water vapor content than the existing cloud
     1381        !--(if it exists). This implies than unphysical fluxes of ice and vapor
     1382        !--occur between the existing cloud and the newly formed cloud (from sedimentation).
     1383        !--Note also that currently, the precipitation scheme assume a maximum
     1384        !--random overlap, meaning that the new formed clouds will not be affected
     1385        !--by vertical wind shear.
     1386        !
     1387        sedfra_abv = MIN(dzsed_abv(i), dz) / dz * cfsed_abv(i)
     1388        IF ( sedfra_abv .GT. eps ) THEN
     1389
     1390          !--Three regions to be defined: (1) clear-sky, (2) cloudy, and
     1391          !--(3) region where it was previously cloudy but now clear-sky because of
     1392          !--ice sedimentation
     1393          !--(1) we use the clear-sky PDF to determine the humidity and increase the
     1394          !--cloud fraction / sublimate falling ice
     1395          !--(2) we do not increase cloud fraction and add the falling ice to the cloud
     1396          !--(3) we increase the cloud fraction but use in-cloud water vapor as
     1397          !--background vapor
     1398          sedfra2 = MIN(cfsed(i), cfsed_abv(i)) &
     1399                  * MAX(0., MIN(dz, dzsed_abv(i)) - dzsed(i)) / dz
     1400          sedfra3 = MIN(cfsed(i), cfsed_abv(i)) &
     1401                  * MIN(MIN(dz, dzsed_abv(i)), dzsed(i)) / dz
     1402          sedfra1 = sedfra_abv - sedfra3 - sedfra2
     1403
     1404          qiceincld = qised_abv(i) / sedfra_abv
     1405
     1406          !--For region (1)
     1407          IF ( ( sedfra1 .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN
     1408
     1409            IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN
     1410              tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
     1411              qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) )
     1412            ELSE
     1413              qvapinclr_lim = qsat(i) - qiceincld
     1414            ENDIF
     1415
     1416            rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
     1417            pdf_x = qvapinclr_lim / qsatl(i) * 100.
     1418            pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
     1419            pdf_x = qsat(i) / qsatl(i) * 100.
     1420            pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
     1421            IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
     1422              pdf_fra_above_lim = 0.
     1423              pdf_q_above_lim = 0.
     1424            ELSEIF ( pdf_y .LT. -10. ) THEN
     1425              pdf_fra_above_lim = clrfra(i)
     1426              pdf_q_above_lim = qclr(i)
     1427            ELSE
     1428              pdf_y = EXP( pdf_y )
     1429              pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
     1430              pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
     1431              pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
     1432              pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
     1433                                + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
     1434            ENDIF
     1435
     1436            IF ( pdf_fra_above_lim .GT. eps ) THEN
     1437              sedfra1 = sedfra1 * pdf_fra_above_lim / clrfra(i)
     1438              dcf_sed(i)  = dcf_sed(i)  + sedfra1
     1439              dqvc_sed(i) = dqvc_sed(i) + sedfra1 * pdf_q_above_lim / pdf_fra_above_lim
     1440              dqi_sed(i)  = dqi_sed(i)  + sedfra1 * qiceincld
     1441            ENDIF
     1442          ENDIF ! ( sedfra1 .GT. eps .AND. ( clrfra(i) .GT. eps )
     1443
     1444          !--For region (2)
     1445          dqi_sed(i) = dqi_sed(i) + sedfra2 * qiceincld
     1446
     1447          !--For region (3)
     1448          IF ( sedfra3 .GT. eps ) THEN
     1449            dcf_sed(i)  = dcf_sed(i)  + sedfra3
     1450            dqvc_sed(i) = dqvc_sed(i) + sedfra3 * qvc(i) / cldfra(i)
     1451            dqi_sed(i)  = dqi_sed(i)  + sedfra3 * qiceincld
    13421452          ENDIF
    1343 
    1344           rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100.
    1345           pdf_x = qvapinclr_lim / qsatl(i) * 100.
    1346           pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i)
    1347           pdf_x = qsat(i) / qsatl(i) * 100.
    1348           pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i)
    1349           IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows
    1350             pdf_fra_above_lim = 0.
    1351             pdf_q_above_lim = 0.
    1352           ELSEIF ( pdf_y .LT. -10. ) THEN
    1353             pdf_fra_above_lim = clrfra(i)
    1354             pdf_q_above_lim = qclr(i)
    1355           ELSE
    1356             pdf_y = EXP( pdf_y )
    1357             pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )
    1358             pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i)
    1359             pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i)
    1360             pdf_q_above_lim = ( pdf_e3 * clrfra(i) &
    1361                               + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100.
    1362           ENDIF
    1363 
    1364           IF ( pdf_fra_above_lim .GT. eps ) THEN
    1365             dcf_sed(i)  = clrfra_sed * pdf_fra_above_lim / clrfra(i)
    1366             dqvc_sed(i) = dcf_sed(i) * pdf_q_above_lim / pdf_fra_above_lim
    1367           ENDIF
    1368           !--We include the sedimentated ice:
    1369           dqi_sed(i)  = qiceincld &    !--We include the sedimentated ice:
    1370                       * ( dcf_sed(i) & !--the part that falls in the newly formed cloud and
    1371                         + cldfra(i) )  !--the part that falls in the existing cloud
    1372 
    1373         ELSE
    1374 
    1375           dqi_sed(i) = qice_sedim(i)
    1376 
    1377         ENDIF ! ( clrfra_sed .GT. eps .AND. ( clrfra(i) .GT. eps )
     1453        ENDIF ! qised_abv(i) .GT. 0.
     1454        !PRINT *, 'A', cldfra(i), dcf_sed(i), clrfra(i)
     1455        !PRINT *, 'B', qcld(i) - qvc(i), qvc(i), dqvc_sed(i), dqi_sed(i)
    13781456
    13791457        !--Add tendencies
    … …  
    13831461        clrfra(i)  = MAX(0., clrfra(i) - dcf_sed(i))
    13841462        !--We re-include sublimated sedimentated ice into clear sky water vapor
    1385         qclr(i)    = qclr(i) - dqvc_sed(i) + qice_sedim(i) - dqi_sed(i)
    1386 
    1387       ENDIF ! icesed_flux(i) .GT. 0.
     1463        qclr(i)    = qclr(i) - dqvc_sed(i) + qised_abv(i) - dqi_sed(i)
     1464
     1465      ENDIF ! ok_ice_sedim
    13881466
    13891467

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

    r5643 r5691  
    228228  !$OMP THREADPRIVATE(ok_plane_contrail)
    229229
    230   LOGICAL, SAVE, PROTECTED :: ok_precip_lincontrails=.TRUE.  ! if True, linear contrails can precipitate
    231   !$OMP THREADPRIVATE(ok_precip_lincontrails)
     230  LOGICAL, SAVE, PROTECTED :: ok_precip_contrails=.TRUE.     ! if True, contrails can be autoconverted to snow
     231  !$OMP THREADPRIVATE(ok_precip_contrails)
    232232
    233233  REAL, SAVE, PROTECTED :: aspect_ratio_lincontrails=.1      ! [-] aspect ratio of linear contrails
    … …  
    260260  REAL, SAVE, PROTECTED :: initial_height_contrails=200.     ! [m] initial height of the linear contrails formed
    261261  !$OMP THREADPRIVATE(initial_height_contrails)
     262
     263  REAL, SAVE, PROTECTED :: fallice_linear_contrails=1.       ! [m/s] Ice fallspeed velocity in linear contrails
     264  !$OMP THREADPRIVATE(fallice_linear_contrails)
     265
     266  REAL, SAVE, PROTECTED :: fallice_cirrus_contrails=1.       ! [m/s] Ice fallspeed velocity in cirrus contrails
     267  !$OMP THREADPRIVATE(fallice_cirrus_contrails)
    262268
    263269  REAL, SAVE, PROTECTED :: aviation_coef=1.                  ! [-] scaling factor for aviation emissions and flown distance
    … …  
    512518    CALL getin_p('coef_shear_lscp',coef_shear_lscp)
    513519    ! for aviation
    514     CALL getin_p('ok_precip_lincontrails',ok_precip_lincontrails)
     520    CALL getin_p('ok_precip_contrails',ok_precip_contrails)
    515521    CALL getin_p('aspect_ratio_lincontrails',aspect_ratio_lincontrails)
    516522    coef_mixing_lincontrails=coef_mixing_lscp
    … …  
    525531    CALL getin_p('initial_width_contrails',initial_width_contrails)
    526532    CALL getin_p('initial_height_contrails',initial_height_contrails)
     533    CALL getin_p('fallice_linear_contrails',fallice_linear_contrails)
     534    CALL getin_p('fallice_cirrus_contrails',fallice_cirrus_contrails)
    527535    CALL getin_p('aviation_coef',aviation_coef)
    528536
    … …  
    615623    WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
    616624    ! for aviation
    617     WRITE(lunout,*) 'lscp_ini, ok_precip_lincontrails:', ok_precip_lincontrails
     625    WRITE(lunout,*) 'lscp_ini, ok_precip_contrails:', ok_precip_contrails
    618626    WRITE(lunout,*) 'lscp_ini, aspect_ratio_lincontrails:', aspect_ratio_lincontrails
    619627    WRITE(lunout,*) 'lscp_ini, coef_mixing_lincontrails:', coef_mixing_lincontrails
    … …  
    626634    WRITE(lunout,*) 'lscp_ini, initial_width_contrails:', initial_width_contrails
    627635    WRITE(lunout,*) 'lscp_ini, initial_height_contrails:', initial_height_contrails
     636    WRITE(lunout,*) 'lscp_ini, fallice_linear_contrails:', fallice_linear_contrails
     637    WRITE(lunout,*) 'lscp_ini, fallice_cirrus_contrails:', fallice_cirrus_contrails
    628638    WRITE(lunout,*) 'lscp_ini, aviation_coef:', aviation_coef
    629639

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

    r5613 r5691  
    126126!--------------------
    127127IF ( ok_ice_sedim ) THEN
     128
     129  zcpeau = RCPD * RVTMP2
     130
    128131  DO i =1, klon
    129132    !--If the sedimentation of ice crystals is activated, the falling ice is sublimated and
    … …  
    131134    IF ( icesed_flux(i) .GT. 0. ) THEN
    132135      qice_sedim = icesed_flux(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime
     136
     137      !--Thermalisation of sedimentated ice
     138      zcpair = RCPD * ( 1. + RVTMP2 * zq(i) )
     139      zt(i) = ( ztupnew(i) * qice_sedim * zcpeau + zcpair * zt(i) ) &
     140            / ( zcpair + qice_sedim * zcpeau )
     141
    133142      !--Vapor is updated after evaporation/sublimation (it is increased)
    134143      zq(i) = zq(i) + qice_sedim
    … …  
    384393REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflclr        !--flux of snow gridbox-mean in clear sky [kg/s/m2]
    385394REAL,    INTENT(INOUT), DIMENSION(klon) :: ziflcld        !--flux of snow gridbox-mean in cloudy air [kg/s/m2]
    386 REAL,    INTENT(OUT),  DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
     395REAL,    INTENT(INOUT), DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
    387396
    388397REAL,    INTENT(OUT),   DIMENSION(klon) :: zradocond      !--condensed water used in the radiation scheme [kg/kg]
    … …  
    429438ziflprev(:) = 0.
    430439
    431 icesed_flux(:) = 0.
     440! AB ICE SEDIM
     441!icesed_flux(:) = 0.
    432442
    433443
    … …  
    691701!--                  ICE SEDIMENTATION
    692702!---------------------------------------------------------
    693 IF ( ok_ice_sedim ) THEN
     703!IF ( ok_ice_sedim ) THEN
     704IF ( .FALSE. ) THEN ! AB ICE SEDIM
    694705
    695706  DO i = 1, klon
    … …  
    908919  ENDDO
    909920
     921ENDIF
     922
     923!--------------------------------------
     924!-- THERMALISATION OF ICE SEDIMENTATION
     925!--------------------------------------
     926
     927!--If the sedimentation of ice crystals is activated, the falling ice is thermzalised
     928!--to the temperature of the gridbox
     929IF ( ok_ice_sedim ) THEN
     930  cpw = RCPD * RVTMP2
     931  DO i = 1, klon
     932    IF ( icesed_flux(i) .GT. 0. ) THEN
     933      qice_sedim = icesed_flux(i) / dhum_to_dflux(i)
     934      !--No condensed water so cp=cp(vapor+dry air)
     935      !-- RVTMP2=rcpv/rcpd-1
     936      cpair = RCPD * ( 1. + RVTMP2 * qvap(i) )
     937      !-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer
     938      temp(i) = ( tempupnew(i) * qice_sedim * cpw + cpair * temp(i) ) &
     939              / ( cpair + qice_sedim * cpw )
     940    ENDIF ! ok_ice_sedim
     941  ENDDO
    910942ENDIF
    911943
    … …  
    14241456REAL,    INTENT(INOUT), DIMENSION(klon) :: snowcld        !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2]
    14251457
    1426 REAL,    INTENT(OUT),  DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
     1458REAL,    INTENT(INOUT), DIMENSION(klon) :: icesed_flux    !--flux of sedimentated ice [kg/s/m2]
    14271459REAL,    INTENT(OUT),   DIMENSION(klon) :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
    14281460REAL,    INTENT(OUT),   DIMENSION(klon) :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
    … …  
    14871519
    14881520qzero(:)    = 0.
    1489 icesed_flux(:) = 0.
     1521! AB ICE SEDIM
     1522!icesed_flux(:) = 0.
    14901523
    14911524dqrcol(:)   = 0.
    … …  
    19601993  !--                  ICE SEDIMENTATION
    19611994  !---------------------------------------------------------
    1962   IF ( ok_ice_sedim ) THEN
     1995  !IF ( ok_ice_sedim ) THEN
     1996  IF ( .FALSE. ) THEN ! AB ICE SEDIM
    19631997
    19641998    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. eps ) ) THEN
Note: See TracChangeset for help on using the changeset viewer.