Ignore:
Timestamp:
May 26, 2025, 5:20:01 PM (6 weeks ago)
Author:
aborella
Message:

Various small adjustments

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

Legend:

Unmodified
Added
Removed

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

    r5641 r5684  
    353353        (delPhase(tracers(:)%gen0Name)     == 'CLDFRA')))
    354354!   IF(COUNT(tracers%iso_iName == 0) - COUNT(delPhase(tracers%name) == 'H2O' .AND. tracers%component == 'lmdz') /= nqtottr) &
    355    IF(COUNT(tracers%iso_iName == 0) - COUNT(tracers(:)%component == 'lmdz' .AND. &
    356        ((delPhase(tracers(:)%name)     == 'H2O') .OR. &
    357         (delPhase(tracers(:)%name)     == 'CLDFRA'))) /= nqtottr) &
    358       CALL abort_physic(modname, 'problem with the computation of nqtottr', 1)
     355   !IF(COUNT(tracers%iso_iName == 0) - COUNT(tracers(:)%component == 'lmdz' .AND. &
     356   !    ((delPhase(tracers(:)%name)     == 'H2O') .OR. &
     357   !     (delPhase(tracers(:)%name)     == 'CLDFRA'))) /= nqtottr) &
     358   !   CALL abort_physic(modname, 'problem with the computation of nqtottr', 1)
    359359
    360360   !--- Compute indices for water
    … …  
    370370   iqtc = strIdx(tracers(:)%name, 'CIRCONTWATER')
    371371   !--Two ways of declaring tracers - the way below should be deleted in the future
    372    IF (icf.EQ.0) icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
     372   IF (icf.EQ.0)  icf  = strIdx(tracers(:)%name, addPhase('H2O', 'f'))
    373373   IF (iqvc.EQ.0) iqvc = strIdx(tracers(:)%name, addPhase('H2O', 'c'))
    374374   IF (icfl.EQ.0) icfl = strIdx(tracers(:)%name, addPhase('H2O', 'z'))

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

    r5641 r5684  
    599599    REAL, INTENT(IN)    :: pplay(klon, klev) ! mid-layer pressure [Pa]
    600600    REAL, INTENT(IN)    :: temp(klon, klev) ! temperature [K]
    601     REAL, INTENT(OUT)   :: flight_dist(klon,klev) ! Aviation distance flown within the mesh [m/s/mesh]
    602     REAL, INTENT(OUT)   :: flight_h2o(klon,klev)  ! Aviation H2O emitted within the mesh [kgH2O/s/mesh]
     601    REAL, INTENT(OUT)   :: flight_dist(klon,klev) ! Aviation distance flown concentration [m/s/m3]
     602    REAL, INTENT(OUT)   :: flight_h2o(klon,klev)  ! Aviation H2O emitted concentration [kgH2O/s/m3]
    603603
    604604    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

    r5644 r5684  
    2525     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
    2626     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
    27      cfl_seri, cfc_seri, qtl_seri, qtc_seri,flight_dist,&
    28      flight_h2o, qice_lincont, qice_circont, Tcritcont, &
    29      qcritcont, potcontfraP, potcontfraNP,              &
     27     cfl_seri, cfc_seri, qtl_seri, qtc_seri,            &
     28     qice_lincont, qice_circont, flight_dist,           &
     29     flight_h2o, qradice_lincont, qradice_circont,      &
     30     Tcritcont, qcritcont, potcontfraP, potcontfraNP,   &
    3031     cloudth_sth,                                       &
    3132     cloudth_senv, cloudth_sigmath, cloudth_sigmaenv,   &
    … …  
    132133USE phys_local_var_mod, ONLY : dcfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix
    133134USE phys_local_var_mod, ONLY : dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix
     135USE geometry_mod, ONLY: longitude_deg, latitude_deg
    134136
    135137IMPLICIT NONE
    … …  
    253255  ! for contrails and aviation
    254256
    255   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_lincont   !--condensed water in linear contrails used in the radiation scheme [kg/kg]
    256   REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_circont   !--condensed water in contrail cirrus used in the radiation scheme [kg/kg]
     257  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_lincont   !--condensed water in linear contrails [kg/kg]
     258  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_circont   !--condensed water in contrail cirrus [kg/kg]
     259  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qradice_lincont!--condensed water in linear contrails used in the radiation scheme [kg/kg]
     260  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qradice_circont!--condensed water in contrail cirrus used in the radiation scheme [kg/kg]
    257261  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
    258262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
    … …  
    333337  ! for condensation and ice supersaturation
    334338  REAL, DIMENSION(klon) :: qvc, qvcl, shear
    335   REAL :: delta_z
     339  REAL :: delta_z, deepconv_coef
    336340  ! for contrails
    337341  REAL, DIMENSION(klon) :: lincontfra, circontfra, qlincont, qcircont
    … …  
    673677
    674678          DO i = 1, klon
    675             pt_pron_clds(i) = ( ( ( ( zt(i) .LE. temp_nowater ) .OR. ok_weibull_warm_clouds ) &
    676                 .AND. ( .NOT. ok_no_issr_strato .OR. ( stratomask(i,k) .EQ. 0. ) ) ) &
    677                 .AND. ( cfcon(i,k) .LT. ( 1. - eps ) ) )
     679            pt_pron_clds(i) = ( cfcon(i,k) .LT. ( 1. - eps ) )
    678680          ENDDO
     681          IF ( .NOT. ok_weibull_warm_clouds ) THEN
     682            DO i = 1, klon
     683              pt_pron_clds(i) = pt_pron_clds(i) .AND. ( zt(i) .LE. temp_nowater )
     684            ENDDO
     685          ENDIF
     686          IF ( ok_no_issr_strato ) THEN
     687            DO i = 1, klon
     688              pt_pron_clds(i) = pt_pron_clds(i) .AND. ( stratomask(i,k) .EQ. 0. )
     689            ENDDO
     690          ENDIF
    679691
    680692          totfra_in(:) = 1.
    … …  
    704716                !--Else if deep convection is strengthening, it consumes the existing cloud
    705717                !--fraction (which does not at this moment represent deep convection)
    706                 cldfra_in(i) = cldfra_in(i) * ( 1. &
    707                              - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
    708                 qvc_in(i)    = qvc_in(i)    * ( 1. &
    709                              - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
    710                 qice_in(i)   = qice_in(i)   * ( 1. &
    711                              - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) ) )
     718                deepconv_coef = 1. - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) )
     719                cldfra_in(i) = cldfra_in(i) * deepconv_coef
     720                qvc_in(i)    = qvc_in(i)    * deepconv_coef
     721                qice_in(i)   = qice_in(i)   * deepconv_coef
     722                IF ( ok_plane_contrail ) THEN
     723                  !--If contrails are activated, their fraction is also reduced when deep
     724                  !--convection is active
     725                  cfl_seri(i,k) = cfl_seri(i,k) * deepconv_coef
     726                  qtl_seri(i,k) = qtl_seri(i,k) * deepconv_coef
     727                  cfc_seri(i,k) = cfc_seri(i,k) * deepconv_coef
     728                  qtc_seri(i,k) = qtc_seri(i,k) * deepconv_coef
     729                ENDIF
    712730              ENDIF
    713 
    714               !--Barriers
    715               cldfra_in(i) = MAX(0., MIN(totfra_in(i), cldfra_in(i)))
    716               qvc_in(i)    = MAX(0., MIN(qtot_in(i), qvc_in(i)))
    717               qice_in(i)   = MAX(0., MIN(qtot_in(i) - qvc_in(i), qice_in(i)))
    718731
    719732              !--Calculate the shear value (input for condensation and ice supersat)
    … …  
    11321145    IF ( ok_plane_contrail ) THEN
    11331146      !--Contrails fraction is left unchanged, but contrails water has changed
    1134       !--We alse compute the ice content that will be seen by radiation (qice_lincont/circont)
     1147      !--We alse compute the ice content that will be seen by radiation (qradice_lincont/circont)
    11351148      IF ( ok_precip_lincontrails ) THEN
    11361149        DO i = 1, klon
    11371150          IF ( zoliqi(i) .GT. 0. ) THEN
    1138             qice_lincont(i,k) = zradocond(i) * qlincont(i)
     1151            qradice_lincont(i,k) = zradocond(i) * qlincont(i)
    11391152            qlincont(i) = zqs(i) * lincontfra(i) + zoliqi(i) * qlincont(i)
    11401153          ELSE
    1141             qice_lincont(i,k) = 0.
     1154            qradice_lincont(i,k) = 0.
    11421155            lincontfra(i) = 0.
    11431156            qlincont(i) = 0.
    … …  
    11541167          zradocond(i) = zradocond(i) + qice_cont
    11551168          zradoice(i) = zradoice(i) + qice_cont
    1156           qice_lincont(i,k) = qice_cont
     1169          qradice_lincont(i,k) = qice_cont
    11571170        ENDDO
    11581171      ENDIF
    … …  
    11601173      DO i = 1, klon
    11611174        IF ( zoliqi(i) .GT. 0. ) THEN
    1162           qice_circont(i,k) = zradocond(i) * qcircont(i)
     1175          qradice_circont(i,k) = zradocond(i) * qcircont(i)
    11631176          qcircont(i) = zqs(i) * circontfra(i) + zoliqi(i) * qcircont(i)
    11641177        ELSE
    1165           qice_circont(i,k) = 0.
     1178          qradice_circont(i,k) = 0.
    11661179          circontfra(i) = 0.
    11671180          qcircont(i) = 0.
    … …  
    12611274    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
    12621275    CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
    1263    
     1276
    12641277    !--AB Write diagnostics and tracers for ice supersaturation
    12651278    IF ( ok_plane_contrail ) THEN
     1279      DO i = 1, klon
     1280        IF ( zoliq(i) .LE. 0. ) THEN
     1281          lincontfra(i) = 0.
     1282          circontfra(i) = 0.
     1283          qlincont(i) = 0.
     1284          qcircont(i) = 0.
     1285        ENDIF
     1286      ENDDO
    12661287      cfl_seri(:,k) = lincontfra(:)
    12671288      cfc_seri(:,k) = circontfra(:)
    12681289      qtl_seri(:,k) = qlincont(:)
    12691290      qtc_seri(:,k) = qcircont(:)
     1291      qice_lincont(:,k) = qlincont(:) - zqs(:) * lincontfra(:)
     1292      qice_circont(:,k) = qcircont(:) - zqs(:) * circontfra(:)
    12701293    ENDIF
    12711294
    … …  
    12931316        ENDIF
    12941317
     1318        !--Deep convection clouds properties are not advected
     1319        IF ( ptconv(i,k) .AND. pt_pron_clds(i) .AND. ok_nodeep_lscp ) THEN
     1320          cf_seri(i,k) = MAX(0., cf_seri(i,k) - cfcon(i,k))
     1321          qvc_seri(i,k) = MAX(0., qvc_seri(i,k) - qvcon_old(i,k) * cfcon(i,k))
     1322          zoliq(i) = MAX(0., zoliq(i) - qccon_old(i,k) * cfcon(i,k))
     1323          zoliqi(i) = MAX(0., zoliqi(i) - qccon_old(i,k) * cfcon(i,k))
     1324        ENDIF
    12951325        !--Deep convection clouds properties are removed from radiative properties
    12961326        !--outputed from lscp (NB. rneb and radocond are only used for the radiative

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

    r5643 r5684  
    255255!
    256256! for dissipation
    257 REAL, DIMENSION(klon) :: temp_diss, qsati_diss
    258 REAL :: pdf_shape, qiceincld_min
     257REAL, DIMENSION(klon) :: temp_diss, qsati_diss, qiceincld_min
     258REAL :: pdf_shape
    259259!
    260260! for condensation
    … …  
    292292temp_diss(:) = temp(:) + cooling_rate_ice_thresh * dtime
    293293CALL calc_qsat_ecmwf(klon, temp_diss, qzero, pplay, RTT, 2, .FALSE., qsati_diss, dqsat_tmp)
     294!--Additionally to a minimum in cloud water vapor, we impose a minimum
     295!--on the in-cloud ice water content. It is calculated following
     296!--Marti and Mauersberger (1993), see also Schiller et al. (2008)
     297qiceincld_min(:) = qsati_diss(:) - qsat(:)
    294298
    295299!
    … …  
    475479          lincontfra(i) = lincontfra(i) * mixed_fraction
    476480          circontfra(i) = circontfra(i) * mixed_fraction
    477           qlincont(i)   = qlincont(i)   * mixed_fraction
    478           qcircont(i)   = qcircont(i)   * mixed_fraction
    479481        ENDIF
    480482        mixed_fraction = qlincont(i) + qcircont(i)
    481483        IF ( mixed_fraction .GT. qcld(i) ) THEN
    482484          mixed_fraction = qcld(i) / mixed_fraction
    483           lincontfra(i) = lincontfra(i) * mixed_fraction
    484           circontfra(i) = circontfra(i) * mixed_fraction
    485           qlincont(i)   = qlincont(i)   * mixed_fraction
    486           qcircont(i)   = qcircont(i)   * mixed_fraction
     485          qlincont(i) = qlincont(i) * mixed_fraction
     486          qcircont(i) = qcircont(i) * mixed_fraction
    487487        ENDIF
    488 
    489488
    490489        dcfl_ini(i) = 0.
    … …  
    519518      !--If there is a linear contrail
    520519      IF ( lincontfra(i) .GT. eps ) THEN
    521 
    522520        !--The contrail is always adjusted to saturation
    523521        qiceincld = ( qlincont(i) / lincontfra(i) - qsat(i) )
    524522        !--If the ice water content is too low, the cloud is purely sublimated
    525         IF ( qiceincld .LT. qiceincld_min ) THEN
     523        IF ( qiceincld .LT. qiceincld_min(i) ) THEN
    526524          dcfl_sub(i) = - lincontfra(i)
    527525          dqil_sub(i) = - qiceincld * lincontfra(i)
    … …  
    532530          qclr(i) = qclr(i) - dqtl_sub(i)
    533531        ENDIF   ! qiceincld .LT. eps
    534 
    535532        !--We remove contrails from the main class
    536533        cldfra(i) = MAX(0., cldfra(i) - lincontfra(i))
    … …  
    541538      !--If there is a contrail cirrus
    542539      IF ( circontfra(i) .GT. eps ) THEN
    543 
    544540        !--The contrail is always adjusted to saturation
    545541        qiceincld = ( qcircont(i) / circontfra(i) - qsat(i) )
    546542        !--If the ice water content is too low, the cloud is purely sublimated
    547         IF ( qiceincld .LT. qiceincld_min ) THEN
     543        IF ( qiceincld .LT. qiceincld_min(i) ) THEN
    548544          dcfc_sub(i) = - circontfra(i)
    549545          dqic_sub(i) = - qiceincld * circontfra(i)
    … …  
    554550          qclr(i) = qclr(i) - dqtc_sub(i)
    555551        ENDIF   ! qiceincld .LT. eps
    556 
    557552        !--We remove contrails from the main class
    558553        cldfra(i) = MAX(0., cldfra(i) - circontfra(i))
    … …  
    601596            ELSE
    602597              !--If the cloud is initially subsaturated
    603               !!--Exact explicit integration (qvc exact, qice explicit)
    604               !!--Same but depo_coef_cirrus = 1
    605               !tauinv_depsub = qiceincld**(1./3.) * kappa_depsub
    606               !qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) * EXP( - dtime * tauinv_depsub )
    607598              !--Exact explicit integration (qice exact, qvc explicit)
    608599              !--The barrier is set so that the resulting vapor in cloud
    … …  
    634625          !--    DISSIPATION OF THE CLOUD    --
    635626          !------------------------------------
    636           !--Additionally to a minimum in cloud water vapor, we impose a minimum
    637           !--on the in-cloud ice water content. It is calculated following
    638           !--Marti and Mauersberger (1993), see also Schiller et al. (2008)
    639           qiceincld_min = qsati_diss(i) - qsat(i)
    640627
    641628          !--If the dissipation process must be activated
    642           IF ( ( qvapincld_new + qiceincld_min ) .GT. qvapincld ) THEN
     629          IF ( ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) ) .GT. qvapincld ) THEN
    643630            !--Gamma distribution starting at qvapincld
    644631            pdf_shape = nu_iwc_pdf_lscp / qiceincld
    645             pdf_y = pdf_shape * ( qvapincld_new + qiceincld_min - qvapincld )
     632            pdf_y = pdf_shape * ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) - qvapincld )
    646633            pdf_e1 = GAMMAINC ( nu_iwc_pdf_lscp      , pdf_y )
    647634            pdf_e2 = GAMMAINC ( nu_iwc_pdf_lscp + 1. , pdf_y )
    … …  
    673660            clrfra(i) = MIN(totfra_in(i), clrfra(i) - dcf_sub(i))
    674661            qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i)
    675           ENDIF ! qvapincld_new .GT. qvapincld
     662          ENDIF ! ( MIN(qsat(i), qvapincld_new) + qiceincld_min(i) ) .GT. qvapincld
    676663
    677664
    … …  
    716703        !--  tuning coef * (clear sky area**0.25) * (function of temperature)
    717704        pdf_e1 = beta_pdf_lscp * ( clrfra(i) * cell_area(i) )**0.25 &
    718                 * MAX( temp(i) - temp_thresh_pdf_lscp, 0. )
     705                * MAX( MAX(205., MIN(250., temp(i))) - temp_thresh_pdf_lscp, 0. )
    719706        IF ( rhl_clr .GT. 50. ) THEN
    720707          pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp
    … …  
    722709          pdf_std = pdf_e1 * rhl_clr / 50.
    723710        ENDIF
    724         pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. )
     711        pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * &
     712                    MAX( temp_nowater - MAX(205., MIN(250., temp(i))), 0. )
    725713        pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3
    726714        pdf_alpha(i) = MIN(10., pdf_alpha(i)) !--Avoid overflows
    … …  
    830818        !-- clouds_perim = P * N_cld_mix
    831819        !--
    832         !--We assume that the ratio between b and a is a function of
    833         !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
    834         !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
    835         !--are spherical.
    836         !-- b / a = bovera = MAX(0.1, cldfra)
    837         !bovera = MAX(0.1, cldfra(i))
    838820        bovera = aspect_ratio_cirrus
    839821        !--P / a is a function of b / a only, that we can calculate
    … …  
    994976        !-- clouds_perim = P * N_cld_mix
    995977        !--
    996         !--We assume that the ratio between b and a is a function of
    997         !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
    998         !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
    999         !--are spherical.
    1000         !-- b / a = bovera = MAX(0.1, cldfra)
    1001978        bovera = aspect_ratio_lincontrails
    1002979        !--P / a is a function of b / a only, that we can calculate
    … …  
    10781055          dcfc_mix(i) = dcfc_mix(i) + clrfra_mix
    10791056          dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * qclr(i) / clrfra(i)
    1080           dqic_mix(i) = dqic_mix(i) + clrfra_mix * sigma_mix &
    1081                       * ( qclr(i) / clrfra(i) - qsat(i) )
     1057          dqic_mix(i) = dqic_mix(i) + clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) )
    10821058        ELSE
    10831059          !--We then calculate the clear sky part where the humidity is lower than
    … …  
    11451121
    11461122        ENDIF
    1147       ENDIF ! ( lincontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
     1123      ENDIF ! ( lincontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )
    11481124
    11491125      IF ( ( circontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps ) ) THEN
    … …  
    11671143        !-- clouds_perim = P * N_cld_mix
    11681144        !--
    1169         !--We assume that the ratio between b and a is a function of
    1170         !--cldfra such that it is 1 for cldfra = 1 and it is low for little cldfra, because
    1171         !--if cldfra is low the clouds are linear, and if cldfra is high, the clouds
    1172         !--are spherical.
    1173         !-- b / a = bovera = MAX(0.1, cldfra)
    11741145        bovera = aspect_ratio_cirrus
    11751146        !--P / a is a function of b / a only, that we can calculate
    … …  
    12461217          dcfc_mix(i) = dcfc_mix(i) + clrfra_mix
    12471218          dqtc_mix(i) = dqtc_mix(i) + clrfra_mix * qclr(i) / clrfra(i)
    1248           dqic_mix(i) = dqic_mix(i) + clrfra_mix * sigma_mix &
    1249                       * ( qclr(i) / clrfra(i) - qsat(i) )
     1219          dqic_mix(i) = dqic_mix(i) + clrfra_mix * ( qclr(i) / clrfra(i) - qsat(i) )
    12501220        ELSE
    12511221          !--We then calculate the clear sky part where the humidity is lower than
    … …  
    13011271
    13021272        ENDIF
    1303       ENDIF ! ( circontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT.  eps )
    1304 
    1305       !--We balance the mixing tendencies between the different cloud classes
    1306       mixed_fraction = dcf_mix(i) + dcfl_mix(i) + dcfc_mix(i)
    1307       IF ( mixed_fraction .GT. clrfra(i) ) THEN
    1308         mixed_fraction = clrfra(i) / mixed_fraction
    1309         dcf_mix(i)  = dcf_mix(i)  * mixed_fraction
    1310         dqvc_mix(i) = dqvc_mix(i) * mixed_fraction
    1311         dqi_mix(i)  = dqi_mix(i)  * mixed_fraction
    1312         dcfl_mix(i) = dcfl_mix(i) * mixed_fraction
    1313         dqtl_mix(i) = dqtl_mix(i) * mixed_fraction
    1314         dqil_mix(i) = dqil_mix(i) * mixed_fraction
    1315         dcfc_mix(i) = dcfc_mix(i) * mixed_fraction
    1316         dqtc_mix(i) = dqtc_mix(i) * mixed_fraction
    1317         dqic_mix(i) = dqic_mix(i) * mixed_fraction
    1318       ENDIF
    1319 
    1320       IF ( lincontfra(i) .GT. eps ) THEN
    1321         !--Add tendencies
    1322         lincontfra(i) = lincontfra(i) + dcfl_mix(i)
    1323         qlincont(i)   = qlincont(i) + dqtl_mix(i)
    1324         clrfra(i)  = clrfra(i) - dcfl_mix(i)
    1325         qclr(i)    = qclr(i) - dqtl_mix(i)
    1326       ENDIF
    1327 
    1328       IF ( circontfra(i) .GT. eps ) THEN
    1329         !--Add tendencies
    1330         circontfra(i) = circontfra(i) + dcfc_mix(i)
    1331         qcircont(i)   = qcircont(i) + dqtc_mix(i)
    1332         clrfra(i)  = clrfra(i) - dcfc_mix(i)
    1333         qclr(i)    = qclr(i) - dqtc_mix(i)
     1273      ENDIF ! ( circontfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps )
     1274
     1275      IF ( ( lincontfra(i) + circontfra(i) ) .GT. eps ) THEN
     1276        !--We balance the mixing tendencies between the different cloud classes
     1277        mixed_fraction = dcf_mix(i) + dcfl_mix(i) + dcfc_mix(i)
     1278        IF ( mixed_fraction .GT. clrfra(i) ) THEN
     1279          mixed_fraction = clrfra(i) / mixed_fraction
     1280          dcf_mix(i)  = dcf_mix(i)  * mixed_fraction
     1281          dqvc_mix(i) = dqvc_mix(i) * mixed_fraction
     1282          dqi_mix(i)  = dqi_mix(i)  * mixed_fraction
     1283          dcfl_mix(i) = dcfl_mix(i) * mixed_fraction
     1284          dqtl_mix(i) = dqtl_mix(i) * mixed_fraction
     1285          dqil_mix(i) = dqil_mix(i) * mixed_fraction
     1286          dcfc_mix(i) = dcfc_mix(i) * mixed_fraction
     1287          dqtc_mix(i) = dqtc_mix(i) * mixed_fraction
     1288          dqic_mix(i) = dqic_mix(i) * mixed_fraction
     1289        ENDIF
     1290
     1291        IF ( lincontfra(i) .GT. eps ) THEN
     1292          !--Add tendencies
     1293          lincontfra(i) = lincontfra(i) + dcfl_mix(i)
     1294          qlincont(i)   = qlincont(i) + dqtl_mix(i)
     1295          clrfra(i)  = clrfra(i) - dcfl_mix(i)
     1296          qclr(i)    = qclr(i) - dqtl_mix(i)
     1297        ENDIF
     1298
     1299        IF ( circontfra(i) .GT. eps ) THEN
     1300          !--Add tendencies
     1301          circontfra(i) = circontfra(i) + dcfc_mix(i)
     1302          qcircont(i)   = qcircont(i) + dqtc_mix(i)
     1303          clrfra(i)  = clrfra(i) - dcfc_mix(i)
     1304          qclr(i)    = qclr(i) - dqtc_mix(i)
     1305        ENDIF
    13341306      ENDIF
    13351307
    … …  
    14521424      ENDIF
    14531425
     1426
    14541427      !-------------------------------------------
    14551428      !--       FINAL BARRIERS AND OUTPUTS      --
    … …  
    15091482
    15101483      !--Convert existing contrail fraction into "natural" cirrus cloud fraction
    1511       IF ( ( cldfra(i) .GE. ( totfra_in(i) - eps ) ) .OR. ( lincontfra(i) .LE. eps ) ) THEN
     1484      IF ( ( clrfra(i) .LE. eps ) .OR. ( lincontfra(i) .LE. eps ) ) THEN
    15121485        contrails_conversion_factor = 1.
    15131486      ELSE
    … …  
    15181491            !--cannot exist. The exponent is set so that this only happens for
    15191492            !--very cloudy gridboxes
    1520             * ( 1. - cldfra(i) / totfra_in(i) )**0.1 )
     1493            * ( clrfra(i) / totfra_in(i) )**0.1 )
    15211494      ENDIF
    15221495      dcfl_cir(i) = - contrails_conversion_factor * lincontfra(i)
    15231496      dqtl_cir(i) = - contrails_conversion_factor * qlincont(i)
    15241497
    1525       dcfl_ini(i) = MIN(MIN(dcfl_ini(i), issrfra(i)), totfra_in(i) - cldfra(i))
    1526       dqtl_ini(i) = MIN(MIN(dqtl_ini(i), qissr(i)), qtot_in(i) - qcld(i))
     1498      dcfl_ini(i) = MIN(dcfl_ini(i), issrfra(i))
     1499      dqtl_ini(i) = MIN(dqtl_ini(i), qissr(i))
    15271500
    15281501      !--Add tendencies
     1502      cldfra(i)  = cldfra(i) + dcfl_ini(i)
     1503      qcld(i)    = qcld(i) + dqtl_ini(i)
     1504      qvc(i)     = qvc(i) + dcfl_ini(i) * qsat(i)
    15291505      issrfra(i) = issrfra(i) - dcfl_ini(i)
    15301506      qissr(i)   = qissr(i) - dqtl_ini(i)
    … …  
    15321508      qclr(i)    = qclr(i) - dqtl_ini(i)
    15331509
    1534       cldfra(i)  = cldfra(i) + dcfl_ini(i)
    1535       qcld(i)    = qcld(i) + dqtl_ini(i)
    1536       qvc(i)     = qvc(i) + dcfl_ini(i) * qsat(i)
    15371510      lincontfra(i) = lincontfra(i) + dcfl_cir(i) + dcfl_ini(i)
    15381511      qlincont(i)   = qlincont(i) + dqtl_cir(i) + dqtl_ini(i)
    15391512      circontfra(i) = circontfra(i) - dcfl_cir(i)
    15401513      qcircont(i)   = qcircont(i) - dqtl_cir(i)
     1514
     1515
     1516      !-------------------------------------------
     1517      !--       FINAL BARRIERS AND OUTPUTS      --
     1518      !-------------------------------------------
     1519
     1520      IF ( (lincontfra(i) .LT. eps) .OR. (qlincont(i) .LT. (qsat(i) * lincontfra(i))) ) THEN
     1521        cldfra(i) = cldfra(i) - lincontfra(i)
     1522        qcld(i) = qcld(i) - qlincont(i)
     1523        qvc(i) = qvc(i) - qsat(i) * lincontfra(i)
     1524        lincontfra(i) = 0.
     1525        qlincont(i) = 0.
     1526      ENDIF
     1527
     1528      IF ( (circontfra(i) .LT. eps) .OR. (qcircont(i) .LT. (qsat(i) * circontfra(i))) ) THEN
     1529        cldfra(i) = cldfra(i) - circontfra(i)
     1530        qcld(i) = qcld(i) - qcircont(i)
     1531        qvc(i) = qvc(i) - qsat(i) * circontfra(i)
     1532        circontfra(i) = 0.
     1533        qcircont(i) = 0.
     1534      ENDIF
     1535     
     1536      IF ( cldfra(i) .LT. eps ) THEN
     1537        !--If the cloud is too small, it is sublimated.
     1538        cldfra(i) = 0.
     1539        qcld(i)   = 0.
     1540        qvc(i)    = 0.
     1541        lincontfra(i) = 0.
     1542        qlincont(i)   = 0.
     1543        circontfra(i) = 0.
     1544        qcircont(i)   = 0.
     1545        qincld(i) = qsat(i)
     1546      ELSE
     1547        qincld(i) = qcld(i) / cldfra(i)
     1548      ENDIF
     1549
     1550      cldfra(i) = MAX(cldfra(i), lincontfra(i) + circontfra(i))
     1551      qcld(i) = MAX(qcld(i), qlincont(i) + qcircont(i))
     1552      qvc(i) = MAX(qvc(i), qsat(i) * ( lincontfra(i) + circontfra(i) ))
    15411553
    15421554      !--Diagnostics
    … …  
    15591571      dqtc_mix(i) = dqtc_mix(i) / dtime
    15601572
    1561       !-------------------------------------------
    1562       !--       FINAL BARRIERS AND OUTPUTS      --
    1563       !-------------------------------------------
    1564 
    1565       cldfra(i) = MIN(cldfra(i), totfra_in(i))
    1566       qcld(i) = MIN(qcld(i), qtot_in(i))
    1567       qvc(i) = MIN(qvc(i), qcld(i))
    1568 
    1569       IF ( cldfra(i) .LT. eps ) THEN
    1570         !--If the cloud is too small, it is sublimated.
    1571         cldfra(i) = 0.
    1572         qcld(i)   = 0.
    1573         qvc(i)    = 0.
    1574         lincontfra(i) = 0.
    1575         qlincont(i)   = 0.
    1576         circontfra(i) = 0.
    1577         qcircont(i)   = 0.
    1578         qincld(i) = qsat(i)
    1579       ELSE
    1580         qincld(i) = qcld(i) / cldfra(i)
    1581       ENDIF ! cldfra .LT. eps
    1582 
    1583       IF ( (lincontfra(i) .LT. eps) .OR. (qlincont(i) .LT. (qsat(i) * lincontfra(i))) ) THEN
    1584         lincontfra(i) = 0.
    1585         qlincont(i) = 0.
    1586       ENDIF
    1587 
    1588       IF ( (circontfra(i) .LT. eps) .OR. (qcircont(i) .LT. (qsat(i) * circontfra(i))) ) THEN
    1589         circontfra(i) = 0.
    1590         qcircont(i) = 0.
    1591       ENDIF
    1592 
    15931573    ENDIF ! keepgoing
    15941574  ENDDO ! loop on klon

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

    r5643 r5684  
    687687      REAL, SAVE, ALLOCATABLE :: qice_lincont(:,:), qice_circont(:,:)
    688688      !$OMP THREADPRIVATE(qice_lincont, qice_circont)
     689      REAL, SAVE, ALLOCATABLE :: qradice_lincont(:,:), qradice_circont(:,:)
     690      !$OMP THREADPRIVATE(qradice_lincont, qradice_circont)
    689691      REAL, SAVE, ALLOCATABLE :: dcfl_cir(:,:), dqtl_cir(:,:)
    690692      !$OMP THREADPRIVATE(dcfl_cir, dqtl_cir)
    … …  
    12721274      ALLOCATE(potcontfraP(klon,klev), potcontfraNP(klon,klev))
    12731275      ALLOCATE(qice_lincont(klon,klev), qice_circont(klon,klev))
     1276      ALLOCATE(qradice_lincont(klon,klev), qradice_circont(klon,klev))
    12741277      ALLOCATE(dcfl_cir(klon,klev), dqtl_cir(klon,klev))
    12751278      ALLOCATE(dcfl_ini(klon,klev), dqil_ini(klon,klev), dqtl_ini(klon,klev))
    … …  
    16941697      DEALLOCATE(d_q_avi, flight_dist, flight_h2o)
    16951698      DEALLOCATE(Tcritcont, qcritcont, potcontfraP, potcontfraNP)
    1696       DEALLOCATE(qice_lincont, qice_circont, dcfl_cir, dqtl_cir, dcfl_ini, dqil_ini)
     1699      DEALLOCATE(qice_lincont, qice_circont, qradice_lincont, qradice_circont)
     1700      DEALLOCATE(dcfl_cir, dqtl_cir, dcfl_ini, dqil_ini)
    16971701      DEALLOCATE(dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, dcfl_mix, dqil_mix, dqtl_mix)
    16981702      DEALLOCATE(dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix)

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

    r5643 r5684  
    334334       !-- LSCP - aviation and contrails variables
    335335       cfl_seri, cfc_seri, qtl_seri, qtc_seri, d_cfl_dyn, d_cfc_dyn, d_qtl_dyn, d_qtc_dyn, &
    336        d_q_avi, flight_dist, flight_h2o, qice_lincont, qice_circont, &
     336       d_q_avi, flight_dist, flight_h2o, &
     337       qice_lincont, qice_circont, qradice_lincont, qradice_circont, &
    337338       Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    338339       cldfra_nocont, cldtau_nocont, cldemi_nocont, cldh_nocont, &
    … …  
    13951396       IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) &
    13961397          CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
    1397        CALL init_read_aviation_emissions
     1398       IF ( ok_plane_h2o .OR. ok_plane_contrail ) CALL init_read_aviation_emissions
    13981399
    13991400       print*, '================================================='
    … …  
    39883989         dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
    39893990         dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
    3990          cfl_seri, cfc_seri, qtl_seri, qtc_seri, flight_dist, flight_h2o, &
    3991          qice_lincont, qice_circont, Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
     3991         cfl_seri, cfc_seri, qtl_seri, qtc_seri, qice_lincont, qice_circont, &
     3992         flight_dist, flight_h2o, qradice_lincont, qradice_circont, &
     3993         Tcritcont, qcritcont, potcontfraP, potcontfraNP, &
    39923994         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
    39933995         qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
    … …  
    45414543               zfice, dNovrN, ptconv, rnebcon, clwcon, &
    45424544               !--AB contrails
    4543                cfl_seri, cfc_seri, qice_lincont, qice_circont, cldfra_nocont, &
     4545               cfl_seri, cfc_seri, qradice_lincont, qradice_circont, cldfra_nocont, &
    45444546               cldtau_nocont, cldemi_nocont, conttau, contemi, cldh_nocont, contcov, &
    45454547               fiwp_nocont, fiwc_nocont, ref_ice_nocont)
Note: See TracChangeset for help on using the changeset viewer.