Index: /LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90	(revision 4558)
+++ /LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90	(revision 4559)
@@ -12,6 +12,6 @@
   !$OMP THREADPRIVATE(seuil_neb)
 
-  INTEGER, SAVE :: ninter=5                     ! number of iterations to calculate autoconversion to precipitation
-  !$OMP THREADPRIVATE(ninter)
+  INTEGER, SAVE :: niter_lscp=5                 ! number of iterations to calculate autoconversion to precipitation
+  !$OMP THREADPRIVATE(niter_lscp)
 
   INTEGER,SAVE :: iflag_evap_prec=1             ! precipitation evaporation flag. 0: nothing, 1: "old way", 
@@ -41,7 +41,4 @@
   !$OMP THREADPRIVATE(ok_radocond_snow)
 
-  LOGICAL, SAVE :: ok_debug_autoconversion=.true.   ! removes a bug in the autoconversion process
-  !$OMP THREADPRIVATE(ok_debug_autoconversion)
-
   REAL, SAVE :: t_glace_min=258.0                ! lower-bound temperature parameter for cloud phase determination
   !$OMP THREADPRIVATE(t_glace_min)
@@ -77,4 +74,7 @@
   !$OMP THREADPRIVATE(iflag_pdf)
 
+  INTEGER, SAVE :: iflag_autoconversion=0        ! autoconversion option
+  !$OMP THREADPRIVATE(iflag_autoconversion)
+
   LOGICAL, SAVE :: reevap_ice=.false.            ! no liquid precip for T< threshold
   !$OMP THREADPRIVATE(reevap_ice)
@@ -92,4 +92,10 @@
   !$OMP THREADPRIVATE(cld_tau_con)
 
+  REAL, SAVE :: cld_expo_lsc=2.                  ! liquid autoconversion threshold exponent, stratiform rain
+  !$OMP THREADPRIVATE(cld_expo_lsc)
+
+  REAL, SAVE :: cld_expo_con=2.                  ! liquid autoconversion threshold exponent, convective rain
+  !$OMP THREADPRIVATE(cld_expo_con)
+
   REAL, SAVE :: ffallv_lsc=1.                    ! tuning coefficient crystal fall velocity, stratiform 
   !$OMP THREADPRIVATE(ffallv_lsc)
@@ -104,5 +110,9 @@
   !$OMP THREADPRIVATE(coef_eva_i)
 
-
+  REAL cice_velo                                 ! factor in the ice fall velocity formulation
+  PARAMETER (cice_velo=1.645)
+
+  REAL dice_velo                                 ! exponent in the ice fall velocity formulation
+  PARAMETER (dice_velo=0.16)
 
 
@@ -122,4 +132,6 @@
    REAL, INTENT(IN)      :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in
    REAL, INTENT(IN)      ::  RVTMP2_in, RTT_in, RD_in, RG_in
+   character (len=20) :: modname='lscp_ini_mod'
+   character (len=80) :: abort_message
 
 
@@ -136,5 +148,5 @@
 
 
-    CALL getin_p('ninter',ninter)
+    CALL getin_p('niter_lscp',niter_lscp)
     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
     CALL getin_p('seuil_neb',seuil_neb)
@@ -142,5 +154,4 @@
     CALL getin_p('iflag_mpc_bl',iflag_mpc_bl)
     CALL getin_p('ok_radocond_snow',ok_radocond_snow) 
-    CALL getin_p('ok_debug_autoconversion',ok_debug_autoconversion)    
     CALL getin_p('t_glace_max',t_glace_max)
     CALL getin_p('t_glace_min',t_glace_min)
@@ -159,4 +170,6 @@
     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)
@@ -164,9 +177,10 @@
     coef_eva_i=coef_eva
     CALL getin_p('coef_eva_i',coef_eva_i)
-
-
-
-
-    WRITE(lunout,*) 'lscp, ninter:', ninter
+    CALL getin_p('iflag_autoconversion',iflag_autoconversion)
+
+
+
+
+    WRITE(lunout,*) 'lscp, niter_lscp:', niter_lscp
     WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec
     WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb
@@ -174,5 +188,4 @@
     WRITE(lunout,*) 'lscp, iflag_mpc_bl:', iflag_mpc_bl
     WRITE(lunout,*) 'lscp, ok_radocond_snow:', ok_radocond_snow
-    WRITE(lunout,*) 'lscp, ok_debug_autoconversion:', ok_debug_autoconversion
     WRITE(lunout,*) 'lscp, t_glace_max:', t_glace_max
     WRITE(lunout,*) 'lscp, t_glace_min:', t_glace_min
@@ -191,8 +204,11 @@
     WRITE(lunout,*) 'lscp, cld_tau_lsc', cld_tau_lsc
     WRITE(lunout,*) 'lscp, cld_tau_con', cld_tau_con
+    WRITE(lunout,*) 'lscp, cld_expo_lsc', cld_expo_lsc
+    WRITE(lunout,*) 'lscp, cld_expo_con', cld_expo_con
     WRITE(lunout,*) 'lscp, ffallv_lsc', ffallv_lsc
     WRITE(lunout,*) 'lscp, ffallv_con', ffallv_con
     WRITE(lunout,*) 'lscp, coef_eva', coef_eva
     WRITE(lunout,*) 'lscp, coef_eva_i', coef_eva_i
+    WRITE(lunout,*) 'lscp, iflag_autoconversion', iflag_autoconversion
 
 
@@ -201,7 +217,17 @@
 
     ! check for precipitation sub-time steps
-    IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
+    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'
+    ENDIF
+
+    ! check consistency between numerical resolution of autoconversion 
+    ! and other options
+   
+    IF (iflag_autoconversion .EQ. 2) 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) 
+        ENDIF 
     ENDIF
 
Index: /LMDZ6/trunk/libf/phylmd/lscp_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4558)
+++ /LMDZ6/trunk/libf/phylmd/lscp_mod.F90	(revision 4559)
@@ -93,9 +93,10 @@
 USE ice_sursat_mod, ONLY : ice_sursat
 
-USE lscp_ini_mod, ONLY : seuil_neb, ninter, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
-USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, ok_debug_autoconversion
+USE lscp_ini_mod, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
+USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
 USE lscp_ini_mod, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min 
 USE lscp_ini_mod, ONLY : coef_eva, coef_eva_i,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
-USE lscp_ini_mod, ONLY : iflag_bergeron, iflag_fisrtilp_qsat
+USE lscp_ini_mod, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
+USE lscp_ini_mod, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc
 USE lscp_ini_mod, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
 
@@ -190,5 +191,5 @@
 
   REAL qsl(klon), qsi(klon)
-  REAL zct, zcl
+  REAL zct, zcl,zexpo
   REAL ctot(klon,klev)
   REAL ctot_vol(klon,klev)
@@ -232,7 +233,8 @@
   REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
   REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
-  REAL velo(klon,klev), vr
+  REAL velo(klon,klev), vr, ffallv
   REAL qlmpc, qimpc, rnebmpc
   REAL radocondi(klon,klev), radocondl(klon,klev)
+  REAL effective_zneb
 
   INTEGER i, k, n, kk
@@ -922,4 +924,9 @@
     ENDIF
 
+
+    ! Autoconversion
+    ! -------------------------------------------------------------------------------
+
+
     ! Initialisation of zoliq and radocond variables
 
@@ -933,14 +940,12 @@
             iwc(i)   = 0. 
             zneb(i)  = MAX(rneb(i,k),seuil_neb)
-            radocond(i,k) = zoliq(i)/REAL(ninter+1)
+            radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
             radicefrac(i,k) = zfice(i)
-            radocondi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
-            radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
+            radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
+            radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
     ENDDO
 
-    ! Autoconversion
-    ! -------------------------------------------------------------------------------
-     
-    DO n = 1, ninter
+       
+    DO n = 1, niter_lscp
 
         DO i=1,klon
@@ -968,18 +973,16 @@
                 ELSE
 
-                    IF (ok_debug_autoconversion) THEN ! if condition to be removed after test phase
-
-                    ! water quantity to remove: zchau (Sundqvist, 1978)
-                    ! same thing for the ice: zfroi (Zender & Kiehl, 1997)
                     IF (ptconv(i,k)) THEN ! if convective point
                         zcl=cld_lc_con
                         zct=1./cld_tau_con
-                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i)*velo(i,k)
+                        zexpo=cld_expo_con
+                        ffallv=ffallv_con
                     ELSE
                         zcl=cld_lc_lsc
                         zct=1./cld_tau_lsc
-                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i) &   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
-                            *velo(i,k)
+                        zexpo=cld_expo_lsc
+                        ffallv=ffallv_lsc
                     ENDIF
+
 
                     ! if vertical heterogeneity is taken into account, we use
@@ -988,26 +991,8 @@
                     ! reduces the in-cloud water).
 
+                    ! Liquid water quantity to remove: zchau (Sundqvist, 1978)
+                    ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
+                    !.........................................................
                     IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
-                        zchau = zct   *dtime/REAL(ninter) * zoliql(i) &
-                        *(1.0-EXP(-(zoliql(i)/ctot_vol(i,k)/zcl)**2))
-                    ELSE
-                        zchau = zct   *dtime/REAL(ninter) * zoliql(i) &
-                        *(1.0-EXP(-(zoliql(i)/zneb(i)/zcl)**2))        ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
-                    ENDIF
-
-                    ELSE ! with old bug in autoconversion 
-
-                    ! water quantity to remove: zchau (Sundqvist, 1978)
-                    ! same thing for the ice: zfroi (Zender & Kiehl, 1997)
-                    IF (ptconv(i,k)) THEN ! if convective point
-                        zcl=cld_lc_con
-                        zct=1./cld_tau_con
-                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*zfice(i)
-                    ELSE
-                        zcl=cld_lc_lsc
-                        zct=1./cld_tau_lsc
-                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
-                            *velo(i,k) * zfice(i)
-                    ENDIF
 
                     ! if vertical heterogeneity is taken into account, we use
@@ -1015,15 +1000,39 @@
                     ! surface fraction (which is larger and artificially
                     ! reduces the in-cloud water).
-
-
-                    IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
-                        zchau = zct   *dtime/REAL(ninter) * zoliq(i) &
-                        *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl)**2)) *(1.-zfice(i))
+                        effective_zneb=ctot_vol(i,k)
                     ELSE
-                        zchau = zct   *dtime/REAL(ninter) * zoliq(i) &
-                        *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl)**2)) *(1.-zfice(i)) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
+                        effective_zneb=zneb(i)
+                    ENDIF 
+
+
+                    IF (iflag_autoconversion .EQ. 2) THEN
+                    ! two-steps resolution with niter_lscp=1 sufficient
+                    ! we first treat the second term (with exponential) in an explicit way
+                    ! and then treat the first term (-q/tau) in an exact way
+                        zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct &
+                        *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo))))
+                    ELSE
+                    ! old explicit resolution with subtimesteps
+                        zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) &
+                        *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo))
                     ENDIF
 
-                    ENDIF ! ok_debug_autoconversion
+
+                    ! Ice water quantity to remove (Zender & Kiehl, 1997)
+                    ! dqice/dt=1/rho*d(rho*wice*qice)/dz
+                    !....................................
+                    IF (iflag_autoconversion .EQ. 2) THEN
+                    ! exact resolution, niter_lscp=1 is sufficient but works only
+                    ! with iflag_vice=0
+                       IF (zoliqi(i) .GT. 0.) THEN
+                          zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & 
+                          +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo))
+                       ELSE
+                          zfroi=0.
+                       ENDIF
+                    ELSE
+                    ! old explicit resolution with subtimesteps
+                       zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k) 
+                    ENDIF
 
                     zrain   = MIN(MAX(zchau,0.0),zoliql(i))
@@ -1033,13 +1042,21 @@
                 ENDIF
 
-                zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
-                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
-                zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
+                IF (iflag_autoconversion .GE. 1) THEN
+                   ! debugged version with phase conservation through the autoconversion process
+                   zoliql(i) = MAX(zoliql(i)-1.*zrain  , 0.0)
+                   zoliqi(i) = MAX(zoliqi(i)-1.*zsnow  , 0.0)
+                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
+                ELSE
+                   ! bugged version with phase resetting
+                   zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
+                   zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
+                   zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
+                ENDIF
 
                 ! c_iso: call isotope_conversion (for readibility) or duplicate 
 
-                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(ninter+1)
-                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(ninter+1)
-                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(ninter+1)
+                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
+                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
+                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)
 
             ENDIF ! rneb >0
@@ -1137,4 +1154,5 @@
         ENDDO
 
+
     ENDIF
 
Index: /LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90
===================================================================
--- /LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90	(revision 4558)
+++ /LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90	(revision 4559)
@@ -17,4 +17,5 @@
     
     use lscp_ini_mod, only: iflag_vice, ffallv_con, ffallv_lsc
+    use lscp_ini_mod, only: cice_velo, dice_velo 
 
     IMPLICIT NONE
@@ -31,5 +32,5 @@
 
     INTEGER i
-    REAL logvm,iwcg,tempc,phpa,cvel,dvel,fallv_tun
+    REAL logvm,iwcg,tempc,phpa,fallv_tun
     REAL m2ice, m2snow, vmice, vmsnow
     REAL aice, bice, asnow, bsnow
@@ -60,6 +61,4 @@
         
         velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s
-        dvel=0.2
-        cvel=fallv_tun*65.0*(rho(i)**0.2)*(150./phpa)**0.15
 
     ELSE IF (iflag_vice .EQ. 2) THEN
@@ -92,12 +91,8 @@
         vmsnow=vmsnow*((1000./phpa)**0.35)
         velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s
-        dvel=0.2
-        cvel=velo(i)/((iwcg/1000.*rho(i))**dvel)
 
     ELSE
         ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990
-        velo(i) = fallv_tun*3.29/2.0*((iwcg/1000.)**0.16)
-        dvel=0.16
-        cvel=fallv_tun*3.29/2.0*(rho(i)**0.16)
+        velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo)
     ENDIF
     ENDDO
