!$gpum horizontal klon
MODULE lmdz_call_lscp

   IMPLICIT NONE

CONTAINS

   SUBROUTINE call_lscp(klon, klev, nbsrf, is_ter, is_lic, &
                        abortphy, flag_inhib_tend, itap, &
                        nqo, dtime, missing_val, ok_new_lscp, &
                        paprs, pplay, omega, temp, qt, &
                        zmasse, ptconv, ratqsc, ratqs, ratqs_inter_, sigma_qtherm, &
                        qtc_cv, sigt_cv, detrain_cv, fm_cv, fqd, fqcomp, &
                        sigd, qsat, ratio_ql_qtot, ratio_qi_qtot, &
                        pctsrf, s_pblh, zstd, wake_s, wake_deltaq, &
                        d_t, d_q, d_ql, d_qi, rneb, rneblsvol, &
                        pfraclr, pfracld, &
                        cldfraliq, cldfraliqth, &
                        sigma2_icefracturb, sigma2_icefracturbth, &
                        mean_icefracturb, mean_icefracturbth, &
                        radocond, radicefrac, rain_lsc, snow_lsc, &
                        rain_fall, snow_fall, &
                        frac_impa, frac_nucl, beta, &
                        pfrac_impa, pfrac_nucl, pfrac_1nucl, &
                        prfl, psfl, rhcl, qta, fraca, &
                        tv, pspsk, tla, thl, wth, &
                        distcltop, temp_cltop, &
                        tke, tke_dissip, &
                        fm_therm, entr_therm, detr_therm, &
                        cell_area, &
                        cf, rvc, u, v, &
                        qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
                        dcf_sub, dcf_con, dcf_mix, &
                        dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, &
                        dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
                        Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi, &
                        dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
                        cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv, &
                        qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, &
                        dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim, &
                        dqsmelt, dqsfreez, rh, rhl, rhi)

      USE lmdz_lscp_main, ONLY: lscp
      USE lmdz_lscp_old, ONLY: fisrtilp, fisrtilp_first
      USE lmdz_lscp_subgridvarq, ONLY: ratqs_main, ratqs_main_first
      USE lmdz_lscp_ini, ONLY: qlmax, qimax, RTT
      USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf
      USE add_phys_tend_mod, ONLY: add_phys_tend, prt_enerbil
      USE phys_local_var_mod, ONLY: t_seri, q_seri, ql_seri, qs_seri

      IMPLICIT NONE

!============================================================================
! call_lscp in the main interface between the LMDZ physics monitor physiq_mod
! and the large scale clouds and precipitation (lscp) parameterizations
!
! The aim of lscp routines is to compute the cloud fraction, cloud water 
! contents, precipitation and temperature and water tendencies from the
! 'large scale' (i.e. not from deep convection)
!
! contact: Etienne Vignon etienne.vignon@lmd.ipsl.fr
!============================================================================

!===============================================================================
! VARIABLES DECLARATION
!===============================================================================

! INPUT VARIABLES:
!-----------------

      INTEGER, INTENT(IN)   :: klon, klev       ! number of horizontal grid points and vertical levels
      INTEGER, INTENT(IN)   :: nbsrf, is_ter, is_lic ! number of subgrid tiles and indices for land and landice
      INTEGER, INTENT(IN)   :: abortphy, flag_inhib_tend, itap ! flag and physics time counter for add_phys_tend
      INTEGER, INTENT(IN)   :: nqo             ! number of water species
      LOGICAL, INTENT(IN)   :: ok_new_lscp     ! flag that controls the version of lscp code used

      REAL, INTENT(IN)   :: dtime           ! time step [s]
      REAL, INTENT(IN)   :: missing_val     ! missing value for output

      REAL, DIMENSION(klon, klev + 1), INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: zmasse          ! masse of layers [kg/m2]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: temp            ! temperature (K)
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: ratio_ql_qtot   ! ratio ql/qt at the beginning of physics time step [-]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: ratio_qi_qtot   ! ratio qi/qt at the beginning of physics time step [-]
      REAL, DIMENSION(klon, klev + 1), INTENT(IN)   :: tke             ! turbulent kinetic energy [m2/s2]
      REAL, DIMENSION(klon, klev + 1), INTENT(IN)   :: tke_dissip      ! TKE dissipation [m2/s3]

      LOGICAL, DIMENSION(klon, klev), INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active

      !Inputs associated with thermal plumes

      REAL, DIMENSION(klon, klev), INTENT(IN)   :: tv              ! virtual potential temperature [K]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: qta             ! specific humidity within thermals [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: fraca           ! fraction of thermals within the mesh [-]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: pspsk           ! exner potential (p/100000)**(R/cp)
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: tla             ! liquid potential temperature within thermals [K]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: wth             ! vertical velocity within thermals [m/s]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: entr_therm      ! thermal plume entrainment rate * dz [kg/s/m2]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: detr_therm      ! thermal plume detrainment rate * dz [kg/s/m2]
      REAL, DIMENSION(klon, klev + 1), INTENT(IN)   :: fm_therm        ! convective mass flux of thermals [kg/s/m2]

      ! Inputs associated with subgrid water variability

      REAL, DIMENSION(klon, klev), INTENT(IN)   :: wake_deltaq     ! difference in humidity between wakes and environment [kg/kg]
      REAL, DIMENSION(klon), INTENT(IN)   :: wake_s          ! wake fraction area [-]
      REAL, DIMENSION(klon, nbsrf), INTENT(IN)   :: pctsrf          ! fraction of each subgrid tiles [0-1]
      REAL, DIMENSION(klon), INTENT(IN)   :: s_pblh          ! boundary layer height [m]
      REAL, DIMENSION(klon), INTENT(IN)   :: zstd            ! sub grid orography standard deviation [m]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: ratqsc          ! normalized subgrid distribution of water due to deep convection
      real, dimension(klon, klev), intent(in)   :: qtc_cv          ! total specific moisture in convective saturated draughts  [kg/kg]
      real, dimension(klon, klev), intent(in)   :: sigt_cv         ! surface fraction of convective saturated draughts [-]
      real, dimension(klon, klev), intent(in)   :: detrain_cv      ! deep convection detrainment of specific humidity variance [kg/s/m2*(kg/kg)^2]
      real, dimension(klon, klev), intent(in)   :: fm_cv           ! deep convective mass flux [kg/s/m2]
      real, dimension(klon, klev), intent(in)   :: fqd             ! specific humidity tendency due to convective precip [kg/kg/s]
      real, dimension(klon, klev), intent(in)   :: fqcomp          ! specific humidity tendency due to convective mixed draughts [kg/ks/s]
      real, dimension(klon), intent(in)   :: sigd            ! fractional area covered by unsaturated convective downdrafts [-]

      ! INPUT/OUTPUT variables
      !------------------------
      REAL, DIMENSION(klon, klev), INTENT(INOUT)   :: qsat            ! saturation specific humidity [kg/kg] from/to the physics
      REAL, DIMENSION(klon, klev), INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
      REAL, DIMENSION(klon, klev), INTENT(INOUT)   :: pfrac_impa   ! product of impaction scavenging coeff.
      REAL, DIMENSION(klon, klev), INTENT(INOUT)   :: pfrac_nucl   ! product of nucleation scavenging coeff.
      REAL, DIMENSION(klon, klev), INTENT(INOUT)   :: pfrac_1nucl  ! product of nucleation scavenging coeff. (alpha=1)
      REAL, DIMENSION(klon), INTENT(INOUT)   :: rain_fall    ! total surface rainfall [kg/m2/s]
      REAL, DIMENSION(klon), INTENT(INOUT)   :: snow_fall    ! total surface snowfall [kg/m2/s]

      ! INPUT/OUTPUT condensation and ice supersaturation
      !--------------------------------------------------
      REAL, DIMENSION(klon, klev), INTENT(INOUT):: cf          ! cloud fraction [-]
      REAL, DIMENSION(klon, klev), INTENT(INOUT):: rvc         ! cloudy water vapor to total water vapor ratio [-]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: u           ! eastward wind [m/s]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: v           ! northward wind [m/s]
      REAL, DIMENSION(klon), INTENT(IN)   :: cell_area        ! area of each cell [m2]

      ! INPUT/OUTPUT aviation
      !--------------------------------------------------
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
      REAL, DIMENSION(klon, klev), INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]

      ! OUTPUT variables
      !-------------------

      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: d_t              ! temperature increment [K]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: rneb             ! cloud fraction [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: cldfraliq        ! liquid fraction of cloud fraction [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: cldfraliqth      ! liquid fraction of cloud fraction in thermals [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: sigma2_icefracturbth ! Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: mean_icefracturbth   ! Mean of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
      REAL, DIMENSION(klon), INTENT(OUT)  :: rain_lsc             ! surface large-scale rainfall [kg/s/m2]
      REAL, DIMENSION(klon), INTENT(OUT)  :: snow_lsc             ! surface large-scale snowfall [kg/s/m2]
      REAL, DIMENSION(klon, klev + 1), INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
      REAL, DIMENSION(klon, klev + 1), INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: beta             ! conversion rate of condensed water
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: rh              ! relative humidity [0-1], wrt liq if T>0C, wrt ice if T<0C
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: rhl             ! relative humidity [0-1] wrt liq
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: rhi             ! relative humidity [0-1] wrt ice

      ! for subgrid variability of water

      real, dimension(klon, klev), intent(out) :: ratqs             ! ratqs i.e. factor for subgrid standard deviation of humidit
      real, dimension(klon, klev), intent(out) :: ratqs_inter_      ! interactive ratqs
      real, dimension(klon, klev), intent(out) :: sigma_qtherm      ! standard deviation of humidity in thermals [kg/kg]

      ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)

      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]

! for condensation and ice supersaturation

      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qsatliq        !--saturation specific humidity wrt liquid [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qsatice        !--saturation specific humidity wrt ice [kg/kg]

      ! for contrails and aviation

      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]

      ! for POPRECIP

      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]

      ! for thermals
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
      REAL, DIMENSION(klon, klev), INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment

      ! LOCAL variables
      !-------------------

      INTEGER                                      :: i, k
      REAL, DIMENSION(klon)                        :: rain_num
      REAL, DIMENSION(klon, klev)                  :: ql_estimate_lscp, qi_estimate_lscp
      REAL, DIMENSION(klon, klev)                  :: du0, dv0, dqbs0

      REAL, DIMENSION(klon) :: z_q, z_t, z_p, z_qs
      REAL, DIMENSION(klon) ::  z_qsl, z_qsi, z_dqs, z_dqsl, z_dqsi

!===============================================================================
! Computation of subgrid water variability
!===============================================================================

      !$gpum nocall
      CALL ratqs_main_first(klon, cell_area)

      CALL ratqs_main(klon, klev, nbsrf, is_ter, is_lic, dtime, &
                      pctsrf, s_pblh, zstd, wake_s, wake_deltaq, &
                      paprs, pplay, temp, qt, &
                      qtc_cv, sigt_cv, detrain_cv, fm_cv, fqd, fqcomp, sigd, qsat, &
                      fm_therm, entr_therm, detr_therm, cell_area, &
                      ratqsc, ratqs, ratqs_inter_, sigma_qtherm)

!===============================================================================
! Computation of large scale clouds and precipitation
!===============================================================================

      IF (ok_new_lscp) THEN

         ! recentmost version of the lscp routine. To use by default.

         DO k = 1, klev
            DO i = 1, klon
               ql_estimate_lscp(i, k) = ratio_ql_qtot(i, k)*qt(i, k)
               qi_estimate_lscp(i, k) = ratio_qi_qtot(i, k)*qt(i, k)
            END DO
         END DO

         CALL lscp(klon, klev, dtime, missing_val, paprs, pplay, omega, &
                   temp, qt, ql_estimate_lscp, qi_estimate_lscp, ptconv, ratqs, sigma_qtherm, &
                   d_t, d_q, d_ql, d_qi, rneb, rneblsvol, &
                   pfraclr, pfracld, cldfraliq, cldfraliqth, &
                   sigma2_icefracturb, sigma2_icefracturbth, &
                   mean_icefracturb, mean_icefracturbth, &
                   radocond, radicefrac, rain_lsc, snow_lsc, &
                   frac_impa, frac_nucl, beta, &
                   prfl, psfl, rhcl, &
                   qta, fraca, tv, pspsk, tla, thl, wth, &
                   distcltop, temp_cltop, &
                   tke, tke_dissip, &
                   entr_therm, detr_therm, &
                   cell_area, &
                   cf, rvc, u, v, &
                   qsub, qissr, qcld, subfra, issrfra, gamma_cond, &
                   dcf_sub, dcf_con, dcf_mix, dqi_adj, dqi_sub, dqi_con, dqi_mix, &
                   dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, qsatliq, qsatice, &
                   Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
                   dcf_avi, dqi_avi, dqvc_avi, flight_dist, flight_h2o, &
                   cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv, &
                   qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
                   dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)

      ELSE

         CALL fisrtilp(klon, klev, dtime, paprs, pplay, &
                       temp, qt, ptconv, ratqs, sigma_qtherm, &
                       d_t, d_q, d_ql, d_qi, rneb, rneblsvol, radocond, &
                       rain_lsc, snow_lsc, &
                       pfrac_impa, pfrac_nucl, pfrac_1nucl, &
                       frac_impa, frac_nucl, beta, &
                       prfl, psfl, rhcl, &
                       qta, fraca, tv, pspsk, tla, thl, &
                       cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)

      END IF

!===============================================================================
! Add tendencies
!===============================================================================
      du0(:, :) = 0.
      dv0(:, :) = 0.
      dqbs0(:, :) = 0.
      CALL add_phys_tend(du0, dv0, d_t, d_q, d_ql, d_qi, dqbs0, paprs, &
                         'lsc', abortphy, flag_inhib_tend, itap, 0)
      CALL prt_enerbil('lsc', itap)

!===============================================================================
! Final computations
!===============================================================================

      ! rain and snow are set to 0 when negative
      WHERE (rain_lsc < 0) rain_lsc = 0.
      WHERE (snow_lsc < 0) snow_lsc = 0.

      ! so-called 'numerical rain' is computed when qlnew=ql+dql>qlmax and qinew=qi+dqi>qimax
      ! i.e. when the condensed content is unrealistically large
      ! qlnew is set to qlmax and qinew to qimax

      rain_num(:) = 0.
      DO k = 1, klev
         DO i = 1, klon
            IF (ql_seri(i, k) > qlmax) THEN
               rain_num(i) = rain_num(i) + (ql_seri(i, k) - qlmax)*zmasse(i, k)/dtime
               ql_seri(i, k) = qlmax
            END IF
         END DO
      END DO
      IF (nqo >= 3) THEN
      DO k = 1, klev
         DO i = 1, klon
            IF (qs_seri(i, k) > qimax) THEN
               rain_num(i) = rain_num(i) + (qs_seri(i, k) - qimax)*zmasse(i, k)/dtime
               qs_seri(i, k) = qimax
            END IF
         END DO
      END DO
      END IF

      ! total precipitation
      DO i = 1, klon
         rain_fall(i) = rain_fall(i) + rain_lsc(i)
         snow_fall(i) = snow_fall(i) + snow_lsc(i)
      END DO

      ! relative humidity diagnostics
      DO k = 1, klev

         DO i=1, klon
            z_q(i) = q_seri(i, k)
            z_t(i) = t_seri(i, k)
            z_p(i) = pplay(i, k)
         END DO

         CALL CALC_QSAT_ECMWF(klon, z_t, z_q, z_p, RTT, 0, .false., z_qs, z_dqs)
         CALL CALC_QSAT_ECMWF(klon, z_t, z_q, z_p, RTT, 1, .false., z_qsl, z_dqsl)
         CALL CALC_QSAT_ECMWF(klon, z_t, z_q, z_p, RTT, 2, .false., z_qsi, z_dqsi)

         DO i = 1, klon
            rhl(i, k) = z_q(i)/z_qsl(i)
            rhi(i, k) = z_q(i)/z_qsi(i)
            rh(i, k) = z_q(i)/z_qs(i)
            qsat(i, k) = z_qs(i)
         END DO

      ENDDO



   END SUBROUTINE call_lscp

END MODULE lmdz_call_lscp
