Changeset 4913


Ignore:
Timestamp:
Apr 21, 2024, 1:29:21 PM (7 months ago)
Author:
evignon
Message:

amelioration de poprecip avec:

  • formule du rayon des cristaux
  • qvap du ciel clair pour l'évap des precip

Lea et Audran

Location:
LMDZ6/trunk/libf/phylmd
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90

    r4910 r4913  
    267267  REAL, DIMENSION(klon) :: zlh_solid
    268268  REAL, DIMENSION(klon) :: ztupnew
     269  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
    269270  REAL :: zm_solid         ! for liquid -> solid conversion
    270271  REAL, DIMENSION(klon) :: zrflclr, zrflcld
    … …  
    387388dqsfreez(:,:) = 0.
    388389dqsmelt(:,:)  = 0.
     390zqupnew(:)    = 0.
     391zqvapclr(:)   = 0.
    389392
    390393
    … …  
    413416        zq(i)=qt(i,k)
    414417        IF (.not. iftop) THEN
    415            ztupnew(i)=temp(i,k+1)+d_t(i,k+1)
     418           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
     419           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
     420           !--zqs(i) is the saturation specific humidity in the layer above
     421           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
    416422        ENDIF
    417423        !c_iso init of iso
    … …  
    424430            CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
    425431                              zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
     432                              zqvapclr, zqupnew, &
    426433                              zrfl, zrflclr, zrflcld, &
    427434                              zifl, ziflclr, ziflcld, &

  • LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90

    r4910 r4913  
    146146  !$OMP THREADPRIVATE(ok_poprecip)
    147147
     148  LOGICAL, SAVE, PROTECTED :: ok_corr_vap_evasub=.FALSE.    ! use the corrected version of clear-sky water vapor for the evap / subl processes
     149  !$OMP THREADPRIVATE(ok_corr_vap_evasub)
     150
    148151  REAL, SAVE, PROTECTED :: cld_lc_lsc_snow=2.e-5            ! snow autoconversion coefficient, stratiform. default from  Chaboureau and PInty 2006
    149152  !$OMP THREADPRIVATE(cld_lc_lsc_snow)
    … …  
    178181  REAL, SAVE, PROTECTED :: r_snow=1.E-3                    ! A COMMENTER TODO [m]
    179182  !$OMP THREADPRIVATE(r_snow)
    180 
    181   REAL, SAVE, PROTECTED :: r_ice=50.E-6                    ! A COMMENTER TODO [m]
    182   !$OMP THREADPRIVATE(r_ice)
    183183
    184184  REAL, SAVE, PROTECTED :: tau_auto_snow_min=100.          ! A COMMENTER TODO [s]
    … …  
    301301    CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp)
    302302    CALL getin_p('ok_poprecip',ok_poprecip)
     303    CALL getin_p('ok_corr_vap_evasub',ok_corr_vap_evasub)
    303304    CALL getin_p('rain_int_min',rain_int_min)
    304305    CALL getin_p('gamma_agg',gamma_agg)
    … …  
    358359    WRITE(lunout,*) 'lscp_ini, fl_cor_ebil', fl_cor_ebil
    359360    WRITE(lunout,*) 'lscp_ini, ok_poprecip', ok_poprecip
     361    WRITE(lunout,*) 'lscp_ini, ok_corr_vap_evasub', ok_corr_vap_evasub
    360362    WRITE(lunout,*) 'lscp_ini, rain_int_min:', rain_int_min
    361363    WRITE(lunout,*) 'lscp_ini, gamma_agg:', gamma_agg

  • LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90

    r4898 r4913  
    1818SUBROUTINE poprecip_precld( &
    1919           klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, &
    20            qprecip, precipfracclr, precipfraccld, &
     20           qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, &
    2121           rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub &
    2222           )
    … …  
    2525USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac
    2626USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
     27USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub
    2728USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf
    2829
    … …  
    4748REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfracclr  !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-]
    4849REAL,    INTENT(INOUT), DIMENSION(klon) :: precipfraccld  !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-]
     50
     51REAL,    INTENT(IN),    DIMENSION(klon) :: qvapclrup      !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg]
     52REAL,    INTENT(IN),    DIMENSION(klon) :: qtotupnew      !--total specific humidity IN THE LAYER ABOVE [kg/kg]
    4953
    5054REAL,    INTENT(INOUT), DIMENSION(klon) :: rain           !--flux of rain gridbox-mean coming from the layer above [kg/s/m2]
    … …  
    5963
    6064
    61 
    62 
    6365!--Integer for interating over klon
    6466INTEGER :: i
    … …  
    7072!--Saturation values
    7173REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati
     74!--Vapor in the clear sky
     75REAL :: qvapclr
    7276!--Fluxes tendencies because of evaporation and sublimation
    7377REAL :: dprecip_evasub_max, draineva, dsnowsub, dprecip_evasub_tot
    … …  
    129133ENDIF
    130134
     135! TODO Probleme : on utilise qvap total dans la maille pour l'evap / sub
     136! alors qu'on n'evap / sub que dans le ciel clair
     137! deux options pour cette routine :
     138! - soit on diagnostique le nuage AVANT l'evap / sub et on estime donc
     139! la fraction precipitante ciel clair dans la maille, ce qui permet de travailler
     140! avec des fractions, des fluxs et surtout un qvap dans le ciel clair
     141! - soit on pousse la param de Ludo au bout, et on prend un qvap de k+1
     142! dans le ciel clair, avec un truc comme :
     143!   qvapclr(k) = qvapclr(k+1)/qtot(k+1) * qtot(k)
     144! UPDATE : on code la seconde version. A voir si on veut mettre la premiere version.
     145
    131146
    132147DO i = 1, klon
    133148
    134149  !--If there is precipitation from the layer above
    135   IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN
     150  ! NOTE TODO here we could replace the condition on precipfracclr(i) by a condition
     151  ! such as eps or thresh_precip_frac, to remove the senseless barrier in the formulas
     152  ! of evap / sublim
     153  IF ( ( ( rain(i) + snow(i) ) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN
     154
     155    IF ( ok_corr_vap_evasub ) THEN
     156      !--Corrected version - we use the same water ratio between
     157      !--the clear and the cloudy sky as in the layer above. This
     158      !--extends the assumption that the cloud fraction is the same
     159      !--as the layer above. This is assumed only for the evap / subl
     160      !--process
     161      !--Note that qvap(i) is the total water in the gridbox, and
     162      !--precipfraccld(i) is the cloud fraction in the layer above
     163      qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) )
     164    ELSE
     165      !--Legacy version from Ludo - we use the total specific humidity
     166      !--for the evap / subl process
     167      qvapclr = qvap(i)
     168    ENDIF
    136169
    137170    !--Evaporation of liquid precipitation coming from above
    … …  
    142175    !--which does not need a barrier on rainclr, because included in the formula
    143176    draineva = precipfracclr(i) * ( MAX(0., &
    144              - coef_eva * ( 1. - expo_eva ) * (1. - qvap(i) / qsatl(i)) * dz(i) &
     177             - coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) &
    145178             + ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_eva ) &
    146179               ) )**( 1. / ( 1. - expo_eva ) ) - rainclr(i)
    147     ! TODO veut on vraiment mettre qvap ici, qui est en fait qtot dans la maille ?
    148180             
    149181    !--Evaporation is limited by 0
    … …  
    156188    !--which does not need a barrier on snowclr, because included in the formula
    157189    dsnowsub = precipfracclr(i) * ( MAX(0., &
    158              - coef_sub * ( 1. - expo_sub ) * (1. - qvap(i) / qsati(i)) * dz(i) &
     190             - coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) &
    159191             + ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_sub ) &
    160192             ) )**( 1. / ( 1. - expo_sub ) ) - snowclr(i)
    161     ! TODO veut on vraiment mettre qvap ici, qui est en fait qtot dans la maille ?
    162193
    163194    !--Sublimation is limited by 0
    … …  
    172203    !--It is expressed as a max flux dprecip_evasub_max
    173204   
    174     ! TODO veut on vraiment mettre qvap ici, qui est en fait qtot dans la maille ?
    175     dprecip_evasub_max = MIN(0., ( qvap(i) - qsat(i) ) * precipfracclr(i)) &
     205    dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr(i)) &
    176206                     * dhum_to_dflux(i)
    177207    dprecip_evasub_tot = draineva + dsnowsub
    … …  
    256286                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
    257287                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
    258                           rho_rain, r_rain, r_snow, rho_ice, r_ice,            &
     288                          rho_rain, r_rain, r_snow, rho_ice,                   &
    259289                          tau_auto_snow_min, tau_auto_snow_max,                &
    260290                          thresh_precip_frac, eps,                             &
    … …  
    344374REAL :: nb_snowflake_clr, nb_snowflake_cld
    345375REAL :: capa_snowflake, temp_wetbulb
     376REAL :: rho, r_ice
    346377REAL :: dqsclrmelt, dqscldmelt, dqstotmelt
    347378REAL, DIMENSION(klon) :: qzero, qsat, dqsat
    … …  
    685716    dqscldmelt = 0.
    686717
     718    !--We assume that the snowflakes are spherical
    687719    capa_snowflake = r_snow
    688     ! ATTENTION POUR L'INSTANT C'EST UN WBULB SELON FORMULE ECMWF
    689     ! ATTENTION EST-CE QVAP ?????
    690     ! TODO veut on vraiment mettre qvap ici, on pourrait vouloir qvap dans le nuage / clr sky ?
    691     ! (a mettre dans les deux IF)
    692     temp_wetbulb = temp(i) - ( qsat(i) - qvap(i) ) &
    693                  * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
    694                  - 40.637 * ( temp(i) - 275. ) )
    695    
    696720    !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971)
    697721    air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184   
    698 
    699722    !--rho_snow formula follows Brandes et al. 2007 (JAMC)
    700723    rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922)
    … …  
    702725    !--In clear air
    703726    IF ( ( snowclr(i) .GT. 0. ) .AND.  ( precipfracclr(i) .GT. 0. ) ) THEN
     727      !--Formula for the wet-bulb temperature from ECMWF (IFS)
     728      !--The vapor used is the vapor in the clear sky
     729      temp_wetbulb = temp(i) &
     730                   - ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) &
     731                   * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) &
     732                   - 40.637 * ( temp(i) - 275. ) )
    704733      !--Calculated according to
    705734      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
    … …  
    717746    !--In cloudy air
    718747    IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
     748      !--As the air is saturated, the wet-bulb temperature is equal to the
     749      !--temperature
     750      temp_wetbulb = temp(i)
    719751      !--Calculated according to
    720752      !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow
    … …  
    798830    dqrtotfreez_step1 = 0.
    799831
    800     IF ( ( qice(i) .GT. 0. ) .AND. ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
     832    IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. &
     833         ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN
    801834      dqrclrfreez = 0.
    802835      dqrcldfreez = 0.
    … …  
    822855      !--The assumption above corresponds to dNi = dNr, i.e.,
    823856      !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3)
     857      !--Dry density [kg/m3]
     858      rho = pplay(i) / temp(i) / RD
     859      !--r_ice formula from Sun and Rikus (1999)
     860      r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 &
     861            + 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2.
    824862      dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. )
    825863      dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i))
Note: See TracChangeset for help on using the changeset viewer.