Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/atke_exchange_coeff_mod.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/atke_exchange_coeff_mod.F90	(revision 4585)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/atke_exchange_coeff_mod.F90	(revision 4585)
@@ -0,0 +1,199 @@
+module atke_exchange_coeff_mod
+
+implicit none
+
+contains 
+
+subroutine atke_compute_km_kh(ngrid,nlay, &
+                        wind_u,wind_v,temp,play,pinterf, &
+                        tke,Km_out,Kh_out)
+
+!========================================================================
+! Routine that computes turbulent Km / Kh coefficients with a
+! 1.5 order closure scheme (TKE) with or without stationarity assumption
+!
+! This parameterization has been constructed in the framework of a
+! collective and collaborative workshop, 
+! the so-called 'Atelier TKE (ATKE)' with
+! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos, 
+! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon
+!
+! Main assumptions of the model :
+! (1) dry atmosphere
+! (2) horizontal homogeneity (Dx=Dy=0.)
+!=======================================================================
+
+
+
+USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd
+USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd
+USE atke_turbulence_ini_mod, ONLY : viscom, viscoh
+
+implicit none
+
+
+! Declarations:
+!=============
+
+INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
+INTEGER, INTENT(IN) :: nlay ! number of vertical index  
+
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
+
+
+REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
+
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers
+
+! Local variables
+REAL, DIMENSION(ngrid,nlay+1) :: Km          ! Exchange coefficient for momentum at interface between layers
+REAL, DIMENSION(ngrid,nlay+1) :: Kh          ! Exchange coefficient for heat flux at interface between layers
+REAL, DIMENSION(ngrid,nlay)   :: theta       ! Potential temperature
+REAL, DIMENSION(ngrid,nlay+1) :: l_exchange  ! Length of exchange (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: z_interf    ! Altitude at the interface
+REAL, DIMENSION(ngrid,nlay)   :: z_lay       ! Altitude of layers
+REAL, DIMENSION(ngrid,nlay)   :: dz_interf   ! distance between two consecutive interfaces
+REAL, DIMENSION(ngrid,nlay)   :: dz_lay      ! distance between two layer middles (NB: first and last are half layers)
+REAL, DIMENSION(ngrid,nlay+1) :: Ri          ! Richardson's number (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Prandtl     ! Turbulent Prandtl's number (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Sm          ! Stability function for momentum (at interface)
+REAL, DIMENSION(ngrid,nlay+1) :: Sh          ! Stability function for heat (at interface)
+
+INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
+REAL    :: cn,Ri0,Ri1    ! parameter for Sm stability function and Prandlt
+REAL    :: preff      ! reference pressure for potential temperature calculations
+REAL    :: thetam     ! mean potential temperature at interface
+
+
+! Initializations:
+!================
+
+DO igrid=1,ngrid
+    dz_interf(igrid,1) = 0.0
+    z_interf(igrid,1) = 0.0
+END DO
+
+! Calculation of potential temperature: (if vapor -> todo virtual potential temperature)
+!=====================================
+
+preff=100000.
+! The result should not depend on the choice of preff
+DO ilay=1,nlay
+     DO igrid = 1, ngrid
+        theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd)
+     END DO
+END DO
+
+
+
+! Calculation of altitude of layers' middle and bottom interfaces:
+!=================================================================
+
+DO ilay=2,nlay+1
+    DO igrid=1,ngrid
+        dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
+        z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
+    ENDDO
+ENDDO
+
+DO ilay=1,nlay
+   DO igrid=1,ngrid 
+      z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
+   ENDDO
+ENDDO
+
+
+! Computing the mixing length:
+! so far, we have neglected the effect of local stratification
+!==============================================================
+
+
+DO ilay=2,nlay+1
+    DO igrid=1,ngrid
+        l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
+    ENDDO
+ENDDO
+
+! Computes the gradient Richardson's number and stability functions:
+!===================================================================
+
+! calculation of cn = Sm value at Ri=0
+! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0
+cn=(1./sqrt(cepsilon))**(2/3)
+! calculation of Ri0 such that continuity in slope of Sm at Ri=0
+Ri0=2./rpi*(cinf - cn)*ric/cn
+! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0
+Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope
+
+
+DO ilay=2,nlay
+    DO igrid=1,ngrid
+        dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
+        thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
+        Ri(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) / & 
+        MAX(((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 
+        ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2,1E-10) 
+        
+        IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
+            Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
+            Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut 
+        ELSE ! stable cases
+            Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric))
+            Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope
+            IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN
+               call abort_physic("atke_compute_km_kh", &
+               'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
+            ENDIF
+        END IF
+        
+        Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
+
+    ENDDO
+ENDDO
+
+! Computing the TKE:
+!===================
+IF (iflag_atke == 0) THEN
+
+! stationary solution neglecting the vertical transport of TKE by turbulence
+    DO ilay=2,nlay
+        DO igrid=1,ngrid
+            tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * &
+            (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & 
+            ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) * &
+            (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
+        ENDDO
+    ENDDO
+
+ELSE ! TO DO
+    
+   call abort_physic("atke_compute_km_kh", &
+        'traitement non-stationnaire de la tke pas encore prevu', 1)
+
+END IF
+
+
+! Computing eddy diffusivity coefficients:
+!========================================
+DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
+    DO igrid=1,ngrid
+        ! we add the molecular viscosity to Km,h 
+        Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
+        Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
+    END DO
+END DO
+
+! for output:
+!===========
+Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
+Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
+
+end subroutine atke_compute_km_kh
+
+
+end module atke_exchange_coeff_mod
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/atke_turbulence_ini_mod.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/atke_turbulence_ini_mod.F90	(revision 4585)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/atke_turbulence_ini_mod.F90	(revision 4585)
@@ -0,0 +1,84 @@
+MODULE atke_turbulence_ini_mod
+
+implicit none
+
+save
+
+integer :: iflag_atke, iflag_num_atke
+!$OMP THREADPRIVATE(iflag_atke, iflag_num_atke)
+  real :: kappa = 0.4 ! Von Karman constant
+  !$OMP THREADPRIVATE(kappa)
+  real :: l0, ric, ri0, cinf, cepsilon, pr_slope, pr_asym, pr_neut
+  !$OMP THREADPRIVATE(l0, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut)
+  integer :: lunout,prt_level
+  !$OMP THREADPRIVATE(lunout,prt_level)
+  real :: rg, rd, rpi, rcpd
+  !$OMP THREADPRIVATE(rg, rd, rpi, rcpd)
+
+  real :: viscom, viscoh
+  !$OMP THREADPRIVATE(viscom,viscoh)
+
+
+
+CONTAINS
+
+SUBROUTINE atke_ini(prt_level_in, lunout_in, rg_in, rd_in, rpi_in, rcpd_in)
+
+   USE ioipsl_getin_p_mod, ONLY : getin_p
+
+  integer, intent(in) :: lunout_in,prt_level_in
+  real, intent(in) :: rg_in, rd_in, rpi_in, rcpd_in
+
+
+  lunout=lunout_in
+  prt_level=prt_level_in
+  rd=rd_in
+  rg=rg_in
+  rpi=rpi_in
+  rcpd=rcpd_in
+
+  viscom=1.46E-5
+  viscoh=2.06E-5
+
+  ! flag that controls options in atke_compute_km_kh
+  iflag_atke=0
+  CALL getin_p('iflag_atke',iflag_atke)
+
+  ! flag that controls the numerical treatment of diffusion coeffiient calculation
+  iflag_num_atke=0
+  CALL getin_p('iflag_num_atke',iflag_num_atke)
+
+  ! asymptotic mixing length [m]
+  l0=150.0 
+  CALL getin_p('atke_l0',l0)
+
+  ! critical Richardson number 
+  ric=0.25
+  CALL getin_p('atke_ric',ric)
+
+  ! asymptotic value of Sm for Ri=-Inf
+  cinf=1.5
+  CALL getin_p('atke_cinf',cinf)
+
+  ! constant for tke dissipation calculation
+  cepsilon=16.6 ! default value as in yamada4
+  CALL getin_p('atke_cepsilon',cepsilon)
+
+  ! slope of Pr=f(Ri) for stable conditions
+  pr_slope=5.0 ! default value from Zilitinkevich et al. 2005
+  CALL getin_p('atke_pr_slope',pr_slope)
+
+  ! asymptotic turbulent prandt number value for Ri=-Inf
+  pr_asym=0.4
+  CALL getin_p('atke_pr_asym',pr_asym)
+
+  ! value of turbulent prandtl number in neutral conditions (Ri=0)
+  pr_neut=0.8
+  CALL getin_p('atke_pr_neut',pr_neut)
+
+   
+ RETURN
+
+END SUBROUTINE atke_ini
+
+END MODULE  atke_turbulence_ini_mod
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/blowing_snow_ini_mod.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/blowing_snow_ini_mod.F90	(revision 4585)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/blowing_snow_ini_mod.F90	(revision 4585)
@@ -0,0 +1,75 @@
+module blowing_snow_ini_mod
+
+implicit none
+
+save
+   integer :: prt_level=0,lunout
+   real RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG
+   real coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs
+   real prt_bs, pbst_bs, qbst_bs
+   integer :: niter_bs, iflag_saltation_bs
+
+   !$OMP THREADPRIVATE(prt_level, lunout)
+   !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)
+   !$OMP THREADPRIVATE(coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs)
+   !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs)
+   !$OMP THREADPRIVATE(niter_bs, iflag_saltation_bs)
+
+      contains
+
+      subroutine blowing_snow_ini(prt_level_in,lunout_in, &
+                                  RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in,&
+                                  RVTMP2_in, RTT_in,RD_in,RG_in)
+
+         USE ioipsl_getin_p_mod, ONLY : getin_p
+
+         integer, intent(in) :: prt_level_in,lunout_in
+         real, intent(in) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in
+         real, intent(in) ::  RVTMP2_in, RTT_in, RD_in, RG_in
+
+
+         prt_level=prt_level_in
+         RG=RG_in
+         RD=RD_in
+         RCPD=RCPD_in
+         RLVTT=RLVTT_in
+         RLSTT=RLSTT_in
+         RLMLT=RLMLT_in
+         RTT=RTT_in
+         RG=RG_in
+         RVTMP2=RVTMP2_in
+         lunout=lunout_in
+
+
+         qbst_bs= 0.001
+         CALL getin_p('qbst_bs',qbst_bs)
+
+         pbst_bs= 0.01
+         CALL getin_p('pbst_bs',pbst_bs)
+
+         prt_bs= 0.0003
+         CALL getin_p('prt_bs',prt_bs)
+
+
+         zeta_bs= 3.
+         CALL getin_p('zeta_bs',zeta_bs)
+
+         fallv_bs = 0.1
+         CALL getin_p('fallv_bs',fallv_bs)
+
+         coef_eva_bs =  2.e-5 
+         CALL getin_p('coef_eva_bs',coef_eva_bs)
+
+         expo_eva_bs = 0.5
+         CALL getin_p('expo_eva_bs',expo_eva_bs)
+
+         niter_bs = 5
+         CALL getin_p('niter_bs',niter_bs)
+
+         iflag_saltation_bs=0
+         CALL getin_p('iflag_saltation_bs',iflag_saltation_bs)
+
+
+      end subroutine blowing_snow_ini
+
+end module blowing_snow_ini_mod
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/blowing_snow_sublim_sedim.F90	(revision 4585)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/blowing_snow_sublim_sedim.F90	(revision 4585)
@@ -0,0 +1,212 @@
+subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
+
+!==============================================================================
+! Routine that calculates the evaporation and sedimentation of blowing snow
+! inspired by what is done in lscp_mod
+! Etienne Vignon, October 2022
+!==============================================================================
+
+
+use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs
+use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, niter_bs
+USE lscp_tools_mod, only : calc_qsat_ecmwf
+
+implicit none
+
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Declarations
+!++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+!INPUT
+!=====
+
+integer, intent(in)                     :: ngrid,nlay
+real, intent(in)                        :: dtime
+real, intent(in), dimension(ngrid,nlay) :: temp
+real, intent(in), dimension(ngrid,nlay) :: q
+real, intent(in), dimension(ngrid,nlay) :: qbs
+real, intent(in), dimension(ngrid,nlay) :: pplay
+real, intent(in), dimension(ngrid,nlay+1) :: paprs
+
+
+
+! OUTPUT
+!========
+
+
+real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
+real, intent(out), dimension(ngrid,nlay) :: dq_bs
+real, intent(out), dimension(ngrid,nlay) :: dqbs_bs
+real, intent(out), dimension(ngrid,nlay+1) :: bsfl
+real, intent(out), dimension(ngrid)      :: precip_bs
+
+
+! LOCAL
+!======
+
+
+integer                                  :: k,i,n
+real                                     :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt
+real                                     :: dqsedim,precbs
+real, dimension(ngrid)                   :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
+real, dimension(ngrid)                   :: zqbsi,zmqc, zmair, zdz
+real, dimension(ngrid,nlay)              :: velo, zrho
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++
+! Initialisation
+!++++++++++++++++++++++++++++++++++++++++++++++++++
+
+qzero(:)=0.
+dtemp_bs(:,:)=0.
+dq_bs(:,:)=0.
+dqbs_bs(:,:)=0.
+velo(:,:)=0.
+zt(:)=0.
+zq(:)=0.
+zqbs(:)=0.
+sedim(:)=0.
+
+! begin of top-down loop
+DO k = nlay, 1, -1
+    
+    DO i=1,ngrid
+        zt(i)=temp(i,k)
+        zq(i)=q(i,k)
+        zqbs(i)=qbs(i,k)
+    ENDDO
+
+
+    IF (k.LE.nlay-1) THEN
+
+        ! thermalization of blowing snow precip coming from above    
+        DO i = 1, ngrid
+            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+            ! RVTMP2=rcpv/rcpd-1
+            zcpair=RCPD*(1.0+RVTMP2*zq(i))
+            zcpeau=RCPD*RVTMP2
+            ! zmqc: precipitation mass that has to be thermalized with 
+            ! layer's air so that precipitation at the ground has the
+            ! same temperature as the lowermost layer
+            zmqc(i) = (sedim(i))*dtime/zmair(i)
+            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
+            zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
+             / (zcpair + zmqc(i)*zcpeau)
+        ENDDO
+    ELSE
+        DO i = 1, ngrid
+            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
+            zmqc(i) = 0.
+        ENDDO
+
+    ENDIF
+
+
+
+    ! calulation saturation specific humidity
+    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
+    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
+    ! sublimation calculation
+    ! SUndqvist formula dP/dz=beta*(1-q/qsat)*sqrt(P)
+    
+    DO i = 1, ngrid
+       ! if sedimentation:
+       IF (sedim(i) .GT. 0.) THEN
+
+          IF (zt(i) .GT. RTT) THEN
+             ! if positive celcius temperature, we assume
+             ! that all the blowing snow melt and evaporate
+              zqevti=sedim(i)*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG     
+             ! we ensure that the whole mesh does not exceed saturation wrt liquid
+              zqev0 = MAX(0.0, qsl(i)-zq(i))
+              zqevi = MIN(zqev0,zqevti)
+              ! New solid precipitation fluxes
+              sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
+              ! vapor, temperature, precip fluxes update
+              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              ! melting
+              zt(i) = zt(i)                             &
+                + (sedimn(i)-sedim(i))                      &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
+                * RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+              ! evaporation
+              zt(i) = zt(i)                             &
+                + (sedimn(i)-sedim(i))                      &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
+                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))      
+      
+      
+          ELSE     
+              zqevti = coef_eva_bs*(1.0-zq(i)/qsi(i))*(sedim(i)**expo_eva_bs) &
+                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
+              zqevti = MAX(0.0,MIN(zqevti,sedim(i)))*RG*dtime/(paprs(i,k)-paprs(i,k+1))
+
+              ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
+              zqev0 = MAX(0.0, qsi(i)-zq(i))
+              zqevi = MIN(zqev0,zqevti)
+
+              ! New solid precipitation fluxes
+              sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
+
+
+              ! vapor, temperature, precip fluxes update
+              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime                
+              zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
+              zt(i) = zt(i)                             &
+                + (sedimn(i)-sedim(i))                      &
+                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
+                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
+           ENDIF
+
+
+          ! sedim update
+          sedim(i)=sedimn(i)
+
+
+        ELSE
+           sedim(i)=0.
+        ENDIF ! if sedim > 0
+
+        zqbsi(i)=zqbs(i)
+
+    ENDDO ! loop on ngrid
+
+    ! Now sedimention scheme (alike that in lscp)
+    DO n = 1, niter_bs
+       DO i = 1, ngrid
+         zrho(i,k)  = pplay(i,k) / zt(i) / RD
+         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
+         velo(i,k) = fallv_bs 
+         dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k)   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
+         precbs   = MIN(MAX(dqsedim,0.0),zqbs(i))
+         zqbs(i) = MAX(zqbs(i)-1.*precbs  , 0.0)
+       ENDDO !loop on ngrid
+    ENDDO ! loop on niter_bs
+
+    ! add to non sublimated precip
+    DO i=1,ngrid
+    sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
+    ENDDO
+
+    ! Outputs:
+    DO i = 1, ngrid
+        bsfl(i,k)=sedim(i)
+        dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
+        dq_bs(i,k) = zq(i) - q(i,k)
+        dtemp_bs(i,k) = zt(i) - temp(i,k)
+    ENDDO
+
+
+ENDDO ! vertical loop
+
+
+!surface bs flux
+DO i = 1, ngrid
+  precip_bs(i) = sedim(i)
+ENDDO
+
+
+return
+
+end subroutine blowing_snow_sublim_sedim
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/call_atke_mod.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/call_atke_mod.F90	(revision 4585)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/call_atke_mod.F90	(revision 4585)
@@ -0,0 +1,160 @@
+module call_atke_mod
+
+USE atke_exchange_coeff_mod, ONLY :  atke_compute_km_kh
+
+implicit none
+
+
+contains
+
+subroutine call_atke(dtime,ngrid,nlay,cdrag_uv,cdrag_t,u_surf,v_surf,temp_surf, &
+                        wind_u,wind_v,temp,play,pinterf, &
+                        tke,Km_out,Kh_out)
+
+
+
+
+USE atke_turbulence_ini_mod, ONLY : iflag_num_atke, rg, rd
+
+implicit none
+
+
+! Declarations:
+!=============
+
+REAL, INTENT(IN)    :: dtime ! timestep (s)
+INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
+INTEGER, INTENT(IN) :: nlay ! number of vertical index  
+
+
+REAL, DIMENSION(ngrid), INTENT(IN)       :: cdrag_uv  ! drag coefficient for wind
+REAL, DIMENSION(ngrid), INTENT(IN)       :: cdrag_t   ! drag coefficient for temperature
+
+REAL, DIMENSION(ngrid), INTENT(IN)       :: u_surf    ! x wind velocity at the surface
+REAL, DIMENSION(ngrid), INTENT(IN)       :: v_surf    ! y wind velocity at the surface
+REAL, DIMENSION(ngrid), INTENT(IN)       :: temp_surf ! surface temperature
+
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_u   ! zonal velocity (m/s)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: wind_v   ! meridional velocity (m/s)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
+REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
+
+
+REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
+
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Km_out   ! output: Exchange coefficient for momentum at interface between layers
+REAL, DIMENSION(ngrid,nlay), INTENT(OUT)      :: Kh_out   ! output: Exchange coefficient for heat flux at interface between layers
+
+
+REAL, DIMENSION(ngrid,nlay) :: wind_u_predict, wind_v_predict
+REAL, DIMENSION(ngrid) ::  wind1
+INTEGER i
+
+
+
+call atke_compute_km_kh(ngrid,nlay, &
+                        wind_u,wind_v,temp,play,pinterf, &
+                        tke,Km_out,Kh_out)
+
+if (iflag_num_atke .EQ. 1) then
+
+   !! pay attention that the treatment of the TKE
+   !! has to be adapted when solving the TKE with a prognostic equation
+
+   do i=1,ngrid
+      wind1(i)=sqrt(wind_u(i,1)**2+wind_v(i,1)**2)
+   enddo
+   call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_u,Km_out,u_surf,cdrag_uv,wind_u_predict)  
+   call atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,wind_v,Km_out,v_surf,cdrag_uv,wind_v_predict)
+
+
+   call atke_compute_km_kh(ngrid,nlay, &
+                        wind_u_predict,wind_v_predict,temp,play,pinterf, &
+                        tke,Km_out,Kh_out)
+
+end if
+
+end subroutine call_atke
+
+!----------------------------------------------------------------------------------------
+
+subroutine atke_explicit_prediction(ngrid,nlay,rg,rd,dtime,pinterf,play,temp,wind1,x_in,K_in,x_surf,cdrag,x_predict)
+
+INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
+INTEGER, INTENT(IN) :: nlay ! number of vertical index  
+REAL, INTENT(IN) :: rg,rd,dtime ! gravity, R dry air and timestep
+REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: play      ! pressure middle of layers (Pa)
+REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: temp      ! temperature (K)
+REAL, DIMENSION(ngrid),        INTENT(IN)  :: wind1     ! wind speed first level (m/s)
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)  :: pinterf ! pressure at interfaces(Pa)
+REAL, DIMENSION(ngrid,nlay),   INTENT(IN)  :: x_in    ! variable at the beginning of timestep
+REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)  :: K_in    ! eddy diffusivity coef at the beginning of time step
+REAL, DIMENSION(ngrid),        INTENT(IN)  :: x_surf  ! surface variable
+REAL, DIMENSION(ngrid),        INTENT(IN)  :: cdrag  ! drag coefficient
+REAL, DIMENSION(ngrid,nlay),   INTENT(OUT) :: x_predict  ! variable at the end of time step after explicit prediction
+
+
+integer i,k
+real ml,F1,rho
+real, dimension(ngrid) :: play1,temp1
+real, dimension(ngrid,nlay+1) :: K_big
+
+! computation of K_big
+
+play1(:)=play(:,1)
+temp1(:)=temp(:,1)
+
+! "big K" calculation
+do  k=2,nlay-1
+   do i=1,ngrid
+      rho=pinterf(i,k)/rd/(0.5*(temp(i,k-1)+temp(i,k)))
+      K_big(i,k)=rg*K_in(i,k)/(play(i,k)-play(i,k+1))*(rho**2)
+   enddo
+enddo
+! speficic treatment for k=nlay
+do i=1,ngrid
+   rho=pinterf(i,nlay)/rd/temp(i,nlay)
+   K_big(i,nlay)=rg*K_in(i,nlay)/(2*(play(i,nlay)-pinterf(i,nlay+1)))*(rho**2)
+enddo
+
+
+
+! x_predict calculation for 2<=k<=nlay-1
+do  k=2,nlay-1
+   do i=1,ngrid
+      ml=(pinterf(i,k)-pinterf(i,k+1))/rg
+      x_predict(i,k)=x_in(i,k)-dtime/ml*(-K_big(i,k+1)*x_in(i,k+1) &
+                  + (K_big(i,k)+K_big(i,k+1))*x_in(i,k) &
+                  - K_big(i,k)*x_in(i,k-1))
+   enddo
+enddo
+
+! Specific treatment for k=1 
+do i=1,ngrid
+   ml=(pinterf(i,1)-pinterf(i,2))/rg
+   F1=-play1(i)/rd/temp1(i)*wind1(i)*cdrag(i)*(x_in(i,1)-x_surf(i)) ! attention convention sens du flux
+   x_predict(i,1)=x_in(i,1)-dtime/ml*(-K_big(i,2)*(x_in(i,2) - x_in(i,1))-F1)
+enddo
+
+! Specific treatment for k=nlay
+! flux at the top of the atmosphere=0
+do i=1,ngrid
+   ml=0.5*(pinterf(i,nlay)-pinterf(i,nlay+1))/rg
+   x_predict(i,nlay)=x_in(i,nlay)+dtime/ml*(K_big(i,nlay)*(x_in(i,nlay)-x_in(i,nlay-1)))
+enddo
+
+!K_big(:,1)=0.
+!do  k=1,nlay
+!   do i=1,ngrid
+!      print*, 'youhou', k, x_in(i,k), x_predict(i,k), K_big(i,k)
+!   end do
+!enddo
+
+
+
+end subroutine atke_explicit_prediction
+
+
+
+end module call_atke_mod
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90	(revision 4585)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/call_blowing_snow_sublim_sedim.F90	(revision 4585)
@@ -0,0 +1,38 @@
+subroutine call_blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs, &
+                                  dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
+
+implicit none
+
+!INPUT
+!=====
+
+integer, intent(in)                     :: ngrid,nlay
+real, intent(in)                        :: dtime
+real, intent(in), dimension(ngrid,nlay) :: temp
+real, intent(in), dimension(ngrid,nlay) :: q
+real, intent(in), dimension(ngrid,nlay) :: qbs
+real, intent(in), dimension(ngrid,nlay) :: pplay
+real, intent(in), dimension(ngrid,nlay+1) :: paprs
+
+
+! OUTPUT
+!========
+
+
+real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
+real, intent(out), dimension(ngrid,nlay) :: dq_bs
+real, intent(out), dimension(ngrid,nlay) :: dqbs_bs
+real, intent(out), dimension(ngrid,nlay+1) :: bsfl
+real, intent(out), dimension(ngrid)      :: precip_bs
+
+      
+
+
+call blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs, &
+                                  dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
+
+
+
+
+
+end subroutine call_blowing_snow_sublim_sedim
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/climb_qbs_mod.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/climb_qbs_mod.F90	(revision 4585)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/climb_qbs_mod.F90	(revision 4585)
@@ -0,0 +1,389 @@
+MODULE climb_qbs_mod
+!
+! Module to solve the verctical diffusion of blowing snow; 
+!
+  USE dimphy
+
+  IMPLICIT NONE
+  SAVE 
+  PRIVATE
+  PUBLIC :: climb_qbs_down, climb_qbs_up
+
+  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaqbs
+  !$OMP THREADPRIVATE(gamaqbs)
+  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_QBS, Dcoef_QBS
+  !$OMP THREADPRIVATE(Ccoef_QBS, Dcoef_QBS)
+  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_QBS, Bcoef_QBS
+  !$OMP THREADPRIVATE(Acoef_QBS, Bcoef_QBS)
+  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefqbs
+  !$OMP THREADPRIVATE(Kcoefqbs)
+
+CONTAINS
+!
+!****************************************************************************************
+!
+  SUBROUTINE climb_qbs_down(knon, coefqbs, paprs, pplay, &
+       delp, temp, qbs, dtime, &
+       Ccoef_QBS_out, Dcoef_QBS_out, &
+       Kcoef_qbs_out, gama_qbs_out, &
+       Acoef_QBS_out, Bcoef_QBS_out)
+
+! This routine calculates recursivly the coefficients C and D
+! for the quantity X=[QBS] in equation X(k) = C(k) + D(k)*X(k-1), where k is
+! the index of the vertical layer.
+!
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: knon
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coefqbs
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay 
+    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs 
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: delp  
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: qbs
+    REAL, INTENT(IN)                         :: dtime
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_QBS_out
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_QBS_out
+
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_QBS_out
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_QBS_out
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Kcoef_qbs_out
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_qbs_out
+
+! Local variables
+!****************************************************************************************
+    LOGICAL, SAVE                            :: first=.TRUE.
+    !$OMP THREADPRIVATE(first)
+    REAL, DIMENSION(klon)                    :: psref 
+    REAL                                     :: delz, pkh
+    INTEGER                                  :: k, i, ierr
+! Include
+!****************************************************************************************
+    INCLUDE "YOMCST.h"
+    INCLUDE "compbl.h"    
+
+
+!****************************************************************************************
+! 1)
+! Allocation at first time step only
+!   
+!****************************************************************************************
+
+    IF (first) THEN
+       first=.FALSE.
+       ALLOCATE(Ccoef_QBS(klon,klev), STAT=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_QBS, ierr=', ierr
+       
+       ALLOCATE(Dcoef_QBS(klon,klev), STAT=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_QBS, ierr=', ierr
+       
+       ALLOCATE(Acoef_QBS(klon), Bcoef_QBS(klon), STAT=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_BS and Bcoef_BS, ierr=', ierr
+       
+       ALLOCATE(Kcoefqbs(klon,klev), STAT=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefqbs, ierr=', ierr
+       
+       ALLOCATE(gamaqbs(1:klon,2:klev), STAT=ierr)
+       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaqbs, ierr=', ierr
+       
+    END IF
+
+!****************************************************************************************
+! 2)
+! Definition of the coeficient K 
+!
+!****************************************************************************************
+    Kcoefqbs(:,:) = 0.0
+    DO k = 2, klev
+       DO i = 1, knon
+          Kcoefqbs(i,k) = &
+               coefqbs(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) &
+               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
+       ENDDO
+    ENDDO
+
+!****************************************************************************************
+! 3)
+! Calculation of gama for "Q" and "H"
+!
+!****************************************************************************************
+!   surface pressure is used as reference
+    psref(:) = paprs(:,1) 
+
+!   definition of gama
+    IF (iflag_pbl == 1) THEN
+       gamaqbs(:,:) = 0.0
+ 
+! conversion de gama
+       DO k = 2, klev
+          DO i = 1, knon
+             delz = RD * (temp(i,k-1)+temp(i,k)) / & 
+                    2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
+             pkh  = (psref(i)/paprs(i,k))**RKAPPA
+          
+! convertie gradient verticale de contenu en neige soufflee en difference de neige soufflee entre centre de couches
+             gamaqbs(i,k) = gamaqbs(i,k) * delz    
+          ENDDO
+       ENDDO
+
+    ELSE
+       gamaqbs(:,:) = 0.0
+    ENDIF
+    
+
+!****************************************************************************************    
+! 4)
+! Calculte the coefficients C and D for specific content of blowing snow, qbs
+!
+!****************************************************************************************
+    
+    CALL calc_coef_qbs(knon, Kcoefqbs(:,:), gamaqbs(:,:), delp(:,:), qbs(:,:), &
+         Ccoef_QBS(:,:), Dcoef_QBS(:,:), Acoef_QBS, Bcoef_QBS)
+
+ 
+!****************************************************************************************
+! 5)
+! Return the first layer in output variables
+!
+!****************************************************************************************
+    Acoef_QBS_out = Acoef_QBS
+    Bcoef_QBS_out = Bcoef_QBS
+
+!****************************************************************************************
+! 6)
+! If Pbl is split, return also the other layers in output variables
+!
+!****************************************************************************************
+    IF (mod(iflag_pbl_split,10) .ge.1) THEN
+    DO k= 1, klev
+      DO i= 1, klon
+        Ccoef_QBS_out(i,k) = Ccoef_QBS(i,k)
+        Dcoef_QBS_out(i,k) = Dcoef_QBS(i,k)
+        Kcoef_qbs_out(i,k) = Kcoefqbs(i,k)
+          IF (k.eq.1) THEN
+            gama_qbs_out(i,k)  = 0.
+          ELSE
+            gama_qbs_out(i,k)  = gamaqbs(i,k)
+          ENDIF
+      ENDDO
+    ENDDO
+!!!      
+       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
+!!!
+
+  END SUBROUTINE climb_qbs_down
+!
+!****************************************************************************************
+!
+  SUBROUTINE calc_coef_qbs(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
+!
+! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
+! where X is QQBS, and k the vertical level k=1,klev
+!
+    INCLUDE "YOMCST.h"
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: knon
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
+    REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
+
+! Local variables
+!****************************************************************************************
+    INTEGER                                  :: k, i
+    REAL                                     :: buf
+
+!****************************************************************************************
+! Niveau au sommet, k=klev
+!
+!****************************************************************************************
+    Ccoef(:,:) = 0.0
+    Dcoef(:,:) = 0.0
+
+    DO i = 1, knon
+       buf = delp(i,klev) + Kcoef(i,klev)
+       
+       Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
+       Dcoef(i,klev) = Kcoef(i,klev)/buf
+    END DO
+
+
+!****************************************************************************************
+! Niveau  (klev-1) <= k <= 2
+!
+!****************************************************************************************
+
+    DO k=(klev-1),2,-1
+       DO i = 1, knon
+          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
+          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
+               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
+          Dcoef(i,k) = Kcoef(i,k)/buf
+       END DO
+    END DO
+
+!****************************************************************************************
+! Niveau k=1
+!
+!****************************************************************************************
+
+    DO i = 1, knon
+       buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
+       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
+       Bcoef(i) = -1. * RG / buf
+    END DO
+
+  END SUBROUTINE calc_coef_qbs
+!
+!****************************************************************************************
+!
+  SUBROUTINE climb_qbs_up(knon, dtime, qbs_old, &
+       flx_qbs1, paprs, pplay, &
+       Acoef_QBS_in, Bcoef_QBS_in, &
+       Ccoef_QBS_in, Dcoef_QBS_in, &
+       Kcoef_qbs_in, gama_qbs_in, &
+       flux_qbs, d_qbs)
+! 
+! This routine calculates the flux and tendency of the specific content of blowing snow qbs 
+! The quantity qbs is calculated according to 
+! X(k) = C(k) + D(k)*X(k-1) for X=[qbs], where the coefficients 
+! C and D are known from before and k is index of the vertical layer.
+!   
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                      :: knon
+    REAL, INTENT(IN)                         :: dtime
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: qbs_old
+    REAL, DIMENSION(klon), INTENT(IN)        :: flx_qbs1
+    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
+
+!!! nrlmd le 02/05/2011
+    REAL, DIMENSION(klon), INTENT(IN)        :: Acoef_QBS_in, Bcoef_QBS_in
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Ccoef_QBS_in, Dcoef_QBS_in
+    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef_qbs_in, gama_qbs_in
+!!!
+
+! Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_qbs, d_qbs
+
+! Local variables
+!****************************************************************************************
+    LOGICAL, SAVE                            :: last=.FALSE.
+    REAL, DIMENSION(klon,klev)               :: qbs_new
+    REAL, DIMENSION(klon)                    :: psref         
+    INTEGER                                  :: k, i, ierr
+ 
+! Include
+!****************************************************************************************
+    INCLUDE "YOMCST.h"
+    INCLUDE "compbl.h"    
+
+!****************************************************************************************
+! 1) 
+! Definition of some variables
+    REAL, DIMENSION(klon,klev)               :: zairm
+!
+!****************************************************************************************
+    flux_qbs(:,:) = 0.0
+    d_qbs(:,:)    = 0.0
+
+    psref(1:knon) = paprs(1:knon,1)  
+
+       IF (mod(iflag_pbl_split,10) .ge.1) THEN
+    DO i = 1, knon
+      Acoef_QBS(i)=Acoef_QBS_in(i)
+      Bcoef_QBS(i)=Bcoef_QBS_in(i)
+    ENDDO
+    DO k = 1, klev
+      DO i = 1, knon
+        Ccoef_QBS(i,k)=Ccoef_QBS_in(i,k)
+        Dcoef_QBS(i,k)=Dcoef_QBS_in(i,k)
+        Kcoefqbs(i,k)=Kcoef_qbs_in(i,k)
+          IF (k.gt.1) THEN
+            gamaqbs(i,k)=gama_qbs_in(i,k)
+          ENDIF
+      ENDDO
+    ENDDO
+!!!      
+       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
+!!!
+
+!****************************************************************************************
+! 2)
+! Calculation of QBS
+!
+!****************************************************************************************
+
+!- First layer
+    qbs_new(1:knon,1) = Acoef_QBS(1:knon) + Bcoef_QBS(1:knon)*flx_qbs1(1:knon)*dtime
+!- All the other layers 
+    DO k = 2, klev
+       DO i = 1, knon
+          qbs_new(i,k) = Ccoef_QBS(i,k) + Dcoef_QBS(i,k)*qbs_new(i,k-1)
+       END DO
+    END DO
+!****************************************************************************************
+! 3)
+! Calculation of the flux for QBS
+!
+!****************************************************************************************
+
+!- The flux at first layer, k=1
+    flux_qbs(1:knon,1)=flx_qbs1(1:knon)
+
+!- The flux at all layers above surface
+    DO k = 2, klev
+       DO i = 1, knon
+          flux_qbs(i,k) = (Kcoefqbs(i,k)/RG/dtime) * &
+               (qbs_new(i,k)-qbs_new(i,k-1)+gamaqbs(i,k))
+       END DO
+    END DO
+
+!****************************************************************************************
+! 4)
+! Calculation of tendency for QBS
+!
+!****************************************************************************************
+    DO k = 1, klev
+       DO i = 1, knon
+          d_qbs(i,k) = qbs_new(i,k) - qbs_old(i,k)
+          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
+        END DO
+    END DO
+
+!****************************************************************************************
+! Some deallocations
+!
+!****************************************************************************************
+    IF (last) THEN
+       DEALLOCATE(Ccoef_QBS, Dcoef_QBS,stat=ierr)    
+       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_QBS, Dcoef_QBS, ierr=', ierr
+       DEALLOCATE(Acoef_QBS, Bcoef_QBS,stat=ierr)    
+       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_QBS, Bcoef_QBS, ierr=', ierr
+       DEALLOCATE(gamaqbs,stat=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaqbs, ierr=', ierr
+       DEALLOCATE(Kcoefqbs,stat=ierr)
+       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefqbs, ierr=', ierr
+    END IF
+  END SUBROUTINE climb_qbs_up
+!
+!****************************************************************************************
+!
+END MODULE climb_qbs_mod
+
+ 
+
+
+
+
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/lscp_ini_mod.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/lscp_ini_mod.F90	(revision 4584)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/lscp_ini_mod.F90	(revision 4585)
@@ -5,14 +5,19 @@
   ! PARAMETERS for lscp:
   !--------------------
-  
+ 
+  REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
+  !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG)
+
   REAL, SAVE :: seuil_neb=0.001                 ! cloud fraction threshold: a cloud really exists when exceeded
   !$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", 
                                                 ! 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
+                                                ! >=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)
 
@@ -38,11 +43,89 @@
   !$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)
+
+  REAL, SAVE :: t_glace_max=273.15               ! upper-bound temperature parameter for cloud phase determination
+  !$OMP THREADPRIVATE(t_glace_max)
+
+  REAL, SAVE :: exposant_glace=1.0               ! parameter for cloud phase determination
+  !$OMP THREADPRIVATE(exposant_glace)
+
+  INTEGER, SAVE :: iflag_vice=0                  ! which expression for ice crystall fall velocity
+  !$OMP THREADPRIVATE(iflag_vice)
+
+  INTEGER, SAVE :: iflag_t_glace=0               ! which expression for cloud phase partitioning
+  !$OMP THREADPRIVATE(iflag_t_glace)
+
+  INTEGER, SAVE :: iflag_cloudth_vert=0          ! option for determining cloud fraction and content in convective boundary layers
+  !$OMP THREADPRIVATE(iflag_cloudth_vert)
+
+  INTEGER, SAVE :: iflag_gammasat=0              ! which threshold for homogeneous nucleation below -40oC
+  !$OMP THREADPRIVATE(iflag_gammasat)
+
+  INTEGER, SAVE :: iflag_rain_incloud_vol=0      ! use of volume cloud fraction for rain autoconversion 
+  !$OMP THREADPRIVATE(iflag_rain_incloud_vol)
+
+  INTEGER, SAVE :: iflag_bergeron=0              ! bergeron effect for liquid precipitation treatment  
+  !$OMP THREADPRIVATE(iflag_bergeron)
+
+  INTEGER, SAVE :: iflag_fisrtilp_qsat=0         ! qsat adjustment (iterative) during autoconversion
+  !$OMP THREADPRIVATE(iflag_fisrtilp_qsat)
+
+  INTEGER, SAVE :: iflag_pdf=0                   ! type of subgrid scale qtot pdf
+  !$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)
+
+  REAL, SAVE :: cld_lc_lsc=2.6e-4                ! liquid autoconversion coefficient, stratiform rain
+  !$OMP THREADPRIVATE(cld_lc_lsc)
+
+  REAL, SAVE :: cld_lc_con=2.6e-4                ! liquid autoconversion coefficient, convective rain
+  !$OMP THREADPRIVATE(cld_lc_con)
+
+  REAL, SAVE :: cld_tau_lsc=3600.                ! liquid autoconversion timescale, stratiform rain
+  !$OMP THREADPRIVATE(cld_tau_lsc)
+
+  REAL, SAVE :: cld_tau_con=3600.                ! liquid autoconversion timescale, convective rain
+  !$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)
+
+  REAL, SAVE :: ffallv_con=1.                    ! tuning coefficient crystal fall velocity, convective
+  !$OMP THREADPRIVATE(ffallv_con)
+
+  REAL, SAVE :: coef_eva=2e-5                    ! tuning coefficient liquid precip evaporation
+  !$OMP THREADPRIVATE(coef_eva)
+
+  REAL, SAVE :: coef_eva_i                       ! tuning coefficient ice precip sublimation
+  !$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)
+
+  REAL, SAVE :: dist_liq=300.                    ! typical deph of cloud-top liquid layer in mpcs
+  !$OMP THREADPRIVATE(dist_liq)
+
+  REAL, SAVE    :: tresh_cl=0.0                  ! cloud fraction threshold for cloud top search
+  !$OMP THREADPRIVATE(tresh_cl)
 
 CONTAINS
 
-SUBROUTINE lscp_ini(dtime,ok_ice_sursat)
+SUBROUTINE lscp_ini(dtime,ok_ice_sursat, RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, &
+                    RVTMP2_in, RTT_in,RD_in,RG_in)
 
 
@@ -54,5 +137,23 @@
    LOGICAL, INTENT(IN)   :: ok_ice_sursat  
 
-    CALL getin_p('ninter',ninter)
+   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
+
+
+    RG=RG_in
+    RD=RD_in
+    RCPD=RCPD_in
+    RLVTT=RLVTT_in
+    RLSTT=RLSTT_in
+    RLMLT=RLMLT_in
+    RTT=RTT_in
+    RG=RG_in
+    RVTMP2=RVTMP2_in
+
+
+
+    CALL getin_p('niter_lscp',niter_lscp)
     CALL getin_p('iflag_evap_prec',iflag_evap_prec)
     CALL getin_p('seuil_neb',seuil_neb)
@@ -60,6 +161,35 @@
     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)    
-    WRITE(lunout,*) 'lscp, ninter:', ninter
+    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_cloudth_vert',iflag_cloudth_vert)
+    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('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_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)
+    CALL getin_p('coef_eva',coef_eva)
+    coef_eva_i=coef_eva
+    CALL getin_p('coef_eva_i',coef_eva_i)
+    CALL getin_p('iflag_autoconversion',iflag_autoconversion)
+    CALL getin_p('dist_liq',dist_liq)
+    CALL getin_p('tresh_cl',tresh_cl)
+
+
+
+
+    WRITE(lunout,*) 'lscp, niter_lscp:', niter_lscp
     WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec
     WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb
@@ -67,10 +197,48 @@
     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
+    WRITE(lunout,*) 'lscp, exposant_glace:', exposant_glace
+    WRITE(lunout,*) 'lscp, iflag_vice:', iflag_vice
+    WRITE(lunout,*) 'lscp, iflag_t_glace:', iflag_t_glace
+    WRITE(lunout,*) 'lscp, iflag_cloudth_vert:', iflag_cloudth_vert
+    WRITE(lunout,*) 'lscp, iflag_gammasat:', iflag_gammasat
+    WRITE(lunout,*) 'lscp, iflag_rain_incloud_vol:', iflag_rain_incloud_vol
+    WRITE(lunout,*) 'lscp, iflag_bergeron:', iflag_bergeron
+    WRITE(lunout,*) 'lscp, iflag_fisrtilp_qsat:', iflag_fisrtilp_qsat
+    WRITE(lunout,*) 'lscp, iflag_pdf', iflag_pdf
+    WRITE(lunout,*) 'lscp, reevap_ice', reevap_ice
+    WRITE(lunout,*) 'lscp, cld_lc_lsc', cld_lc_lsc
+    WRITE(lunout,*) 'lscp, cld_lc_con', cld_lc_con
+    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
+    WRITE(lunout,*) 'lscp, dist_liq', dist_liq
+    WRITE(lunout,*) 'lscp, tresh_cl', tresh_cl
+
+
+
+
 
     ! 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/branches/Portage_acc/libf/phylmdiso/thermcell_down.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/thermcell_down.F90	(revision 4584)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/thermcell_down.F90	(revision 4585)
@@ -58,5 +58,6 @@
 
    if ( iflag_thermals_down < 10 ) then
-        stop 'thermcell_down_dq = 0 or >= 10'
+      call abort_physic("thermcell_updown_dq", &
+           'thermcell_down_dq = 0 or >= 10', 1)
    else
         iflag_impl=iflag_thermals_down-10
@@ -128,6 +129,6 @@
          !!!! tentative de prise en compte d'un flux compensatoire montant  !!!!
          if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
-            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
-            stop          
+            call abort_physic("thermcell_updown_dq", 'flux compensatoire '&
+                 // 'montant, cas non traite par thermcell_updown_dq', 1)
             !fthe(ig,ilay)=(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1)
          else
@@ -161,6 +162,6 @@
        do ig=1,ngrid
          if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
-            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
-            stop          
+            call abort_physic("thermcell_updown_dq", 'flux compensatoire ' &
+                 // 'montant, cas non traite par thermcell_updown_dq', 1)
          else
             fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
@@ -191,5 +192,5 @@
          if((fup(ig,ilay)-fdn(ig,ilay)) .lt. 0) then
             write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq dans le cas d une resolution implicite, ilay : ', ilay
-            stop  
+            call abort_physic("thermcell_updown_dq", "", 1)
          else
            mstar_inv=ptimestep/masse(ig,ilay)
@@ -205,7 +206,6 @@
 
    else
-      write(*,*) 'valeur de iflag_impl non prevue'
-      stop
-
+      call abort_physic("thermcell_updown_dq", &
+           'valeur de iflag_impl non prevue', 1)
    endif
 
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/wake_ini_mod.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/wake_ini_mod.F90	(revision 4584)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/wake_ini_mod.F90	(revision 4585)
@@ -25,8 +25,8 @@
   INTEGER, SAVE, PROTECTED                                         :: prt_level
   REAL, SAVE, PROTECTED, DIMENSION(2)                              :: wdens_ref
-  REAL, SAVE, PROTECTED                                            :: stark, coefgw, alpk, pupperbyphs
+  REAL, SAVE, PROTECTED                                            :: stark, coefgw, alpk, wk_pupper
 !>jyg
   REAL, SAVE, PROTECTED                                            :: crep_upper, crep_sol  
-  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, pupperbyphs, crep_upper, crep_sol)
+  !$OMP THREADPRIVATE(stark, wdens_ref, coefgw, alpk, wk_pupper, crep_upper, crep_sol)
 
   REAL, SAVE, PROTECTED                                            :: tau_cv
@@ -153,6 +153,6 @@
   CALL getin_p('alpk',alpk)
   
-  pupperbyphs=0.6
-  CALL getin_p('pupperbyphs',pupperbyphs)
+  wk_pupper=0.6
+  CALL getin_p('wk_pupper',wk_pupper)
 
 
@@ -195,5 +195,5 @@
   WRITE(*,*) 'stark=', stark
   WRITE(*,*) 'alpk=', alpk
-  WRITE(*,*) 'pupperbyphs=', pupperbyphs
+  WRITE(*,*) 'wk_pupper=', wk_pupper
 !jyg<
 !!  WRITE(*,*) 'wdens_ref=', wdens_ref
Index: /LMDZ6/branches/Portage_acc/libf/phylmdiso/yamada_ini_mod.F90
===================================================================
--- /LMDZ6/branches/Portage_acc/libf/phylmdiso/yamada_ini_mod.F90	(revision 4585)
+++ /LMDZ6/branches/Portage_acc/libf/phylmdiso/yamada_ini_mod.F90	(revision 4585)
@@ -0,0 +1,91 @@
+MODULE yamada_ini_mod
+
+! AFAIRE : frif defini deux fois.
+! peut on passer une fonction du module au code.
+
+implicit none
+
+save
+
+  LOGICAL :: new_yamada4
+  INTEGER :: yamada4_num
+  !$OMP THREADPRIVATE(new_yamada4,yamada4_num)
+  LOGICAL :: hboville=.TRUE.
+  REAL :: viscom,viscoh,pbl_lmixmin_alpha
+  !$OMP THREADPRIVATE( hboville,viscom,viscoh,pbl_lmixmin_alpha)
+  real :: ric, yun,ydeux,lmixmin,b1,kap,kapb
+  !$OMP THREADPRIVATE(ric, yun,ydeux,lmixmin,b1,kap,kapb)
+
+  real, dimension(:), allocatable, save :: l0
+  !$OMP THREADPRIVATE(l0)
+  integer :: lunout,prt_level
+  !$OMP THREADPRIVATE(lunout,prt_level)
+
+
+CONTAINS
+
+SUBROUTINE yamada_ini(klon,lunout_in,prt_level_in)
+
+   USE ioipsl_getin_p_mod, ONLY : getin_p
+
+! local
+  integer, intent(in) :: klon,lunout_in,prt_level_in
+  real :: ric0
+  real frif,falpha,fsm,ri
+
+  lunout=lunout_in
+  prt_level=prt_level_in
+
+  b1=16.6
+  kap=0.4
+  kapb=0.4
+
+    ALLOCATE (l0(klon))
+
+   new_yamada4=.false.
+   CALL getin_p('new_yamada4',new_yamada4)
+
+
+    IF (new_yamada4) THEN
+! Corrections et reglages issus du travail de these d'Etienne Vignon.
+       ric=0.143 ! qui donne des valeurs proches des seuils proposes
+                 ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
+                 !                    sm=1.1213 (1.12  dans Y83)
+
+       CALL getin_p('yamada4_ric',ric)
+       ric0=0.19489      ! ric=0.195 originalement, mais produisait sm<0
+       ric=min(ric,ric0) ! Au dela de ric0, sm devient n??gatif
+       yun=1.
+       ydeux=2.
+       hboville=.FALSE.
+       viscom=1.46E-5
+       viscoh=2.06E-5
+       !lmixmin=1.0E-3
+       lmixmin=0.
+       yamada4_num=5
+    ELSE
+       ric=0.195
+       yun=2.
+       ydeux=1.
+       hboville=.TRUE.
+       viscom=0.
+       viscoh=0.
+       lmixmin=1.
+       yamada4_num=0
+    ENDIF
+
+
+!    WRITE(lunout,*)'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
+    CALL getin_p('lmixmin',lmixmin)
+    CALL getin_p('yamada4_hboville',hboville)
+    CALL getin_p('yamada4_num',yamada4_num)
+    pbl_lmixmin_alpha=0.0
+    CALL getin_p('pbl_lmixmin_alpha',pbl_lmixmin_alpha)
+
+
+   
+ RETURN
+
+END SUBROUTINE yamada_ini
+
+END MODULE yamada_ini_mod
