MODULE lmdz_lscp_ini

!----------------------------------------------------------
! This module initializes all the parameters needed for the
! lscp (large scale clouds and precipitation) routines.
!----------------------------------------------------------

   IMPLICIT NONE

   ! PARAMETERS for lscp_main:
   !-------------------------

   REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W
   !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RV, RG, RPI, EPS_W)

   REAL, SAVE, PROTECTED :: seuil_neb = 0.001       ! cloud fraction threshold: a cloud can precipitate when exceeded
   !$OMP THREADPRIVATE(seuil_neb)

   REAL, SAVE, PROTECTED :: min_neb_th = 1e-10      ! a cloud produced by bi-gaussian really exists when exceeded
   !$OMP THREADPRIVATE(min_neb_th)

   REAL, SAVE, PROTECTED :: min_frac_thermals = 1.e-10   ! minimum thermals fraction for use of bigaussian distribution
   !$OMP THREADPRIVATE(min_frac_thermals)

   REAL, SAVE, PROTECTED :: qlmax = 999.           ! maximum possible liquid water specific content in a mesh [kg/kg]
   !$OMP THREADPRIVATE(qlmax)

   REAL, SAVE, PROTECTED :: qimax = 999.           ! maximum possible ice water specific content in a mesh [kg/kg]
   !$OMP THREADPRIVATE(qimax)

   INTEGER, SAVE :: lunout, prt_level            ! Logical unit number and level for standard output
   !$OMP THREADPRIVATE(lunout,prt_level)

   INTEGER, SAVE, PROTECTED :: niter_lscp = 5      ! number of iterations to calculate autoconversion to precipitation
   !$OMP THREADPRIVATE(niter_lscp)

   INTEGER, SAVE, PROTECTED :: iflag_ice_thermo    ! activation of ice thermodynamics
   !$OMP THREADPRIVATE(iflag_ice_thermo)

   INTEGER, SAVE, PROTECTED :: iflag_cld_th        ! controls the properties of clouds in presence of thermals
   !$OMP THREADPRIVATE(iflag_cld_th)

   INTEGER, SAVE, PROTECTED :: iflag_evap_prec = 1 ! precipitation evaporation flag. 0: nothing, 1: "old way",
   ! 2: Max cloud fraction above to calculate the max of reevaporation
   ! >=4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only
   ! pay attention that iflag_evap_prec=4 may lead to unrealistic and overestimated
   ! evaporation. Use 5 instead
   !$OMP THREADPRIVATE(iflag_evap_prec)

   REAL, SAVE, PROTECTED :: t_coup = 234.0         ! temperature threshold which determines the phase
   ! for which the saturation vapor pressure is calculated
   !$OMP THREADPRIVATE(t_coup)
   REAL, SAVE, PROTECTED :: DDT0 = 0.01            ! iteration parameter
   !$OMP THREADPRIVATE(DDT0)

   REAL, SAVE, PROTECTED :: ztfondue = 278.15      ! parameter to calculate melting fraction of precipitation
   !$OMP THREADPRIVATE(ztfondue)

   REAL, SAVE, PROTECTED :: temp_nowater = 235.15  ! temperature below which liquid water no longer exists
   !$OMP THREADPRIVATE(temp_nowater)

   REAL, SAVE, PROTECTED :: a_tr_sca(4)          ! Variables for tracers temporary: alpha parameter for scavenging, 4 possible scavenging processes
   !$OMP THREADPRIVATE(a_tr_sca)

   REAL, SAVE, PROTECTED ::  min_frac_th_cld = 1.e-10   ! minimum thermal fraction to compute a thermal cloud fraction
   !$OMP THREADPRIVATE(min_frac_th_cld)

   LOGICAL, SAVE, PROTECTED :: ok_radocond_snow = .false. ! take into account the mass of ice precip in the cloud ice content seen by radiation
   !$OMP THREADPRIVATE(ok_radocond_snow)

   REAL, SAVE, PROTECTED :: t_glace_min = 258.0    ! lower-bound temperature parameter for cloud phase determination
   !$OMP THREADPRIVATE(t_glace_min)

   REAL, SAVE, PROTECTED :: t_glace_max = 273.15   ! upper-bound temperature parameter for cloud phase determination
   !$OMP THREADPRIVATE(t_glace_max)

   REAL, SAVE, PROTECTED :: exposant_glace = 1.0   ! parameter for cloud phase determination
   !$OMP THREADPRIVATE(exposant_glace)

   INTEGER, SAVE, PROTECTED :: iflag_vice = 0      ! which expression for ice crystall fall velocity
   !$OMP THREADPRIVATE(iflag_vice)

   INTEGER, SAVE, PROTECTED :: iflag_t_glace = 0   ! which expression for cloud phase partitioning
   !$OMP THREADPRIVATE(iflag_t_glace)

   INTEGER, SAVE, PROTECTED :: iflag_gammasat = 0              ! which threshold for homogeneous nucleation below -40oC
   !$OMP THREADPRIVATE(iflag_gammasat)

   INTEGER, SAVE, PROTECTED :: iflag_rain_incloud_vol = 0      ! use of volume cloud fraction for rain autoconversion
   !$OMP THREADPRIVATE(iflag_rain_incloud_vol)

   INTEGER, SAVE, PROTECTED :: iflag_bergeron = 0              ! bergeron effect for liquid precipitation treatment
   !$OMP THREADPRIVATE(iflag_bergeron)

   INTEGER, SAVE, PROTECTED :: iflag_fisrtilp_qsat = 0         ! qsat adjustment (iterative) during autoconversion
   !$OMP THREADPRIVATE(iflag_fisrtilp_qsat)

   INTEGER, SAVE, PROTECTED :: iflag_pdf = 0                   ! type of subgrid scale qtot pdf
   !$OMP THREADPRIVATE(iflag_pdf)

   INTEGER, SAVE, PROTECTED :: iflag_icefrac = 0               ! which phase partitioning function to use
   !$OMP THREADPRIVATE(iflag_icefrac)

   INTEGER, SAVE, PROTECTED :: iflag_autoconversion = 0        ! autoconversion option
   !$OMP THREADPRIVATE(iflag_autoconversion)

   LOGICAL, SAVE, PROTECTED :: ok_bug_phase_lscp = .true.      ! bug on phase partitioning after precipitation processes
   !$OMP THREADPRIVATE(ok_bug_phase_lscp)

   LOGICAL, SAVE, PROTECTED :: reevap_ice = .false.            ! no liquid precip for T< threshold
   !$OMP THREADPRIVATE(reevap_ice)

   REAL, SAVE, PROTECTED :: cld_lc_lsc = 2.6e-4                ! liquid autoconversion coefficient, stratiform rain
   !$OMP THREADPRIVATE(cld_lc_lsc)

   REAL, SAVE, PROTECTED :: cld_lc_con = 2.6e-4                ! liquid autoconversion coefficient, convective rain
   !$OMP THREADPRIVATE(cld_lc_con)

   REAL, SAVE, PROTECTED :: cld_tau_lsc = 3600.                ! liquid autoconversion timescale, stratiform rain
   !$OMP THREADPRIVATE(cld_tau_lsc)

   REAL, SAVE, PROTECTED :: cld_tau_con = 3600.                ! liquid autoconversion timescale, convective rain
   !$OMP THREADPRIVATE(cld_tau_con)

   REAL, SAVE, PROTECTED :: cld_expo_lsc = 2.                  ! liquid autoconversion threshold exponent, stratiform rain
   !$OMP THREADPRIVATE(cld_expo_lsc)

   REAL, SAVE, PROTECTED :: cld_expo_con = 2.                  ! liquid autoconversion threshold exponent, convective rain
   !$OMP THREADPRIVATE(cld_expo_con)

   REAL, SAVE, PROTECTED :: ffallv_lsc = 1.                    ! tuning coefficient crystal fall velocity, stratiform
   !$OMP THREADPRIVATE(ffallv_lsc)

   REAL, SAVE, PROTECTED :: ffallv_con = 1.                    ! tuning coefficient crystal fall velocity, convective
   !$OMP THREADPRIVATE(ffallv_con)

   REAL, SAVE, PROTECTED :: coef_eva = 2e-5                    ! tuning coefficient liquid precip evaporation
   !$OMP THREADPRIVATE(coef_eva)

   REAL, SAVE, PROTECTED :: coef_sub                         ! tuning coefficient ice precip sublimation
   !$OMP THREADPRIVATE(coef_sub)

   REAL, SAVE, PROTECTED :: expo_eva = 0.5                     ! tuning coefficient liquid precip evaporation
   !$OMP THREADPRIVATE(expo_eva)

   REAL, SAVE, PROTECTED :: expo_sub                         ! tuning coefficient ice precip sublimation
   !$OMP THREADPRIVATE(expo_sub)

   REAL, SAVE, PROTECTED :: cice_velo = 1.645                  ! factor in the ice fall velocity formulation. It is half the value of
   ! Heymsfield and Donner 1990 to concur with previous LMDZ versions
   !$OMP THREADPRIVATE(cice_velo)

   REAL, SAVE, PROTECTED ::  dice_velo = 0.16                  ! exponent in the ice fall velocity formulation
   !$OMP THREADPRIVATE(dice_velo)

   REAL, SAVE, PROTECTED :: dist_liq = 300.                    ! typical deph of cloud-top liquid layer in mpcs
   !$OMP THREADPRIVATE(dist_liq)

   REAL, SAVE, PROTECTED  :: tresh_cl = 0.0                  ! cloud fraction threshold for cloud top search
   !$OMP THREADPRIVATE(tresh_cl)

   !--Parameters for condensation and ice supersaturation

   LOGICAL, SAVE, PROTECTED :: ok_ice_supersat = .FALSE.        ! activates the condensation scheme that allows for ice supersaturation
   !$OMP THREADPRIVATE(ok_ice_supersat)

   LOGICAL, SAVE, PROTECTED :: ok_unadjusted_clouds = .FALSE.   ! if True, relax the saturation adjustment assumption for ice clouds
   !$OMP THREADPRIVATE(ok_unadjusted_clouds)

   LOGICAL, SAVE, PROTECTED :: ok_weibull_warm_clouds = .FALSE. ! if True, the weibull condensation scheme replaces the lognormal condensation scheme at positive temperatures
   !$OMP THREADPRIVATE(ok_weibull_warm_clouds)

   INTEGER, SAVE, PROTECTED :: iflag_cloud_sublim_pdf = 4       ! iflag for the distribution of water inside ice clouds
   !$OMP THREADPRIVATE(iflag_cloud_sublim_pdf)

   REAL, SAVE, PROTECTED :: depo_coef_cirrus = .7               ! [-] deposition coefficient for growth of ice crystals in cirrus clouds
   !$OMP THREADPRIVATE(depo_coef_cirrus)

   REAL, SAVE, PROTECTED :: capa_cond_cirrus = .5               ! [-] capacitance factor for growth/sublimation of ice crystals in cirrus clouds
   !$OMP THREADPRIVATE(capa_cond_cirrus)

   REAL, SAVE, PROTECTED :: std_subl_pdf_lscp = 2.              ! [%] standard deviation of the gaussian distribution of water inside ice clouds
   !$OMP THREADPRIVATE(std_subl_pdf_lscp)

   REAL, SAVE, PROTECTED :: mu_subl_pdf_lscp = 1./3.            ! [-] shape factor of the gamma distribution of water inside ice clouds
   !$OMP THREADPRIVATE(mu_subl_pdf_lscp)

   REAL, SAVE, PROTECTED :: beta_pdf_lscp = 1.E-3               ! [SI] tuning coefficient for the standard deviation of the PDF of water vapor in the clear sky region
   !$OMP THREADPRIVATE(beta_pdf_lscp)

   REAL, SAVE, PROTECTED :: temp_thresh_pdf_lscp = 189.         ! [K] factor for the PDF fit of water vapor in UTLS - below this temperature, water vapor is homogeneously distributed in the clear sky region
   !$OMP THREADPRIVATE(temp_thresh_pdf_lscp)

   REAL, SAVE, PROTECTED :: k0_pdf_lscp = 3.01                  ! [-] factor for the PDF fit of water vapor in UTLS
   !$OMP THREADPRIVATE(k0_pdf_lscp)

   REAL, SAVE, PROTECTED :: kappa_pdf_lscp = 0.0192             ! [] factor for the PDF fit of water vapor in UTLS
   !$OMP THREADPRIVATE(kappa_pdf_lscp)

   REAL, SAVE, PROTECTED :: std100_pdf_lscp = 4.08              ! [%] standard deviation at RHliq=100% of the PDF fit of water vapor in UTLS
   !$OMP THREADPRIVATE(std100_pdf_lscp)

   REAL, SAVE, PROTECTED :: a_homofreez = 2.349                 ! [-] factor for the Koop homogeneous freezing fit
   !$OMP THREADPRIVATE(a_homofreez)

   REAL, SAVE, PROTECTED :: b_homofreez = 259.                  ! [K] factor for the Koop homogeneous freezing fit
   !$OMP THREADPRIVATE(b_homofreez)

   REAL, SAVE, PROTECTED :: delta_hetfreez = 1.                 ! [-] value between 0 and 1 to simulate for heterogeneous freezing.
   !$OMP THREADPRIVATE(delta_hetfreez)

   REAL, SAVE, PROTECTED :: coef_mixing_lscp = 1.E-3            ! [-] tuning coefficient for the mixing process
   !$OMP THREADPRIVATE(coef_mixing_lscp)

   REAL, SAVE, PROTECTED :: coef_shear_lscp = 0.72              ! [-] additional coefficient for the shearing process (subprocess of the mixing process)
   !$OMP THREADPRIVATE(coef_shear_lscp)

   REAL, SAVE, PROTECTED :: chi_mixing_lscp = 1.                ! [-] factor for the macro distribution of ISSRs wrt clouds in a gridbox
   !$OMP THREADPRIVATE(chi_mixing_lscp)
   !--End of the parameters for condensation and ice supersaturation

   !--Parameters for poprecip and cloud phase
   LOGICAL, SAVE, PROTECTED :: ok_poprecip = .FALSE.           ! use the processes-oriented formulation of precipitations
   !$OMP THREADPRIVATE(ok_poprecip)

   LOGICAL, SAVE, PROTECTED :: ok_corr_vap_evasub = .FALSE.    ! use the corrected version of clear-sky water vapor for the evap / subl processes
   !$OMP THREADPRIVATE(ok_corr_vap_evasub)

   LOGICAL, SAVE, PROTECTED :: ok_growth_precip_deposition = .FALSE. ! allows growth of snowfall through vapor deposition in supersat. regions
   !$OMP THREADPRIVATE(ok_growth_precip_deposition)

   REAL, SAVE, PROTECTED :: cld_lc_lsc_snow = 2.e-5            ! snow autoconversion coefficient, stratiform. default from  Chaboureau and PInty 2006
   !$OMP THREADPRIVATE(cld_lc_lsc_snow)

   REAL, SAVE, PROTECTED :: cld_lc_con_snow = 2.e-5            ! snow autoconversion coefficient, convective
   !$OMP THREADPRIVATE(cld_lc_con_snow)

   REAL, SAVE, PROTECTED :: rain_int_min = 1.e-5               ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction
   !$OMP THREADPRIVATE(rain_int_min)

   REAL, SAVE, PROTECTED :: thresh_precip_frac = 1.E-6         ! precipitation fraction threshold [-]
   !$OMP THREADPRIVATE(thresh_precip_frac)

   REAL, SAVE, PROTECTED :: capa_crystal = 1.                  ! Crystal capacitance (shape factor) for lscp_icefrac_turb [-]
   !$OMP THREADPRIVATE(capa_crystal)

   LOGICAL, SAVE, PROTECTED :: fix_icnc = .FALSE.             ! fix ice crystal number concentration for cloud phase determination  
   !$OMP THREADPRIVATE(fix_icnc)

   REAL, SAVE, PROTECTED :: icnc0 = 10.                       ! prescribed ice crystal number concentration for cloud phase determination  [L-1]
   !$OMP THREADPRIVATE(icnc0)

   REAL, SAVE, PROTECTED :: naero5 = 0.5                       ! Number concentration of aerosol larger than 0.5 microns [scm-3]
   !$OMP THREADPRIVATE(naero5)

   REAL, SAVE, PROTECTED :: gamma_snwretro = 0.              ! Proportion of snow taken into account in ice retroaction in lscp_icefrac_turb [-]
   !$OMP THREADPRIVATE(gamma_snwretro)

   REAL, SAVE, PROTECTED :: gamma_mixth = 1.                 ! Tuning coeff for mixing with thermals/env in lscp_icefrac_turb [-]
   !$OMP THREADPRIVATE(gamma_mixth)

   REAL, SAVE, PROTECTED :: gamma_taud = 1.                  ! Tuning coeff for Lagrangian decorrelation timescale in lscp_icefrac_turb [-]
   !$OMP THREADPRIVATE(gamma_taud)

   REAL, SAVE, PROTECTED :: gamma_col = 1.                     ! Tuning coefficient for rain collection efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_col)

   REAL, SAVE, PROTECTED :: gamma_agg = 1.                     ! Tuning coefficient for snow aggregation efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_agg)

   REAL, SAVE, PROTECTED :: gamma_rim = 1.                     ! Tuning coefficient for riming efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_rim)

   REAL, SAVE, PROTECTED :: gamma_melt = 1.                    ! Tuning coefficient for snow melting efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_melt)

   REAL, SAVE, PROTECTED :: gamma_freez = 1.                   ! Tuning coefficient for rain collision freezing efficiency (poprecip) [-]
   !$OMP THREADPRIVATE(gamma_freez)

   REAL, SAVE, PROTECTED :: rho_rain = 1000.                   ! Rain density [kg/m3]
   !$OMP THREADPRIVATE(rho_rain)

   REAL, SAVE, PROTECTED :: rho_ice = 920.                     ! Ice crystal density (assuming spherical geometry) [kg/m3]
   !$OMP THREADPRIVATE(rho_ice)

   REAL, SAVE, PROTECTED :: r_rain = 500.E-6                   ! Rain droplets radius (poprecip) [m]
   !$OMP THREADPRIVATE(r_rain)

   REAL, SAVE, PROTECTED :: r_snow = 1.E-3                     ! Ice crystals radius (poprecip) [m]
   !$OMP THREADPRIVATE(r_snow)

   REAL, SAVE, PROTECTED :: tau_auto_snow_min = 1200.           ! Snow autoconversion minimal timescale (when liquid) [s]
   !$OMP THREADPRIVATE(tau_auto_snow_min)

   REAL, SAVE, PROTECTED :: tau_auto_snow_max = 3600.          ! Snow autoconversion minimal timescale (when only ice) [s]
   !$OMP THREADPRIVATE(tau_auto_snow_max)

   REAL, SAVE, PROTECTED :: expo_tau_auto_snow = 0.1          ! Snow autoconversion timescale exponent for icefrac dependency
   !$OMP THREADPRIVATE(expo_tau_auto_snow)

   REAL, SAVE, PROTECTED :: eps = 1.E-10                       ! Treshold 0 [-]
   !$OMP THREADPRIVATE(eps)

   REAL, SAVE, PROTECTED :: alpha_freez = 4.                   ! Slope of exponential for immersion freezing timescale [-]
   !$OMP THREADPRIVATE(alpha_freez)

   REAL, SAVE, PROTECTED :: beta_freez = 0.1                   ! Inv.time immersion freezing [s-1]
   !$OMP THREADPRIVATE(beta_freez)

   REAL, SAVE, PROTECTED :: rain_fallspeed = 4.                ! Rain fall velocity [m/s]
   !$OMP THREADPRIVATE(rain_fallspeed)

   REAL, SAVE, PROTECTED :: rain_fallspeed_clr               ! Rain fall velocity in clear sky [m/s]
   !$OMP THREADPRIVATE(rain_fallspeed_clr)

   REAL, SAVE, PROTECTED :: rain_fallspeed_cld               ! Rain fall velocity in cloudy sky [m/s]
   !$OMP THREADPRIVATE(rain_fallspeed_cld)

   REAL, SAVE, PROTECTED :: snow_fallspeed = 1.                ! Snow fall velocity [m/s]
   !$OMP THREADPRIVATE(snow_fallspeed)

   REAL, SAVE, PROTECTED :: snow_fallspeed_clr               ! Snow fall velocity in clear sky [m/s]
   !$OMP THREADPRIVATE(snow_fallspeed_clr)

   REAL, SAVE, PROTECTED :: snow_fallspeed_cld               ! Snow fall velocity in cloudy sky [m/s]
   !$OMP THREADPRIVATE(snow_fallspeed_cld)
   !--End of the parameters for poprecip

   ! Parameters for cloudth routines
   LOGICAL, SAVE, PROTECTED :: ok_lscp_mergecond = .false.     ! more consistent condensation stratiform and shallow convective clouds
   !$OMP THREADPRIVATE(ok_lscp_mergecond)

   LOGICAL, SAVE, PROTECTED :: ok_bug_deficit = .true.        ! bug in saturation deficit calculation (in old cloudth routines) found out by JYG in 02/2026 true by default.
   !$OMP THREADPRIVATE(ok_bug_deficit)

   INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert = 0          ! option for determining cloud fraction and content in convective boundary layers
   !$OMP THREADPRIVATE(iflag_cloudth_vert)

   INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert_noratqs = 0  ! option to control the width of gaussian distrib in a specific case
   !$OMP THREADPRIVATE(iflag_cloudth_vert_noratqs)

   REAL, SAVE, PROTECTED :: cloudth_ratqsmin = -1.             ! minimum ratqs in cloudth
   !$OMP THREADPRIVATE(cloudth_ratqsmin)

   REAL, SAVE, PROTECTED :: sigma1s_factor = 1.1               ! factor for standard deviation of gaussian distribution of environment
   !$OMP THREADPRIVATE(sigma1s_factor)

   REAL, SAVE, PROTECTED :: sigma2s_factor = 0.09              ! factor for standard deviation of gaussian distribution of thermals
   !$OMP THREADPRIVATE(sigma2s_factor)

   REAL, SAVE, PROTECTED :: sigma1s_power = 0.6                ! exponent for standard deviation of gaussian distribution of environment
   !$OMP THREADPRIVATE(sigma1s_power)

   REAL, SAVE, PROTECTED :: sigma2s_power = 0.5                ! exponent for standard deviation of gaussian distribution of thermals
   !$OMP THREADPRIVATE(sigma2s_power)

   REAL, SAVE, PROTECTED :: vert_alpha = 0.5                   ! tuning coefficient for standard deviation of gaussian distribution of thermals
   !$OMP THREADPRIVATE(vert_alpha)

   REAL, SAVE, PROTECTED :: vert_alpha_th = 0.5                ! tuning coefficient for standard deviation of gaussian distribution of thermals
   !$OMP THREADPRIVATE(vert_alpha_th)
   ! End of parameters for cloudth routines

   ! parameters for subgrid variability (ratqs) routines

   INTEGER, SAVE, PROTECTED :: iflag_ratqs = 1               ! control of ratqs option
   !$OMP THREADPRIVATE(iflag_ratqs)

   REAL, SAVE, PROTECTED :: ratqshaut = 0.3                  ! ratqs [-] value at the top of the atmosphere
   !$OMP THREADPRIVATE(ratqshaut)

   REAL, SAVE, PROTECTED :: ratqsbas = 0.01                  ! ratqs [-] value near the ground
   !$OMP THREADPRIVATE(ratqsbas)

   REAL, SAVE, PROTECTED :: ratqsdp = 20000.                 ! delta pressure [Pa] that controls the atan variation of ratqs with pressure
   !$OMP THREADPRIVATE(ratqsdp)

   REAL, SAVE, PROTECTED :: ratqsp0 = 50000.                 ! pressure [Pa] where the atan variation of ratqs with pressure occurs
   !$OMP THREADPRIVATE(ratqsp0)

   REAL, SAVE, PROTECTED :: tau_ratqs = 1800.                ! timescale [s] used when combining convective and stratiform ratqs
   !$OMP THREADPRIVATE(tau_ratqs)

   REAL, SAVE, PROTECTED :: a_ratqs_cv = 1.                  ! a_ratqs_cv and a_ratqs_wake control the contibution of cold pools to
   ! water subgrid variability
   !$OMP THREADPRIVATE(a_ratqs_cv)

   REAL, SAVE, PROTECTED :: a_ratqs_wake = 3.
   !$OMP THREADPRIVATE(a_ratqs_wake)

   REAL, SAVE, PROTECTED :: tau_var = 700.0                  ! relaxation time for variability in thermals [s]
   !$OMP THREADPRIVATE(tau_var)

   REAL, SAVE, PROTECTED :: fac_tau = 2.
   !$OMP THREADPRIVATE(fac_tau)

   REAL, SAVE, PROTECTED :: tau_cumul = 1000.
   !$OMP THREADPRIVATE(tau_cumul)

   INTEGER, SAVE, PROTECTED :: dqimpl = 1
   !$OMP THREADPRIVATE(dqimpl)

   REAL, ALLOCATABLE, SAVE :: povariance(:, :)
   !$OMP THREADPRIVATE(povariance)

   REAL, ALLOCATABLE, SAVE :: var_conv(:, :)
   !$OMP THREADPRIVATE(var_conv)

   ! Two parameters used for lmdz_lscp_old only
   INTEGER, SAVE, PROTECTED :: iflag_oldbug_fisrtilp = 0, fl_cor_ebil
   !$OMP THREADPRIVATE(iflag_oldbug_fisrtilp,fl_cor_ebil)

CONTAINS

      SUBROUTINE lscp_ini(dtime, klon, klev, iflag_thermals, iflag_cld_th_in,  &
                       iflag_ice_thermo_in, lunout_in, prt_level_in,           &
                       ok_ice_supersat_in, fl_cor_ebil_in, &
                       RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RVTMP2_in, &
                       RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in)

      USE ioipsl_getin_p_mod, ONLY: getin_p

      REAL, INTENT(IN)      :: dtime
      INTEGER, INTENT(IN)   :: klon, klev, iflag_thermals
      INTEGER, INTENT(IN)   :: iflag_cld_th_in, iflag_ice_thermo_in
      INTEGER, INTENT(IN)   :: lunout_in, prt_level_in, fl_cor_ebil_in
      LOGICAL, INTENT(IN)   :: ok_ice_supersat_in
      REAL, INTENT(IN)      :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in
      REAL, INTENT(IN)      :: RVTMP2_in, RTT_in, RD_in, RV_in, RG_in, RPI_in, EPS_W_in

      CHARACTER(len=20) :: modname = 'lscp_ini_mod'
      CHARACTER(len=80) :: abort_message

      INTEGER               :: i, k

      lunout = lunout_in
      prt_level = prt_level_in
      fl_cor_ebil = fl_cor_ebil_in
      ok_ice_supersat = ok_ice_supersat_in
      iflag_cld_th=iflag_cld_th_in
      iflag_ice_thermo=iflag_ice_thermo_in

      RG = RG_in
      RD = RD_in
      RV = RV_in
      RCPD = RCPD_in
      RLVTT = RLVTT_in
      RLSTT = RLSTT_in
      RLMLT = RLMLT_in
      RTT = RTT_in
      RV = RV_in
      RVTMP2 = RVTMP2_in
      RPI = RPI_in
      EPS_W = EPS_W_in

      IF (iflag_thermals == 15 .OR. iflag_thermals == 16) THEN
         dqimpl = -1
      END IF

      CALL getin_p('niter_lscp', niter_lscp)
      CALL getin_p('iflag_evap_prec', iflag_evap_prec)
      CALL getin_p('seuil_neb', seuil_neb)

      CALL getin_p('oliqmax', qlmax)
      CALL getin_p('qlmax', qlmax)

      CALL getin_p('oicemax', qimax)
      CALL getin_p('qimax', qimax)

      CALL getin_p('ok_radocond_snow', ok_radocond_snow)
      CALL getin_p('t_glace_max', t_glace_max)
      CALL getin_p('t_glace_min', t_glace_min)
      CALL getin_p('exposant_glace', exposant_glace)
      CALL getin_p('iflag_vice', iflag_vice)
      CALL getin_p('iflag_t_glace', iflag_t_glace)
      CALL getin_p('iflag_gammasat', iflag_gammasat)
      CALL getin_p('iflag_rain_incloud_vol', iflag_rain_incloud_vol)
      CALL getin_p('iflag_bergeron', iflag_bergeron)
      CALL getin_p('iflag_fisrtilp_qsat', iflag_fisrtilp_qsat)
      CALL getin_p('iflag_pdf', iflag_pdf)
      CALL getin_p('iflag_icefrac', iflag_icefrac)
      CALL getin_p('reevap_ice', reevap_ice)
      CALL getin_p('cld_lc_lsc', cld_lc_lsc)
      CALL getin_p('cld_lc_con', cld_lc_con)
      CALL getin_p('cld_lc_lsc_snow', cld_lc_lsc_snow)
      CALL getin_p('cld_lc_con_snow', cld_lc_con_snow)
      CALL getin_p('cld_tau_lsc', cld_tau_lsc)
      CALL getin_p('cld_tau_con', cld_tau_con)
      CALL getin_p('cld_expo_lsc', cld_expo_lsc)
      CALL getin_p('cld_expo_con', cld_expo_con)
      CALL getin_p('ffallv_lsc', ffallv_lsc)
      CALL getin_p('ffallv_lsc', ffallv_con)

      ! for poprecip and cloud phase
      CALL getin_p('coef_eva', coef_eva)
      coef_sub = coef_eva
      CALL getin_p('coef_eva_i', coef_sub)
      CALL getin_p('coef_sub', coef_sub)
      CALL getin_p('expo_eva', expo_eva)
      expo_sub = expo_eva
      CALL getin_p('expo_sub', expo_sub)
      CALL getin_p('iflag_autoconversion', iflag_autoconversion)
      CALL getin_p('dist_liq', dist_liq)
      CALL getin_p('tresh_cl', tresh_cl)
      CALL getin_p('fix_icnc', fix_icnc)
      CALL getin_p('icnc0', icnc0)
      CALL getin_p('capa_crystal', capa_crystal)
      CALL getin_p('naero5', naero5)
      CALL getin_p('gamma_snwretro', gamma_snwretro)
      CALL getin_p('gamma_taud', gamma_taud)
      CALL getin_p('gamma_mixth', gamma_mixth)
      CALL getin_p('iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp)
      CALL getin_p('temp_nowater', temp_nowater)
      CALL getin_p('ok_bug_phase_lscp', ok_bug_phase_lscp)
      CALL getin_p('ok_poprecip', ok_poprecip)
      CALL getin_p('ok_corr_vap_evasub', ok_corr_vap_evasub)
      CALL getin_p('ok_growth_precip_deposition', ok_growth_precip_deposition)
      CALL getin_p('rain_int_min', rain_int_min)
      CALL getin_p('gamma_agg', gamma_agg)
      CALL getin_p('gamma_col', gamma_col)
      CALL getin_p('gamma_rim', gamma_rim)
      CALL getin_p('gamma_freez', gamma_freez)
      CALL getin_p('gamma_melt', gamma_melt)
      CALL getin_p('tau_auto_snow_max', tau_auto_snow_max)
      CALL getin_p('tau_auto_snow_min', tau_auto_snow_min)
      CALL getin_p('expo_tau_auto_snow', expo_tau_auto_snow)
      CALL getin_p('alpha_freez', alpha_freez)
      CALL getin_p('beta_freez', beta_freez)
      CALL getin_p('r_snow', r_snow)
      CALL getin_p('rain_fallspeed', rain_fallspeed)
      rain_fallspeed_clr = rain_fallspeed
      rain_fallspeed_cld = rain_fallspeed
      CALL getin_p('rain_fallspeed_clr', rain_fallspeed_clr)
      CALL getin_p('rain_fallspeed_cld', rain_fallspeed_cld)
      CALL getin_p('snow_fallspeed', snow_fallspeed)
      snow_fallspeed_clr = snow_fallspeed
      snow_fallspeed_cld = snow_fallspeed
      CALL getin_p('snow_fallspeed_clr', snow_fallspeed_clr)
      CALL getin_p('snow_fallspeed_cld', snow_fallspeed_cld)

      ! for condensation and ice supersaturation
      CALL getin_p('ok_unadjusted_clouds', ok_unadjusted_clouds)
      CALL getin_p('ok_weibull_warm_clouds', ok_weibull_warm_clouds)
      CALL getin_p('iflag_cloud_sublim_pdf', iflag_cloud_sublim_pdf)
      CALL getin_p('depo_coef_cirrus', depo_coef_cirrus)
      CALL getin_p('capa_cond_cirrus', capa_cond_cirrus)
      CALL getin_p('std_subl_pdf_lscp', std_subl_pdf_lscp)
      CALL getin_p('mu_subl_pdf_lscp', mu_subl_pdf_lscp)
      CALL getin_p('beta_pdf_lscp', beta_pdf_lscp)
      CALL getin_p('temp_thresh_pdf_lscp', temp_thresh_pdf_lscp)
      CALL getin_p('k0_pdf_lscp', k0_pdf_lscp)
      CALL getin_p('kappa_pdf_lscp', kappa_pdf_lscp)
      CALL getin_p('std100_pdf_lscp', std100_pdf_lscp)
      CALL getin_p('a_homofreez', a_homofreez)
      CALL getin_p('b_homofreez', b_homofreez)
      CALL getin_p('delta_hetfreez', delta_hetfreez)
      CALL getin_p('coef_mixing_lscp', coef_mixing_lscp)
      CALL getin_p('coef_shear_lscp', coef_shear_lscp)
      CALL getin_p('chi_mixing_lscp', chi_mixing_lscp)

      ! for cloudth routines
      CALL getin_p('ok_lscp_mergecond', ok_lscp_mergecond)
      CALL getin_p('ok_bug_deficit', ok_bug_deficit)
      CALL getin_p('iflag_cloudth_vert', iflag_cloudth_vert)
      CALL getin_p('cloudth_ratqsmin', cloudth_ratqsmin)
      CALL getin_p('cloudth_sigma1s_factor', sigma1s_factor)
      CALL getin_p('cloudth_sigma1s_power', sigma1s_power)
      CALL getin_p('cloudth_sigma2s_factor', sigma2s_factor)
      CALL getin_p('cloudth_sigma2s_power', sigma2s_power)
      CALL getin_p('cloudth_vert_alpha', vert_alpha)
      vert_alpha_th = vert_alpha
      CALL getin_p('cloudth_vert_alpha_th', vert_alpha_th)
      CALL getin_p('iflag_cloudth_vert_noratqs', iflag_cloudth_vert_noratqs)

      ! for subgrid water variability (ratqs) routines
      CALL getin_p('iflag_ratqs', iflag_ratqs)
      CALL getin_p('ratqshaut', ratqshaut)
      CALL getin_p('ratqsbas', ratqsbas)
      CALL getin_p('ratqsdp', ratqsdp)
      CALL getin_p('ratqsp0', ratqsp0)
      CALL getin_p('tau_ratqs', tau_ratqs)
      CALL getin_p('tau_var', tau_var)
      CALL getin_p('fac_tau', fac_tau)
      CALL getin_p('tau_cumul', tau_cumul)
      CALL getin_p('a_ratqs_wake', a_ratqs_wake)
      CALL getin_p('a_ratqs_cv', a_ratqs_cv)

      WRITE (lunout, *) 'lscp_ini, niter_lscp:', niter_lscp
      WRITE (lunout, *) 'lscp_ini, iflag_evap_prec:', iflag_evap_prec
      WRITE (lunout, *) 'lscp_ini, seuil_neb:', seuil_neb
      WRITE (lunout, *) 'lscp_ini, qlmax:', qlmax
      WRITE (lunout, *) 'lscp_ini, qimax:', qimax
      WRITE (lunout, *) 'lscp_ini, ok_radocond_snow:', ok_radocond_snow
      WRITE (lunout, *) 'lscp_ini, t_glace_max:', t_glace_max
      WRITE (lunout, *) 'lscp_ini, t_glace_min:', t_glace_min
      WRITE (lunout, *) 'lscp_ini, exposant_glace:', exposant_glace
      WRITE (lunout, *) 'lscp_ini, iflag_vice:', iflag_vice
      WRITE (lunout, *) 'lscp_ini, iflag_t_glace:', iflag_t_glace
      WRITE (lunout, *) 'lscp_ini, iflag_gammasat:', iflag_gammasat
      WRITE (lunout, *) 'lscp_ini, iflag_rain_incloud_vol:', iflag_rain_incloud_vol
      WRITE (lunout, *) 'lscp_ini, iflag_bergeron:', iflag_bergeron
      WRITE (lunout, *) 'lscp_ini, iflag_fisrtilp_qsat:', iflag_fisrtilp_qsat
      WRITE (lunout, *) 'lscp_ini, iflag_pdf', iflag_pdf
      WRITE (lunout, *) 'lscp_ini, iflag_icefrac', iflag_icefrac
      WRITE (lunout, *) 'lscp_ini, reevap_ice', reevap_ice
      WRITE (lunout, *) 'lscp_ini, cld_lc_lsc', cld_lc_lsc
      WRITE (lunout, *) 'lscp_ini, cld_lc_con', cld_lc_con
      WRITE (lunout, *) 'lscp_ini, cld_lc_lsc_snow', cld_lc_lsc_snow
      WRITE (lunout, *) 'lscp_ini, cld_lc_con_snow', cld_lc_con_snow
      WRITE (lunout, *) 'lscp_ini, cld_tau_lsc', cld_tau_lsc
      WRITE (lunout, *) 'lscp_ini, cld_tau_con', cld_tau_con
      WRITE (lunout, *) 'lscp_ini, cld_expo_lsc', cld_expo_lsc
      WRITE (lunout, *) 'lscp_ini, cld_expo_con', cld_expo_con
      WRITE (lunout, *) 'lscp_ini, ffallv_lsc', ffallv_lsc
      WRITE (lunout, *) 'lscp_ini, ffallv_con', ffallv_con
      WRITE (lunout, *) 'lscp_ini, coef_eva', coef_eva
      WRITE (lunout, *) 'lscp_ini, coef_sub', coef_sub
      WRITE (lunout, *) 'lscp_ini, expo_eva', expo_eva
      WRITE (lunout, *) 'lscp_ini, expo_sub', expo_sub
      WRITE (lunout, *) 'lscp_ini, iflag_autoconversion', iflag_autoconversion
      WRITE (lunout, *) 'lscp_ini, dist_liq', dist_liq
      WRITE (lunout, *) 'lscp_ini, tresh_cl', tresh_cl
      WRITE (lunout, *) 'lscp_ini, capa_crystal', capa_crystal
      WRITE (lunout, *) 'lscp_ini, naero5', naero5
      WRITE (lunout, *) 'lscp_ini, fix_icnc', fix_icnc
      WRITE (lunout, *) 'lscp_ini, icnc0', icnc0
      WRITE (lunout, *) 'lscp_ini, gamma_snwretro', gamma_snwretro
      WRITE (lunout, *) 'lscp_ini, gamma_mixth', gamma_mixth
      WRITE (lunout, *) 'lscp_ini, gamma_taud', gamma_taud
      WRITE (lunout, *) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp
      WRITE (lunout, *) 'lscp_ini, fl_cor_ebil', fl_cor_ebil
      WRITE (lunout, *) 'lscp_ini, temp_nowater', temp_nowater
      WRITE (lunout, *) 'lscp_ini, ok_bug_phase_lscp', ok_bug_phase_lscp
      WRITE (lunout, *) 'lscp_ini, iflag_cld_th', iflag_cld_th
      WRITE (lunout, *) 'lscp_ini, iflag_ice_thermo', iflag_ice_thermo

      ! for poprecip
      WRITE (lunout, *) 'lscp_ini, ok_poprecip', ok_poprecip
      WRITE (lunout, *) 'lscp_ini, ok_corr_vap_evasub', ok_corr_vap_evasub
      WRITE (lunout, *) 'lscp_ini, ok_growth_precip_deposition', ok_growth_precip_deposition
      WRITE (lunout, *) 'lscp_ini, rain_int_min:', rain_int_min
      WRITE (lunout, *) 'lscp_ini, gamma_agg:', gamma_agg
      WRITE (lunout, *) 'lscp_ini, gamma_col:', gamma_col
      WRITE (lunout, *) 'lscp_ini, gamma_rim:', gamma_rim
      WRITE (lunout, *) 'lscp_ini, gamma_freez:', gamma_freez
      WRITE (lunout, *) 'lscp_ini, gamma_melt:', gamma_melt
      WRITE (lunout, *) 'lscp_ini, tau_auto_snow_max:', tau_auto_snow_max
      WRITE (lunout, *) 'lscp_ini, tau_auto_snow_min:', tau_auto_snow_min
      WRITE (lunout, *) 'lscp_ini, expo_tau_auto_snow:', expo_tau_auto_snow
      WRITE (lunout, *) 'lscp_ini, r_snow:', r_snow
      WRITE (lunout, *) 'lscp_ini, alpha_freez:', alpha_freez
      WRITE (lunout, *) 'lscp_ini, beta_freez:', beta_freez
      WRITE (lunout, *) 'lscp_ini, rain_fallspeed_clr:', rain_fallspeed_clr
      WRITE (lunout, *) 'lscp_ini, rain_fallspeed_cld:', rain_fallspeed_cld
      WRITE (lunout, *) 'lscp_ini, snow_fallspeed_clr:', snow_fallspeed_clr
      WRITE (lunout, *) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld
      ! for condensation and ice supersaturation
      WRITE (lunout, *) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat
      WRITE (lunout, *) 'lscp_ini, ok_unadjusted_clouds:', ok_unadjusted_clouds
      WRITE (lunout, *) 'lscp_ini, ok_weibull_warm_clouds:', ok_weibull_warm_clouds
      WRITE (lunout, *) 'lscp_ini, iflag_cloud_sublim_pdf:', iflag_cloud_sublim_pdf
      WRITE (lunout, *) 'lscp_ini, depo_coef_cirrus:', depo_coef_cirrus
      WRITE (lunout, *) 'lscp_ini, capa_cond_cirrus:', capa_cond_cirrus
      WRITE (lunout, *) 'lscp_ini, std_subl_pdf_lscp:', std_subl_pdf_lscp
      WRITE (lunout, *) 'lscp_ini, mu_subl_pdf_lscp:', mu_subl_pdf_lscp
      WRITE (lunout, *) 'lscp_ini, beta_pdf_lscp:', beta_pdf_lscp
      WRITE (lunout, *) 'lscp_ini, temp_thresh_pdf_lscp:', temp_thresh_pdf_lscp
      WRITE (lunout, *) 'lscp_ini, k0_pdf_lscp:', k0_pdf_lscp
      WRITE (lunout, *) 'lscp_ini, kappa_pdf_lscp:', kappa_pdf_lscp
      WRITE (lunout, *) 'lscp_ini, std100_pdf_lscp:', std100_pdf_lscp
      WRITE (lunout, *) 'lscp_ini, a_homofreez:', a_homofreez
      WRITE (lunout, *) 'lscp_ini, b_homofreez:', b_homofreez
      WRITE (lunout, *) 'lscp_ini, delta_hetfreez', delta_hetfreez
      WRITE (lunout, *) 'lscp_ini, coef_mixing_lscp:', coef_mixing_lscp
      WRITE (lunout, *) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp
      WRITE (lunout, *) 'lscp_ini, chi_mixing_lscp:', chi_mixing_lscp
      ! for cloudth routines
      WRITE (lunout, *) 'lscp_ini, ok_lscp_mergecond:', ok_lscp_mergecond
      WRITE (lunout, *) 'lscp_ini, ok_bug_deficit:', ok_bug_deficit
      WRITE (lunout, *) 'lscp_ini, iflag_cloudth_vert:', iflag_cloudth_vert
      WRITE (lunout, *) 'lscp_ini, cloudth_ratqsmin:', cloudth_ratqsmin
      WRITE (lunout, *) 'lscp_ini, cloudth_sigma1s_factor:', sigma1s_factor
      WRITE (lunout, *) 'lscp_ini, cloudth_sigma1s_power:', sigma1s_power
      WRITE (lunout, *) 'lscp_ini, cloudth_sigma2s_factor:', sigma2s_factor
      WRITE (lunout, *) 'lscp_ini, cloudth_sigma2s_power:', sigma2s_power
      WRITE (lunout, *) 'lscp_ini, cloudth_vert_alpha:', vert_alpha
      WRITE (lunout, *) 'lscp_ini, cloudth_vert_alpha_th:', vert_alpha_th
      WRITE (lunout, *) 'lscp_ini, iflag_cloudth_vert_noratqs:', iflag_cloudth_vert_noratqs
      ! for subgrid water variability (ratqs) routines
      WRITE (lunout, *) 'lscp_ini, iflag_ratqs:', iflag_ratqs
      WRITE (lunout, *) 'lscp_ini, ratqshaut:', ratqshaut
      WRITE (lunout, *) 'lscp_ini, ratqsbas:', ratqsbas
      WRITE (lunout, *) 'lscp_ini, ratqsdp:', ratqsdp
      WRITE (lunout, *) 'lscp_ini, ratqsp0', ratqsp0
      WRITE (lunout, *) 'lscp_ini, tau_ratqs:', tau_ratqs
      WRITE (lunout, *) 'lscp_ini, tau_var:', tau_var
      WRITE (lunout, *) 'lscp_ini, fac_tau:', fac_tau
      WRITE (lunout, *) 'lscp_ini, tau_cumulu:', tau_cumul
      WRITE (lunout, *) 'lscp_ini, a_ratqs_wake:', a_ratqs_wake
      WRITE (lunout, *) 'lscp_ini, a_ratqs_cv:', a_ratqs_cv

      ! check consistency for cloud phase partitioning options

      IF ((iflag_icefrac .GE. 2) .AND. (.NOT. ok_lscp_mergecond)) THEN
         abort_message = 'in lscp, iflag_icefrac .GE. 2 works only if ok_lscp_mergecond=.TRUE.'
         CALL abort_physic(modname, abort_message, 1)
      END IF

      ! check for precipitation sub-time steps
      IF (ABS(dtime/REAL(niter_lscp) - 360.0) .GT. 0.001) THEN
         WRITE (lunout, *) 'lscp: it is not expected, see Z.X.Li', dtime
         WRITE (lunout, *) 'I would prefer a 6 min sub-timestep'
      END IF

      ! check consistency between numerical resolution of autoconversion
      ! and other options

      IF ((iflag_autoconversion .EQ. 2) .AND. .NOT. ok_poprecip) THEN
         IF ((iflag_vice .NE. 0) .OR. (niter_lscp .GT. 1)) THEN
            abort_message = 'in lscp, iflag_autoconversion=2 requires iflag_vice=0 and niter_lscp=1'
            CALL abort_physic(modname, abort_message, 1)
         END IF
      END IF

      IF (ok_weibull_warm_clouds .AND. .NOT. ok_ice_supersat) THEN
         abort_message = 'in lscp, ok_weibull_warm_clouds=y needs ok_ice_supersat=y'
         CALL abort_physic(modname, abort_message, 1)
      END IF

      IF (ok_unadjusted_clouds .AND. .NOT. ok_ice_supersat) THEN
         abort_message = 'in lscp, ok_unadjusted_clouds=y needs ok_ice_supersat=y'
         CALL abort_physic(modname, abort_message, 1)
      END IF

      ! Temporary initialisation
      a_tr_sca(1) = -0.5
      a_tr_sca(2) = -0.5
      a_tr_sca(3) = -0.5
      a_tr_sca(4) = -0.5

      ! Initialisation of variances of subgrid variability of water
      ! for the prognostic parameterisation
      ! Works only in 1D so far, but should become a state variable
      ! of the model and an argument of the routines

      allocate (povariance(klon, klev))
      allocate (var_conv(klon, klev))

      IF (klon .EQ. 1) THEN
         DO k = 1, klev
            DO i = 1, klon
               povariance(i, k) = 0.00000000001
               var_conv(i, k) = 0.00000000001
            END DO
         END DO
      ELSE
         PRINT *, 'variance of subgrid variability of water should be conserved in restart files'
      END IF

      RETURN

   END SUBROUTINE lscp_ini

END MODULE lmdz_lscp_ini
