Changeset 4832


Ignore:
Timestamp:
Feb 23, 2024, 8:07:46 PM (2 months ago)
Author:
evignon
Message:

modification du processus de freezing et debogage de poprecip
suite aux tests sur les cas 1D

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

Legend:

Unmodified
Added
Removed

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

    r4830 r4832  
    161161  !$OMP THREADPRIVATE(rho_rain)
    162162
     163  REAL, SAVE, PROTECTED :: rho_ice=920.                    ! A COMMENTER TODO [kg/m3]
     164  !$OMP THREADPRIVATE(rho_ice)
     165
    163166  REAL, SAVE, PROTECTED :: rho_snow                        ! A COMMENTER TODO [kg/m3]
    164167  !$OMP THREADPRIVATE(rho_snow)
    … …  
    169172  REAL, SAVE, PROTECTED :: r_snow=1.E-3                    ! A COMMENTER TODO [m]
    170173  !$OMP THREADPRIVATE(r_snow)
     174
     175  REAL, SAVE, PROTECTED :: r_ice=50.E-6                    ! A COMMENTER TODO [m]
     176  !$OMP THREADPRIVATE(r_ice)
    171177
    172178  REAL, SAVE, PROTECTED :: tau_auto_snow_min=100.          ! A COMMENTER TODO [s]

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

    r4830 r4832  
    136136    !--multiplying by dz = - dP / g / rho (line 3-4)
    137137    !--formula from Sundqvist 1988, Klemp & Wilhemson 1978
    138 
    139138    draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) &
    140139             * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva &
    … …  
    197196
    198197    !--Add tendencies
    199     rainclr(i) = rainclr(i) + draineva
    200     snowclr(i) = snowclr(i) + dsnowsub
     198
     199    rainclr(i) = MAX(0., rainclr(i) + draineva)
     200    snowclr(i) = MAX(0., snowclr(i) + dsnowsub)
    201201
    202202    !--If there is no more precip fluxes, the precipitation fraction in clear
    … …  
    248248                          cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, &
    249249                          thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, &
    250                           rho_rain, rho_snow, r_rain, r_snow,                  &
     250                          rho_rain, rho_snow, r_rain, r_snow, rho_ice, r_ice,  &
    251251                          tau_auto_snow_min, tau_auto_snow_max,                &
    252252                          thresh_precip_frac, eps, air_thermal_conduct,        &
    … …  
    340340REAL :: dqrfreez_max
    341341REAL :: tau_freez
    342 REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez
     342REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2
    343343REAL :: coef_freez
    344344REAL :: dqifreez        !--loss of ice cloud content due to collection of ice from rain [kg/kg/s]
    … …  
    371371  dqifreez= 0.
    372372
    373   !--hum_to_flux: coef to convert a specific quantity to a flux
     373  !--hum_to_flux: coef to convert a delta in specific quantity to a flux
    374374  !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt
    375375  hum_to_flux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime
    … …  
    398398  !--evaporated (see barrier at the end of poprecip_precld)
    399399  !--NB. precipfraccld(i) is here the cloud fraction of the layer above
    400   precipfractot = 1. - ( 1. - precipfractot ) * &
     400  !precipfractot = 1. - ( 1. - precipfractot ) * &
     401  !               ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
     402  !             / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
     403
     404
     405  IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN
     406    precipfractot = 1.
     407  ELSE
     408    precipfractot = 1. - ( 1. - precipfractot ) * &
    401409                 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) &
    402                / ( 1. - MIN( precipfraccld(i), 1. - eps ) )
    403 
     410               / ( 1. - precipfraccld(i) )
     411  ENDIF
    404412
    405413  !--precipfraccld(i) is here the cloud fraction of the layer above
    … …  
    451459  precipfraccld(i) = precipfraccld(i) + dprecipfraccld
    452460  precipfracclr(i) = precipfracclr(i) + dprecipfracclr
    453   rainclr(i) = rainclr(i) + drainclr
    454   snowclr(i) = snowclr(i) + dsnowclr
    455   raincld(i) = raincld(i) - drainclr
    456   snowcld(i) = snowcld(i) - dsnowclr
     461
     462  rainclr(i) = MAX( 0. , rainclr(i) + drainclr )
     463  snowclr(i) = MAX( 0. , snowclr(i) + dsnowclr )
     464  raincld(i) = MAX( 0. , raincld(i) - drainclr )
     465  snowcld(i) = MAX( 0. , snowcld(i) - dsnowclr )
    457466 
    458 
    459467  !--If vertical heterogeneity is taken into account, we use
    460468  !--the "true" volume fraction instead of a modified
    … …  
    513521      !--Add tendencies
    514522      qliq(i) = qliq(i) + dqlcol
    515       raincld(i) = raincld(i) - dqlcol * hum_to_flux(i)
     523      raincld(i) =  MAX( 0. , raincld(i) - dqlcol * hum_to_flux(i) )
    516524
    517525      !--Diagnostic tendencies
    … …  
    535543      !--Add tendencies
    536544      qice(i) = qice(i) + dqiagg
    537       snowcld(i) = snowcld(i) - dqiagg * hum_to_flux(i)
     545      snowcld(i) = MAX( 0. , snowcld(i) - dqiagg * hum_to_flux(i) )
    538546
    539547      !--Diagnostic tendencies
    … …  
    599607    qliq(i) = qliq(i) + dqlauto
    600608    qice(i) = qice(i) + dqiauto
    601     raincld(i) = raincld(i) - dqlauto * hum_to_flux(i)
    602     snowcld(i) = snowcld(i) - dqiauto * hum_to_flux(i)
     609    raincld(i) = MAX( 0. , raincld(i) - dqlauto * hum_to_flux(i) )
     610    snowcld(i) = MAX( 0. , snowcld(i) - dqiauto * hum_to_flux(i) )
    603611
    604612    !--Diagnostic tendencies
    … …  
    631639      !--Add tendencies
    632640      qliq(i) = qliq(i) + dqlrim
    633       snowcld(i) = snowcld(i) - dqlrim * hum_to_flux(i)
     641      snowcld(i) = MAX( 0. , snowcld(i) - dqlrim * hum_to_flux(i) )
    634642
    635643      !--Temperature adjustment with the release of latent
    … …  
    699707      !--Barrier to limit the melting flux to the clr snow flux in the mesh
    700708      dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / hum_to_flux(i))
    701    ENDIF
     709    ENDIF
    702710
    703711    !--In cloudy air
    … …  
    726734      dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt
    727735      dqstotmelt = dqsmelt_max
     736
    728737    ENDIF
    729738
    730739    !--Add tendencies
    731     rainclr(i) = rainclr(i) - dqsclrmelt * hum_to_flux(i)
    732     raincld(i) = raincld(i) - dqscldmelt * hum_to_flux(i)
    733     snowclr(i) = snowclr(i) + dqsclrmelt * hum_to_flux(i)
    734     snowcld(i) = snowcld(i) + dqscldmelt * hum_to_flux(i)
     740
     741    rainclr(i) = MAX( 0. , rainclr(i) - dqsclrmelt * hum_to_flux(i) )
     742    raincld(i) = MAX( 0. , raincld(i) - dqscldmelt * hum_to_flux(i) )
     743    snowclr(i) = MAX( 0. , snowclr(i) + dqsclrmelt * hum_to_flux(i) )
     744    snowcld(i) = MAX( 0. , snowcld(i) + dqscldmelt * hum_to_flux(i) )
     745
    735746
    736747    !--Temperature adjustment with the release of latent
    … …  
    749760  !--                  FREEZING
    750761  !---------------------------------------------------------
    751   !--Process through which rain freezes into snow. This is
    752   !--parameterized as an exponential decrease of the rain
    753   !--water content.
    754   !--The formula is homemade.
     762  !--Process through which rain freezes into snow.
     763  !-- We parameterize it as a 2 step process:
     764  !-- first: freezing following collision with ice crystals
     765  !-- second: immersion freezing following (inspired by Bigg 1953)
     766  !-- the latter is parameterized as an exponential decrease of the rain
     767  !-- water content with a homemade formulya
    755768  !--This is based on a caracteritic time of freezing, which
    756769  !--exponentially depends on temperature so that it is
    … …  
    758771  !--0 for RTT (=273.15 K).
    759772  !--NB.: this process needs a temperature adjustment
    760 
    761773  !--dqrfreez_max: maximum rain freezing so that temperature
    762774  !--              stays lower than 273 K [kg/kg]
    … …  
    768780  IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN
    769781
    770     !--Computed according to
    771     !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
    772     dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
     782     
     783    !--1st step: freezing following collision with ice crystals
     784    !--Sub-process of freezing which quantifies the collision between
     785    !--ice crystals in suspension and falling rain droplets.
     786    !--The rain droplets freeze, becoming graupel, and carrying
     787    !--the ice crystal (which acted as an ice nucleating particle).
     788    !--The formula is adapted from the riming formula.
     789    !-- it works only in the cloudy part
     790
     791    dqifreez = 0.
     792    dqrtotfreez_step1 = 0.
     793
     794    IF ((qice(i) .GT. 0.) .AND. (raincld(i) .GT. 0.)) THEN
     795      dqrclrfreez = 0.
     796      dqrcldfreez = 0.
     797
     798      !--Computed according to
     799      !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
     800      dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
    773801                         * ( 1. + RVTMP2 * qtot(i) ))
    774802 
    775     tau_freez = 1. / ( beta_freez &
    776               * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
    777 
    778     !--Initialisation
    779     dqrclrfreez = 0.
    780     dqrcldfreez = 0.
    781 
    782     !--In clear air
    783     IF ( rainclr(i) .GT. 0. ) THEN
    784       !--Exact solution of dqrain/dt = -qrain/tau_freez
    785       dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
    786     ENDIF
    787 
    788     !--In cloudy air
    789     IF ( raincld(i) .GT. 0. ) THEN
    790       !--Exact solution of dqrain/dt = -qrain/tau_freez
    791       dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
    792 
    793       !---------------------------------------------------------
    794       !--    FREEZING OF RAIN WITH CONTACT OF ICE CRYSTALS
    795       !---------------------------------------------------------
    796803      !--Sub-process of freezing which quantifies the collision between
    797804      !--ice crystals in suspension and falling rain droplets.
    … …  
    800807      !--The formula is adapted from the riming formula.
    801808
    802       !--The sticking efficacy is perfect.
     809      !--The collision efficiency is assumed unity
    803810      Eff_rain_ice = 1.
    804811      coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice
    805       !-- ATTENTION Double implicit version? + barriers if needed
    806812      !--Exact version, which does not need a barrier because of
    807       !--the exponential decrease
     813      !--the exponential decrease.
    808814      dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. )
    809815
    810       !--We add the two freezing processes in cloud
    811       dqrcldfreez = dqrcldfreez - dqifreez
    812     ENDIF
    813 
    814     !--Barriers
     816      !--We add the part of rain water that freezes, limited by a temperature barrier
     817      !-- this quantity is calculated assuming that the number of drop that freeze correspond to the number
     818      !-- of crystals collected (and assuming uniform distributions of ice crystals and rain drops)
     819      dqrcldfreez = MAX(dqifreez*rho_rain*r_rain/(rho_ice*r_ice),-raincld(i)/hum_to_flux(i))
     820      dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max)
     821      dqrtotfreez_step1 = dqrcldfreez
     822
     823      !--Add tendencies
     824      qice(i) = qice(i) + dqifreez
     825      raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) )
     826      snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) - dqifreez * hum_to_flux(i) )
     827      temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD &
     828                      / ( 1. + RVTMP2 * qtot(i) )
     829
     830    ENDIF
     831   
     832    !-- Second step immersion freezing of rain drops
     833    !-- with a homemade timeconstant depending on temperature
     834
     835    dqrclrfreez = 0.
     836    dqrcldfreez = 0.
     837    dqrtotfreez_step2 = 0.
     838    !--Computed according to
     839    !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q
     840
     841    dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD &
     842                         * ( 1. + RVTMP2 * qtot(i) ))
     843
     844   
     845    tau_freez = 1. / ( beta_freez &
     846              * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )
     847
     848
     849    !--In clear air
     850    IF ( rainclr(i) .GT. 0. ) THEN
     851      !--Exact solution of dqrain/dt = -qrain/tau_freez
     852      dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
     853    ENDIF
     854
     855    !--In cloudy air
     856    IF ( raincld(i) .GT. 0. ) THEN
     857      !--Exact solution of dqrain/dt = -qrain/tau_freez
     858      dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )
     859    ENDIF
     860
     861    !--temperature barrie step 2
    815862    !--It is activated if the total is LOWER than the max
    816863    !--because everything is negative
    817     dqrtotfreez = dqrclrfreez + dqrcldfreez
    818     IF ( dqrtotfreez .LT. dqrfreez_max ) THEN
     864    dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez
     865 
     866    IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN
    819867      !--We redistribute the max freezed rain keeping
    820868      !--the clear/cloud partition of the freezing rain
    821       dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez
    822       dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez
    823       dqrtotfreez = dqrfreez_max
     869      dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2
     870      dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2
     871      dqrtotfreez_step2 = dqrfreez_max
    824872    ENDIF
    825873
    826874
    827875    !--Add tendencies
    828     rainclr(i) = rainclr(i) + dqrclrfreez * hum_to_flux(i)
    829     raincld(i) = raincld(i) + dqrcldfreez * hum_to_flux(i)
    830     snowclr(i) = snowclr(i) - dqrclrfreez * hum_to_flux(i)
    831     snowcld(i) = snowcld(i) - dqrcldfreez * hum_to_flux(i)
     876    rainclr(i) = MAX( 0. , rainclr(i) + dqrclrfreez * hum_to_flux(i) )
     877    raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) )
     878    snowclr(i) = MAX( 0. , snowclr(i) - dqrclrfreez * hum_to_flux(i) )
     879    snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) )
     880
    832881
    833882    !--Temperature adjustment with the uptake of latent
    834883    !--heat because of freezing
    835     temp(i) = temp(i) - dqrtotfreez * RLMLT / RCPD &
     884    temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD &
    836885                      / ( 1. + RVTMP2 * qtot(i) )
    837886
     887    dqrtotfreez=dqrtotfreez_step1+dqrtotfreez_step2         
    838888    !--Diagnostic tendencies
    839889    dqrfreez(i) = dqrtotfreez / dtime
    840     dqsfreez(i) = - dqrtotfreez / dtime
     890    dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
    841891
    842892  ENDIF
Note: See TracChangeset for help on using the changeset viewer.