Ignore:
Timestamp:
Mar 19, 2025, 3:53:17 PM (3 months ago)
Author:
aborella
Message:

Added resuspension of ice crystals in poprecip + contrails do not precipitate anymore + added different options to control the behavior of contrails

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

Legend:

Unmodified
Added
Removed

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

    r5575 r5579  
    7373SUBROUTINE contrails_formation_evolution( &
    7474      dtime, pplay, temp, qsat, qsatl, gamma_cond, rcont_seri, flight_dist, &
    75       cldfra, qvc, dz, pdf_loc, pdf_scale, pdf_alpha, &
     75      cldfra, qvc, pdf_loc, pdf_scale, pdf_alpha, &
    7676      Tcritcont, qcritcont, potcontfraP, potcontfraNP, contfra, &
    7777      dcontfra_cir, dcf_avi, dqvc_avi, dqi_avi &
    … …  
    100100REAL, INTENT(IN)  :: cldfra       ! cloud fraction [-]
    101101REAL, INTENT(IN)  :: qvc          ! gridbox-mean vapor in the cloud [kg/kg]
    102 REAL, INTENT(IN)  :: dz           ! cell width [m]
    103102REAL, INTENT(IN)  :: pdf_loc      ! location parameter of the clear sky PDF [%]
    104103REAL, INTENT(IN)  :: pdf_scale    ! scale parameter of the clear sky PDF [%]
    … …  
    207206dqvc_avi = 0.
    208207IF ( potcontfraP .GT. eps ) THEN
    209   contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA(dz)
     208  contrail_cross_section = CONTRAIL_CROSS_SECTION_ONERA()
    210209  contfra_new = MIN(1., flight_dist * dtime * contrail_cross_section)
    211210  dcf_avi = potcontfraP * contfra_new
    … …  
    610609
    611610!**********************************************************************************
    612 FUNCTION contrail_cross_section_onera(dz)
    613 
    614 USE lmdz_lscp_ini, ONLY: initial_width_contrails
     611FUNCTION contrail_cross_section_onera()
     612
     613USE lmdz_lscp_ini, ONLY: initial_width_contrails, initial_height_contrails
    615614
    616615IMPLICIT NONE
    617616
    618617!
    619 ! Input
    620 !
    621 REAL :: dz ! cell width [m]
    622 !
    623618! Output
    624619!
    … …  
    628623!
    629624
    630 contrail_cross_section_onera = initial_width_contrails * dz
     625contrail_cross_section_onera = initial_width_contrails * initial_height_contrails
    631626
    632627END FUNCTION contrail_cross_section_onera
    … …  
    642637    USE lmdz_xios
    643638    USE print_control_mod, ONLY: lunout
     639    USE lmdz_lscp_ini, ONLY: EI_H2O_aviation
    644640    IMPLICIT NONE
    645641
    … …  
    649645    ! Local variable
    650646    !----------------------------------------------------
    651     REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:)
     647    REAL, ALLOCATABLE :: flight_dist_mpi(:,:,:), flight_h2o_mpi(:,:,:)
    652648    INTEGER :: ierr
    653649
    … …  
    669665        ALLOCATE(flight_dist_mpi(klon_mpi, nleva,1), STAT=ierr)
    670666        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_dist_mpi',1)
     667        ALLOCATE(flight_h2o_mpi(klon_mpi, nleva,1), STAT=ierr)
     668        IF (ierr /= 0) CALL abort_physic('read_aviation_emissions', 'problem to allocate flight_h2o_mpi',1)
    671669        CALL xios_recv_field("KMFLOWN_interp", flight_dist_mpi(:,:,1))
     670        !CALL xios_recv_field("KGH2O_interp", flight_h2o_mpi(:,:,1))
     671        flight_h2o_mpi(:,:,:) = 0.
    672672        ! Get number of vertical levels and level values
    673673        CALL xios_get_axis_attr( "aviation_lev", value=aviation_lev(:))
    … …  
    677677    ! (klon_mpi,klon) = (200,50) avec 80 MPI, 4 OMP, nbp40
    678678    CALL scatter_omp(flight_dist_mpi, flight_dist_read)
     679    CALL scatter_omp(flight_h2o_mpi, flight_h2o_read)
    679680    CALL bcast_omp(aviation_lev)
    680681
    … …  
    687688    ! flight_dist_read(klon,nleva) -> flight_dist(klon, klev)
    688689    ! flight_h2o_read(klon,nleva) -> flight_h2o(klon, klev)
    689     USE print_control_mod, ONLY: lunout
     690
    690691    USE lmdz_lscp_ini, ONLY: RD, RG
     692    USE lmdz_lscp_ini, ONLY: aviation_coef
    691693
    692694    IMPLICIT NONE
    … …  
    760762            flight_dist(i,k) = flight_dist(i,k) / dz
    761763            flight_h2o(i,k) = flight_h2o(i,k) / dz
     764           
     765            !--Enhancement factor
     766            flight_dist(i,k) = flight_dist(i,k) * aviation_coef
     767            flight_h2o(i,k) = flight_h2o(i,k) * aviation_coef
    762768        ENDDO
    763769    ENDDO

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

    r5575 r5579  
    3131     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
    3232     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
    33      dqsmelt, dqsfreez)
     33     dqsmelt, dqsfreez, dcfres, dqsres, dqvcres)
    3434
    3535!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    … …  
    264264  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
    265265  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
     266  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcfres         !--cloud fraction tendency due to resuspension of ice crystals [kg/kg/s]
     267  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsres         !--snow tendency due to resuspension of ice crystals [kg/kg/s]
     268  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvcres        !--cloud water vapor tendency due to resuspension of ice crystals [kg/kg/s]
    266269
    267270  ! for thermals
    … …  
    275278  ! LOCAL VARIABLES:
    276279  !----------------
    277   REAL,DIMENSION(klon) :: qice_ini
     280  REAL, DIMENSION(klon) :: cldfra_in, qvc_in, qliq_in, qice_in
    278281  REAL, DIMENSION(klon,klev) :: ctot
    279282  REAL, DIMENSION(klon,klev) :: ctot_vol
    … …  
    319322  REAL :: delta_z
    320323  ! for contrails
     324  REAL, DIMENSION(klon) :: zqcont
    321325  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrail)
    322326  ! Constants used for calculating ratios that are advected (using a parent-child
    … …  
    483487DO k = klev, 1, -1
    484488
    485     qice_ini = qi_seri(:,k)
     489    cldfra_in(:) = cf_seri(:,k)
     490    qvc_in(:) = rvc_seri(:,k) * q_seri(:,k)
     491    qliq_in(:) = ql_seri(:,k)
     492    qice_in(:) = qi_seri(:,k)
    486493
    487494    IF (k.LE.klev-1) THEN
    … …  
    522529                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
    523530                        zqvapclr, zqupnew, &
    524                         cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
     531                        cldfra_in, qvc_in, qliq_in, qice_in, &
    525532                        zrfl, zrflclr, zrflcld, &
    526533                        zifl, ziflclr, ziflcld, &
    527                         dqreva(:,k), dqssub(:,k) &
     534                        dqreva(:,k), dqssub(:,k), &
     535                        dqsres(:,k), dcfres(:,k), dqvcres(:,k) &
    528536                        )
    529537
    … …  
    746754                        klon, dtime, missing_val, &
    747755                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
    748                         cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
     756                        cldfra_in, qvc_in, qliq_in, qice_in, &
    749757                        shear, tke_dissip(:,k), cell_area, stratomask(:,k), &
    750758                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
    … …  
    771779                  IF (iflag_icefrac .GE. 2) THEN
    772780                     ! phase partitioning depending on temperature. activates here in the iteration process if iflag_icefrac > 2
    773                      CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
     781                     CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_in, ziflcld, zqn, &
    774782                     rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
    775783                  ELSE
    … …  
    858866        ! Partition function depending on tke for non shallow-convective clouds, erase previous estimation
    859867        IF (iflag_icefrac .GE. 1) THEN
    860            CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
     868           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_in, ziflcld, zqn, &
    861869           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
    862870        ENDIF
    … …  
    947955    ENDDO
    948956
     957    IF (ok_plane_contrail) THEN
     958      !--Contrails do not precipitate. We remove then from the variables temporarily
     959      IF (rneb(i,k) .GT. eps) THEN
     960        zqcont(i) = zcond(i) * zfice(i) * contfra(i,k) / rneb(i,k)
     961      ELSE
     962        zqcont(i) = 0.
     963      ENDIF
     964      rneb(i,k) = rneb(i,k) - contfra(i,k)
     965      zoliqi(i) = zoliqi(i) - zqcont(i)
     966    ENDIF
     967
    949968    !================================================================
    950969    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
    … …  
    979998
    980999    ENDIF ! ok_poprecip
     1000   
     1001    IF (ok_plane_contrail) THEN
     1002      !--Contrails are reintroduced in the variables
     1003      rneb(i,k) = rneb(i,k) + contfra(i,k)
     1004      zoliqi(i) = zoliqi(i) + zqcont(i)
     1005    ENDIF
    9811006
    9821007    ! End of precipitation processes after cloud formation

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

    r5575 r5579  
    9595SUBROUTINE condensation_ice_supersat( &
    9696      klon, dtime, missing_val, pplay, paprsdn, paprsup, &
    97       cf_seri, rvc_seri, ql_seri, qi_seri, shear, pbl_eps, cell_area, stratomask, &
     97      cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, stratomask, &
    9898      temp, qtot, qsat, gamma_cond, ratqs, keepgoing, &
    9999      cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, &
    … …  
    141141REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsdn       ! pressure at the lower interface [Pa]
    142142REAL,     INTENT(IN)   , DIMENSION(klon) :: paprsup       ! pressure at the upper interface [Pa]
    143 REAL,     INTENT(IN)   , DIMENSION(klon) :: cf_seri       ! cloud fraction [-]
    144 REAL,     INTENT(IN)   , DIMENSION(klon) :: rvc_seri      ! gridbox-mean water vapor in cloud [kg/kg]
    145 REAL,     INTENT(IN)   , DIMENSION(klon) :: ql_seri       ! specific liquid water content [kg/kg]
    146 REAL,     INTENT(IN)   , DIMENSION(klon) :: qi_seri       ! specific ice water content [kg/kg]
     143REAL,     INTENT(IN)   , DIMENSION(klon) :: cldfra_in     ! cloud fraction [-]
     144REAL,     INTENT(IN)   , DIMENSION(klon) :: qvc_in        ! gridbox-mean water vapor in cloud [kg/kg]
     145REAL,     INTENT(IN)   , DIMENSION(klon) :: qliq_in       ! specific liquid water content [kg/kg]
     146REAL,     INTENT(IN)   , DIMENSION(klon) :: qice_in       ! specific ice water content [kg/kg]
    147147REAL,     INTENT(IN)   , DIMENSION(klon) :: shear         ! vertical shear [s-1]
    148148REAL,     INTENT(IN)   , DIMENSION(klon) :: pbl_eps       ! TKE dissipation [m2/s3]
    … …  
    304304        !--and the condensation process is slightly adapted
    305305        !--This can happen only if ok_weibull_warm_clouds = .TRUE.
    306         cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
    307         qcld(i)   = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i)))
    308         qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
     306        cldfra(i) = MAX(0., MIN(1., cldfra_in(i)))
     307        qcld(i)   = MAX(0., MIN(qtot(i), qliq_in(i) + qice_in(i) + qvc_in(i)))
     308        qvc(i)    = MAX(0., MIN(qcld(i), qvc_in))
    309309        ok_warm_cloud = .TRUE.
    310310      ELSE
    … …  
    312312        !--are consistent. In some rare cases, i.e. the cloud water vapor
    313313        !--can be greater than the total water in the gridbox
    314         cldfra(i) = MAX(0., MIN(1., cf_seri(i)))
    315         qcld(i)   = MAX(0., MIN(qtot(i), rvc_seri(i) * qtot(i) + qi_seri(i)))
    316         qvc(i)    = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i)))
     314        cldfra(i) = MAX(0., MIN(1., cldfra_in(i)))
     315        qcld(i)   = MAX(0., MIN(qtot(i), qvc_in(i) + qice_in(i)))
     316        qvc(i)    = MAX(0., MIN(qcld(i), qvc_in(i)))
    317317        ok_warm_cloud = .FALSE.
    318318      ENDIF
    … …  
    788788            dtime, pplay(i), temp(i), qsat(i), qsatl(i), gamma_cond(i), &
    789789            rcont_seri(i), flight_dist(i), cldfra(i), qvc(i), &
    790             dz, pdf_loc, pdf_scale, pdf_alpha, &
     790            pdf_loc, pdf_scale, pdf_alpha, &
    791791            Tcritcont(i), qcritcont(i), potcontfraP(i), potcontfraNP(i), contfra(i), &
    792792            dcontfra_cir(i), dcf_avi(i), dqvc_avi(i), dqi_avi(i) &

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

    r5577 r5579  
    245245  REAL, SAVE, PROTECTED :: initial_width_contrails=200.      ! [m] initial width of the linear contrails formed
    246246  !$OMP THREADPRIVATE(initial_width_contrails)
     247
     248  REAL, SAVE, PROTECTED :: initial_height_contrails=200.     ! [m] initial height of the linear contrails formed
     249  !$OMP THREADPRIVATE(initial_height_contrails)
     250
     251  REAL, SAVE, PROTECTED :: aviation_coef=1.                  ! [-] scaling factor for aviation emissions and flown distance
     252  !$OMP THREADPRIVATE(aviation_coef)
    247253  !--End of the parameters for aviation
    248254
    … …  
    337343  REAL, SAVE, PROTECTED :: snow_fallspeed_cld               ! Snow fall velocity in cloudy sky [m/s]
    338344  !$OMP THREADPRIVATE(snow_fallspeed_cld)
     345
     346  INTEGER, SAVE, PROTECTED :: ok_precip_resuspension=.FALSE.! Flag to activate the resuspension of ice crystals
     347  !$OMP THREADPRIVATE(ok_precip_resuspension)
     348
     349  REAL, SAVE, PROTECTED :: snow_thresh_resuspension=1.e-5   ! [kg/m2/s] Threshold flux below which part of the snow flux will be resuspended
     350  !$OMP THREADPRIVATE(snow_thresh_resuspension)
    339351  !--End of the parameters for poprecip
    340352
    … …  
    455467    CALL getin_p('snow_fallspeed_clr',snow_fallspeed_clr)
    456468    CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld)
     469    CALL getin_p('ok_precip_resuspension',ok_precip_resuspension)
     470    CALL getin_p('snow_thresh_resuspension',snow_thresh_resuspension)
    457471    ! for condensation and ice supersaturation
    458472    CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds)
    … …  
    488502    CALL getin_p('linear_contrails_lifetime',linear_contrails_lifetime)
    489503    CALL getin_p('initial_width_contrails',initial_width_contrails)
     504    CALL getin_p('initial_height_contrails',initial_height_contrails)
     505    CALL getin_p('aviation_coef',aviation_coef)
    490506
    491507
    … …  
    550566    WRITE(lunout,*) 'lscp_ini, snow_fallspeed_clr:', snow_fallspeed_clr
    551567    WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
     568    WRITE(lunout,*) 'lscp_ini, ok_precip_resuspension:', ok_precip_resuspension
     569    WRITE(lunout,*) 'lscp_ini, snow_thresh_resuspension:', snow_thresh_resuspension
    552570    ! for condensation and ice supersaturation
    553571    WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
    … …  
    581599    WRITE(lunout,*) 'lscp_ini, linear_contrails_lifetime:', linear_contrails_lifetime
    582600    WRITE(lunout,*) 'lscp_ini, initial_width_contrails:', initial_width_contrails
     601    WRITE(lunout,*) 'lscp_ini, initial_height_contrails:', initial_height_contrails
     602    WRITE(lunout,*) 'lscp_ini, aviation_coef:', aviation_coef
    583603
    584604

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

    r5431 r5579  
    714714           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
    715715           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
    716            cldfra, rvc_seri, qliq, qice, &
    717            rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
     716           cldfra, qvc, qliq, qice, &
     717           rain, rainclr, raincld, snow, snowclr, snowcld, &
     718           dqreva, dqssub, dqsres, dcfres, dqvcres &
    718719           )
    719720
    720 USE lmdz_lscp_ini, ONLY : prt_level, lunout
    721721USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
    722722USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
    723723USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds
    724724USE lmdz_lscp_ini, ONLY : eps, temp_nowater
     725USE lmdz_lscp_ini, ONLY : ok_precip_resuspension, snow_thresh_resuspension
    725726USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
    726727
    … …  
    749750REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
    750751
    751 REAL,    INTENT(IN),    DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
    752 REAL,    INTENT(IN),    DIMENSION(klon) :: rvc_seri       !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg]
    753 REAL,    INTENT(IN),    DIMENSION(klon) :: qliq           !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
    754 REAL,    INTENT(IN),    DIMENSION(klon) :: qice           !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
     752REAL,    INTENT(INOUT), DIMENSION(klon) :: cldfra         !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-]
     753REAL,    INTENT(INOUT), DIMENSION(klon) :: qvc            !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg]
     754REAL,    INTENT(INOUT), DIMENSION(klon) :: qliq           !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
     755REAL,    INTENT(INOUT), DIMENSION(klon) :: qice           !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg]
    755756
    756757
    … …  
    764765REAL,    INTENT(OUT),   DIMENSION(klon) :: dqreva         !--rain tendency due to evaporation [kg/kg/s]
    765766REAL,    INTENT(OUT),   DIMENSION(klon) :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
     767REAL,    INTENT(OUT),   DIMENSION(klon) :: dqsres         !--snow tendency due to resuspension of ice crystals [kg/kg/s]
     768REAL,    INTENT(OUT),   DIMENSION(klon) :: dcfres         !--cloud fraction tendency due to resuspension of ice crystals [kg/kg/s]
     769REAL,    INTENT(OUT),   DIMENSION(klon) :: dqvcres        !--cloud water vapor tendency due to resuspension of ice crystals [kg/kg/s]
    766770
    767771
    … …  
    787791!--Specific heat constant
    788792REAL :: cpair, cpw
     793!--For resuspension of ice crystals
     794REAL :: dsnowres, qiceinprecip
    789795
    790796!--Initialisation
    … …  
    794800dqrevap   = 0.
    795801dqssubl   = 0.
     802IF ( ok_ice_supersat .AND. ok_precip_resuspension ) THEN
     803  dqsres(:) = 0.
     804  dcfres(:) = 0.
     805  dqvcres(:) = 0.
     806ENDIF
    796807
    797808!-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt
    … …  
    816827
    817828IF (iftop) THEN
    818 
     829 
    819830  DO i = 1, klon
    820831    qprecip(i) = 0.
    … …  
    894905      !--water vapor in the cloud and cloud fraction
    895906      IF ( ok_unadjusted_clouds .AND. ( cldfra(i) .GT. eps ) ) THEN
    896         qvapcld = rvc_seri(i) * qvap(i) / cldfra(i)
     907        qvapcld = qvc(i) / cldfra(i)
    897908      ENDIF
    898909      !--We can diagnose completely the water vapor in clear sky, because all
    … …  
    901912      !--Note that qvap(i) is the total water in the gridbox
    902913      IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN
    903         qvapclr = ( qvap(i) - qice(i) - qliq(i) - rvc_seri(i) * qvap(i) ) / ( 1. - cldfra(i) )
     914        qvapclr = ( qvap(i) - qice(i) - qliq(i) - qvc(i) ) / ( 1. - cldfra(i) )
    904915      ENDIF
    905916    ELSEIF ( ok_corr_vap_evasub ) THEN
    … …  
    10701081    raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva)
    10711082    snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub)
     1083   
     1084
     1085    !--If ok_precip_to_clouds is set to .TRUE., some of the precipitations can go back to clouds
     1086    !--Only the snow flux can go back to cloud, and only if there is no rain flux
     1087    !--(mixed-phase clouds)
     1088    !--The ok_ice_supersaturation is used to ensure that cloud fraction is advected
     1089    !--but maybe it is not needed
     1090    IF ( ok_ice_supersat .AND. ok_precip_resuspension ) THEN
     1091
     1092      IF ( (precipfracclr_tmp(i) .GT. 0.) .AND. (rainclr_tmp(i) .LT. eps) ) THEN
     1093        dsnowres = snowclr_tmp(i) * EXP( - snowclr_tmp(i) / precipfracclr_tmp(i) / snow_thresh_resuspension )
     1094        dqsres(i) = dqsres(i) - dsnowres / dhum_to_dflux(i)
     1095        !--The following line determines the in-cloud ice water content of the newly formed cloud
     1096        !--It is a linear combination between the in-precip ice water content (left term) and the
     1097        !--in-existing-cloud ice water content (right term). This is done so that the existing cloud
     1098        !--is not too modified, i.e. unphysical ice fluxes going from the existing cloud to the new
     1099        !--cloud. If the existing cloud is already large, the newly formed cloud fraction is reduced
     1100        !--to match the in-cloud ice water content of the existing cloud.
     1101        !--NB. this is done on qi, not on qvc. Not sure what impact that implies
     1102        qiceinprecip = - dqsres(i) / precipfracclr_tmp(i) * ( 1. - cldfra(i) ) + qice(i)
     1103        dcfres(i) = - dqsres(i) / qiceinprecip
     1104        dqvcres(i) = qvapclr * dcfres(i)
     1105        !--Add tendencies
     1106        snowcld_tmp(i) = snowcld_tmp(i) + dsnowres
     1107      ENDIF
     1108
     1109      IF ( (precipfraccld_tmp(i) .GT. 0.) .AND. (raincld_tmp(i) .LT. eps) ) THEN
     1110        dsnowres = snowcld_tmp(i) * EXP( - snowcld_tmp(i) / precipfraccld_tmp(i) / snow_thresh_resuspension )
     1111        dqsres(i) = dqsres(i) - dsnowres / dhum_to_dflux(i)
     1112        !--Add tendencies
     1113        snowcld_tmp(i) = snowcld_tmp(i) + dsnowres
     1114      ENDIF
     1115
     1116      !--Add tendencies
     1117      cldfra(i) = cldfra(i) + dcfres(i)
     1118      qvc(i) = qvc(i) + dqvcres(i)
     1119      qice(i) = qice(i) - dqsres(i)
     1120    ENDIF
    10721121
    10731122

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

    r5575 r5579  
    717717      REAL, SAVE, ALLOCATABLE :: dqsfreez(:,:)
    718718      !$OMP THREADPRIVATE(dqsfreez)
     719      REAL, SAVE, ALLOCATABLE :: dcfres(:,:)
     720      !$OMP THREADPRIVATE(dcfres)
     721      REAL, SAVE, ALLOCATABLE :: dqsres(:,:)
     722      !$OMP THREADPRIVATE(dqsres)
     723      REAL, SAVE, ALLOCATABLE :: dqvcres(:,:)
     724      !$OMP THREADPRIVATE(dqvcres)
    719725
    720726! variables for stratospheric aerosol
    … …  
    12381244      ALLOCATE(dqrauto(klon,klev), dqrcol(klon,klev), dqrmelt(klon,klev), dqrfreez(klon,klev))
    12391245      ALLOCATE(dqsauto(klon,klev), dqsagg(klon,klev), dqsrim(klon,klev), dqsmelt(klon,klev), dqsfreez(klon,klev))
     1246      ALLOCATE(dcfres(klon,klev), dqsres(klon,klev), dqvcres(klon,klev))
    12401247
    12411248IF (CPPKEY_STRATAER) THEN
    … …  
    16451652      DEALLOCATE(dqrauto, dqrcol, dqrmelt, dqrfreez)
    16461653      DEALLOCATE(dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
     1654      DEALLOCATE(dcfres, dqsres, dqvcres)
    16471655
    16481656IF (CPPKEY_STRATAER) THEN

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

    r5575 r5579  
    16111611  TYPE(ctrl_out), SAVE :: o_dqsfreez = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
    16121612    'dqsfreez', 'LS snow tendency due to freezing', 'kg/kg/s', (/ ('', i=1, 10) /))
     1613  TYPE(ctrl_out), SAVE :: o_dcfres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     1614    'dcfres', 'LS cloud fraction tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
     1615  TYPE(ctrl_out), SAVE :: o_dqsres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     1616    'dqsres', 'LS snow tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
     1617  TYPE(ctrl_out), SAVE :: o_dqvcres = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), &
     1618    'dqvcres', 'LS cloud water vapor tendency due to resuspension', 'kg/kg/s', (/ ('', i=1, 10) /))
    16131619  TYPE(ctrl_out), SAVE :: o_rhum = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), &
    16141620    'rhum', 'Relative humidity', '-', (/ ('', i=1, 10) /))

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

    r5575 r5579  
    145145         o_qrainlsc, o_qsnowlsc, o_dqreva, o_dqrauto, o_dqrcol, o_dqrmelt, o_dqrfreez, &
    146146         o_dqssub, o_dqsauto, o_dqsagg, o_dqsrim, o_dqsmelt, o_dqsfreez, &
     147         o_dcfres, o_dqsres, o_dqvcres, &
    147148         o_duphy, o_dtphy, o_dqphy, o_dqphy2d, o_dqlphy, o_dqlphy2d, &
    148149         o_dqsphy, o_dqsphy2d, o_dqbsphy, o_dqbsphy2d, o_albe_srf, o_z0m_srf, o_z0h_srf, &
    … …  
    268269         o_SAD_sulfate, o_reff_sulfate, o_sulfmmr, o_nd_mode, o_sulfmmr_mode
    269270
    270     USE lmdz_lscp_ini, ONLY: ok_poprecip
     271    USE lmdz_lscp_ini, ONLY: ok_poprecip, ok_precip_resuspension
    271272
    272273    USE phys_output_ctrlout_mod, ONLY: o_heat_volc, o_cool_volc !NL
    … …  
    383384         dqrauto,dqrcol,dqrmelt,dqrfreez, &
    384385         dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
     386         dcfres, dqsres, dqvcres, &
    385387         d_t_dyn,  &
    386388         d_q_dyn,  d_ql_dyn, d_qs_dyn, d_qbs_dyn,  &
    … …  
    20932095           CALL histwrite_phy(o_dqsfreez, dqsfreez)
    20942096           CALL histwrite_phy(o_dqsrim, dqsrim)
     2097           IF ( ok_precip_resuspension ) THEN
     2098            CALL histwrite_phy(o_dcfres, dcfres)
     2099            CALL histwrite_phy(o_dqsres, dqsres)
     2100            CALL histwrite_phy(o_dqvcres, dqvcres)
     2101           ENDIF
    20952102           ELSE
    20962103            CALL histwrite_phy(o_dqreva, dqreva)

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

    r5575 r5579  
    187187       ! proprecip
    188188       qraindiag, qsnowdiag, &
    189        dqreva, dqssub, &
     189       dqreva, dqssub, dcfres, dqsres, dqvcres, &
    190190       dqrauto,dqrcol,dqrmelt,dqrfreez, &
    191191       dqsauto,dqsagg,dqsrim,dqsmelt,dqsfreez, &
    … …  
    39233923         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
    39243924         qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
    3925          dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
     3925         dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez, &
     3926         dcfres, dqsres, dqvcres)
    39263927
    39273928
Note: See TracChangeset for help on using the changeset viewer.