Changeset 2721 for trunk/LMDZ.GENERIC/libf/phystd/rain_generic.F90

Timestamp:

Jun 21, 2022, 11:05:45 AM (2 years ago)

Author:

aslmd

Message:

new routine generic_rain integrated in physiq_mod

File:

• trunk/LMDZ.GENERIC/libf/phystd/rain_generic.F90 (modified) (7 diffs)

trunk/LMDZ.GENERIC/libf/phystd/rain_generic.F90

@@ -1 @@
-subroutine rain_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain_generic,dqsrain_generic,dqssnow_generic,reevap_precip,rneb)
+subroutine rain_generic(ngrid,nlayer,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dq_rain_generic_vap,dq_rain_generic_cld,dqsrain_generic,dqssnow_generic,reevap_precip,rneb)
 
 
    use ioipsl_getin_p_mod, only: getin_p
    use generic_cloud_common_h
-   use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater ! T_h2O_ice_clouds,rhowater  for   precip_scheme_generic >=2
-   use radii_mod, only: h2o_cloudrad ! precip_scheme_generic >=2
+   use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater
+   ! T_h2O_ice_clouds,rhowater  are only used for precip_scheme_generic >=2
+   use radii_mod, only: h2o_cloudrad ! only used for precip_scheme_generic >=2
    use tracer_h
    use comcstfi_mod, only: g, r, cpp
+   use generic_tracer_index_mod, only: generic_tracer_index
    implicit none
   
@@ -36 +38 @@
    real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendency on tracers
    real,intent(out) :: d_t(ngrid,nlayer) ! temperature tendency (K/s)
-   real,intent(out) :: dqrain_generic(ngrid,nlayer,nq) ! tendency of generic condensable tracers precipitation (kg/kg.s-1)
-   real,intent(out) :: dqsrain_generic(ngrid)  ! rain flux at the surface (kg.m-2.s-1)
-   real,intent(out) :: dqssnow_generic(ngrid)  ! snow flux at the surface (kg.m-2.s-1)
-   real,intent(out) :: reevap_precip(ngrid)  ! re-evaporation flux of precipitation integrated over the atmospheric column (kg.m-2.s-1)
+   real,intent(out) :: dq_rain_generic_vap(ngrid,nlayer,nq) ! tendency of GCS tracers precipitation (kg/kg.s-1) - vapor
+   real,intent(out) :: dq_rain_generic_cld(ngrid,nlayer,nq) ! tendency of GCS tracers precipitation (kg/kg.s-1) - cloud
+   real,intent(out) :: dqsrain_generic(ngrid,nq)  ! rain flux at the surface (kg.m-2.s-1)
+   real,intent(out) :: dqssnow_generic(ngrid,nq)  ! snow flux at the surface (kg.m-2.s-1)
+   real,intent(out) :: reevap_precip(ngrid,nq)  ! re-evaporation flux of precipitation integrated over the atmospheric column (kg.m-2.s-1)
    real,intent(in) :: rneb(ngrid,nlayer) ! cloud fraction
 
@@ -48 +51 @@
    real d_ql(ngrid,nlayer)       ! liquid GCS / ice increment
 
-   integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice of generic_tracer
-
-   real, save :: RCPD
-   real, save :: metallicity ! metallicity of planet --- is not used here but necessary to call function Psat_generic
+   integer igcm_generic_vap, igcm_generic_ice ! index of the vap and ice fo GCS
+
+   real, save :: RCPD ! equal to cpp
+
+   real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic
+   !$OMP THREADPRIVATE(metallicity)
 
    ! Subroutine options
@@ -72 +77 @@
    !$OMP THREADPRIVATE(precip_scheme_generic,rainthreshold_generic,cloud_sat_generic,precip_timescale_generic,Cboucher_generic,c1)
 
-   integer,parameter :: ninter=5
-
-   logical,save :: evap_prec_generic ! Does the rain evaporate?
+   integer,parameter :: ninter=5 ! only used for precip_scheme_generic >=2
+
+   logical,save :: evap_prec_generic ! Does the rain evaporate ?
    REAL,SAVE :: evap_coeff_generic ! multiplication evaporation constant. 1. gives the model of gregory et al.
    !$OMP THREADPRIVATE(evap_prec_generic evap_coeff_generic)
 
-   ! for simple scheme
+   ! for simple scheme : precip_scheme_generic=1
    real,parameter :: t_crit=218.0 ! need to be modified to be generic for all species
    real lconvert
@@ -109 +114 @@
    !$OMP THREADPRIVATE(firstcall)
 
+   ! to call only one time the ice/vap pair of a tracer
+   logical one_call_ice_vap_generic
+
    ! Online functions
    real fallv, fall2v, zzz ! falling speed of ice crystals
@@ -114 +122 @@
    fall2v (zzz) =10.6 * ((zzz)**0.31)  !for use with radii
 
-   ! Let's loop on tracers 
+   ! Let's loop on tracers
 
    do iq=1,nq
-      if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0) &
-                                 .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"co2") .eq. 0)) then
+
+      one_call_ice_vap_generic = .false.
+
+      call generic_tracer_index(nq,iq,igcm_generic_vap,igcm_generic_ice,one_call_ice_vap_generic)
+
+      if(one_call_ice_vap_generic) then ! to call only one time the ice/vap pair of a tracer
+
          ! Let's get the index of our tracers (we look for igcm _generic_vap and igcm_generic_ice)
          igcm_generic_vap=iq
@@ -135 +148 @@
          end if
          if (igcm_generic_ice .eq. -1) then
-                     write(*,*) "ERROR : You set a vap traceur but you forgot to set the corresponding ice traceur, &
-                     or the pair vap/ice is not written one after another in traceur.def"
+            write(*,*) "ERROR : You set a vap traceur but you forgot to set the corresponding ice traceur, &
+            or the pair vap/ice is not written one after another in traceur.def"
          endif
-      end if ! if((is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0))
-
-      m=constants_mass(iq)
-      delta_vapH=constants_delta_vapH(iq)
-      Tref=constants_Tref(iq)
-      Pref=constants_Pref(iq)
-      epsi=constants_epsi(iq)
-      RLVTT=constants_RLVTT(iq)
-
-      RCPD = cpp
       
-   end do ! do iq=1,nq loop on tracers 
-
-   if (firstcall) then
-
-      metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
-      call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
-
-      i_vap_generic=igcm_generic_vap
-      i_ice_generic=igcm_generic_ice
-      
-      write(*,*) "rain: i_ice_generic=",i_ice_generic
-      write(*,*) "      i_vap_generic=",i_vap_generic
-
-      PRINT*, 'in rain_generic.F, ninter=', ninter
-      PRINT*, 'in rain_generic.F, evap_prec_generic=', evap_prec_generic
-
-      write(*,*) "Precipitation scheme to use?"
-      precip_scheme_generic=1 ! default value
-      call getin_p("precip_scheme_generic",precip_scheme_generic)
-      write(*,*) " precip_scheme_generic = ",precip_scheme_generic
-
-      if (precip_scheme_generic.eq.1) then
-         write(*,*) "rainthreshold_generic in simple scheme?"
-         rainthreshold_generic=0. ! default value
-         call getin_p("rainthreshold_generic",rainthreshold_generic)
-         write(*,*) " rainthreshold_generic = ",rainthreshold_generic
-      
-      !else 
-      !   write(*,*) "precip_scheme_generic = ", precip_scheme_generic
-      !   write(*,*) "this precip_scheme_generic has not been implemented yet,"
-      !   write(*,*) "only precip_scheme_generic = 1 has been implemented"
+         m=constants_mass(iq)
+         delta_vapH=constants_delta_vapH(iq)
+         Tref=constants_Tref(iq)
+         Pref=constants_Pref(iq)
+         epsi=constants_epsi(iq)
+         RLVTT=constants_RLVTT(iq)
+
+         RCPD = cpp
          
-      else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then
-         
-         write(*,*) "cloud GCS saturation level in non simple scheme?"
-         cloud_sat_generic=2.6e-4   ! default value
-         call getin_p("cloud_sat_generic",cloud_sat_generic)
-         write(*,*) " cloud_sat_generic = ",cloud_sat_generic
-      
-         write(*,*) "precipitation timescale in non simple scheme?"
-         precip_timescale_generic=3600.  ! default value
-         call getin_p("precip_timescale_generic",precip_timescale_generic)
-         write(*,*) " precip_timescale_generic = ",precip_timescale_generic
-
-      else if (precip_scheme_generic.eq.4) then
-         
-         write(*,*) "multiplicative constant in Boucher 95 precip scheme"
-         Cboucher_generic=1.   ! default value
-         call getin_p("Cboucher_generic",Cboucher_generic)
-         write(*,*) " Cboucher_generic = ",Cboucher_generic     
-         
-         c1=1.00*1.097/rhowater*Cboucher_generic*Kboucher  
-
-      endif
-
-      write(*,*) "re-evaporate precipitations?"
-      evap_prec_generic=.true. ! default value
-      call getin_p("evap_prec_generic",evap_prec_generic)
-      write(*,*) " evap_prec_generic = ",evap_prec_generic
-
-      if (evap_prec_generic) then
-         write(*,*) "multiplicative constant in reevaporation"
-         evap_coeff_generic=1.   ! default value
-         call getin_p("evap_coeff_generic",evap_coeff_generic)
-         write(*,*) " evap_coeff_generic = ",evap_coeff_generic 
-      end if
-
-      firstcall = .false.
-
-   endif ! of if (firstcall)
-
-   ! GCM -----> subroutine variables
-   do k = 1, nlayer
-      do i = 1, ngrid
-
-         zt(i,k)   = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
-         q(i,k)    = pq(i,k,i_vap_generic)+pdq(i,k,i_vap_generic)*ptimestep
-         ql(i,k)   = pq(i,k,i_ice_generic)+pdq(i,k,i_ice_generic)*ptimestep
-
-         !q(i,k)    = pq(i,k,i_vap_generic)!+pdq(i,k,i_vap_generic)
-         !ql(i,k)   = pq(i,k,i_ice_generic)!+pdq(i,k,i_ice_generic)
-
-         if(q(i,k).lt.0.)then ! if this is not done, we don't conserve GCS
-            q(i,k)=0.
-         endif
-         if(ql(i,k).lt.0.)then
-            ql(i,k)=0.
-         endif
-      enddo
-   enddo
-
-   ! Initialise the outputs
-   d_t(1:ngrid,1:nlayer) = 0.0
-   d_q(1:ngrid,1:nlayer) = 0.0
-   d_ql(1:ngrid,1:nlayer) = 0.0
-   precip_rate(1:ngrid) = 0.0
-   precip_rate_tmp(1:ngrid) = 0.0
-
-   ! calculate saturation mixing ratio
-   do k = 1, nlayer
-      do i = 1, ngrid
-         p_tmp = pplay(i,k)
-         call Psat_generic(zt(i,k),p_tmp,metallicity,psat_tmp,zqs(i,k)) ! metallicit --- is not used here but necessary to call function Psat_generic
-         call Lcpdqsat_generic(zt(i,k),p_tmp,psat_tmp,zqs(i,k),dqsat(i,k),dlnpsat(i,k))
-         call Tsat_generic(p_tmp,Tsat(i,k))
-
-      enddo
-   enddo
-
-   ! get column / layer conversion factor
-   do k = 1, nlayer
-      do i = 1, ngrid
-         dmass(i,k)=(pplev(i,k)-pplev(i,k+1))/g
-      enddo
-   enddo
-
-   ! Vertical loop (from top to bottom)
-   ! We carry the rain with us and calculate that added by warm/cold precipitation
-   ! processes and that subtracted by evaporation at each level.
-   do k = nlayer, 1, -1
-
-      if (evap_prec_generic) then ! note no rneb dependence!
-         do i = 1, ngrid
-            if (precip_rate(i) .GT.0.) then
-
-               if(zt(i,k).gt.Tsat(i,k))then
-                  ! treat the case where all liquid water should boil
-                  zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*dmass(i,k)/RLVTT/ptimestep,precip_rate(i)) ! on évapore
-                  precip_rate(i)=MAX(precip_rate(i)-zqev,0.)
-                  d_q(i,k)=zqev/dmass(i,k)*ptimestep ! quantité évaporée
-                  d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
-               else
-                  zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass(i,k)/(ptimestep*(1.d0+dqsat(i,k)))
-                  !max evaporation to reach saturation implictly accounting for temperature reduction
-                  zqevt= MAX (0.0, evap_coeff_generic*2.0e-5*(1.0-q(i,k)/zqs(i,k))    & !default was 2.e-5
-                        *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*R) ! BC modif here
-                  zqev  = MIN (zqev, zqevt)
-                  zqev  = MAX (zqev, 0.0)
-                  precip_rate_tmp(i)= precip_rate(i) - zqev
-                  precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0)
-
-                  d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmass(i,k)*ptimestep
-                  d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
-                  precip_rate(i)  = precip_rate_tmp(i)
-               end if
-            endif ! of if (precip_rate(i) .GT.0.)
+
+         if (firstcall) then
+
+            metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic
+            call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic
+
+            i_vap_generic=igcm_generic_vap
+            i_ice_generic=igcm_generic_ice
+            
+            write(*,*) "rain: i_ice_generic=",i_ice_generic
+            write(*,*) "      i_vap_generic=",i_vap_generic
+
+            write(*,*) "re-evaporate precipitations?"
+            evap_prec_generic=.true. ! default value
+            call getin_p("evap_prec_generic",evap_prec_generic)
+            write(*,*) " evap_prec_generic = ",evap_prec_generic
+
+            if (evap_prec_generic) then
+               write(*,*) "multiplicative constant in reevaporation"
+               evap_coeff_generic=1.   ! default value
+               call getin_p("evap_coeff_generic",evap_coeff_generic)
+               write(*,*) " evap_coeff_generic = ",evap_coeff_generic   
+            end if
+
+            write(*,*) "Precipitation scheme to use?"
+            precip_scheme_generic=1 ! default value
+            call getin_p("precip_scheme_generic",precip_scheme_generic)
+            write(*,*) " precip_scheme_generic = ",precip_scheme_generic
+
+            if (precip_scheme_generic.eq.1) then
+               write(*,*) "rainthreshold_generic in simple scheme?"
+               rainthreshold_generic=0. ! default value
+               call getin_p("rainthreshold_generic",rainthreshold_generic)
+               write(*,*) " rainthreshold_generic = ",rainthreshold_generic
+            
+            !else 
+            !   write(*,*) "precip_scheme_generic = ", precip_scheme_generic
+            !   write(*,*) "this precip_scheme_generic has not been implemented yet,"
+            !   write(*,*) "only precip_scheme_generic = 1 has been implemented"
+               
+            else if (precip_scheme_generic.eq.2.or.precip_scheme_generic.eq.3) then
+               
+               write(*,*) "cloud GCS saturation level in non simple scheme?"
+               cloud_sat_generic=2.6e-4   ! default value
+               call getin_p("cloud_sat_generic",cloud_sat_generic)
+               write(*,*) " cloud_sat_generic = ",cloud_sat_generic
+            
+               write(*,*) "precipitation timescale in non simple scheme?"
+               precip_timescale_generic=3600.  ! default value
+               call getin_p("precip_timescale_generic",precip_timescale_generic)
+               write(*,*) " precip_timescale_generic = ",precip_timescale_generic
+
+            else if (precip_scheme_generic.eq.4) then
+               
+               write(*,*) "multiplicative constant in Boucher 95 precip scheme"
+               Cboucher_generic=1.   ! default value
+               call getin_p("Cboucher_generic",Cboucher_generic)
+               write(*,*) " Cboucher_generic = ",Cboucher_generic       
+               
+               c1=1.00*1.097/rhowater*Cboucher_generic*Kboucher  
+
+            endif
+
+            PRINT*, 'in rain_generic.F, ninter=', ninter
+
+            firstcall = .false.
+
+         endif ! of if (firstcall)
+
+         ! GCM -----> subroutine variables
+         do k = 1, nlayer
+            do i = 1, ngrid
+
+               zt(i,k)   = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
+               q(i,k)    = pq(i,k,i_vap_generic)+pdq(i,k,i_vap_generic)*ptimestep
+               ql(i,k)   = pq(i,k,i_ice_generic)+pdq(i,k,i_ice_generic)*ptimestep
+
+               !q(i,k)    = pq(i,k,i_vap_generic)!+pdq(i,k,i_vap_generic)
+               !ql(i,k)   = pq(i,k,i_ice_generic)!+pdq(i,k,i_ice_generic)
+
+               if(q(i,k).lt.0.)then ! if this is not done, we don't conserve GCS
+                  q(i,k)=0. ! vap
+               endif
+               if(ql(i,k).lt.0.)then
+                  ql(i,k)=0. ! ice
+               endif
+            enddo
          enddo
-      endif ! of if (evap_prec_generic)
-
-      zoliq(1:ngrid) = 0.0
-
-
-      if(precip_scheme_generic.eq.1)then
-
-         do i = 1, ngrid
-            t_tmp = zt(i,k)
-            if (t_tmp .ge. T_h2O_ice_liq) then
-               lconvert=rainthreshold_generic
-            elseif (t_tmp .gt. t_crit) then
-               lconvert=rainthreshold_generic*(1.- t_crit/t_tmp)
-               lconvert=MAX(0.0,lconvert)              
-            else
-               lconvert=0.
-            endif
-
-
-            if (ql(i,k).gt.1.e-9) then
-               zneb(i)  = MAX(rneb(i,k), seuil_neb)
-               if ((ql(i,k)/zneb(i)).gt.lconvert)then ! precipitate!
-                  d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
-                  precip_rate(i)   = precip_rate(i) - d_ql(i,k)*dmass(i,k)/ptimestep
-               endif
-            endif
+
+         ! Initialise the outputs
+         d_t(1:ngrid,1:nlayer) = 0.0
+         d_q(1:ngrid,1:nlayer) = 0.0
+         d_ql(1:ngrid,1:nlayer) = 0.0
+         precip_rate(1:ngrid) = 0.0
+         precip_rate_tmp(1:ngrid) = 0.0
+
+         ! calculate saturation mixing ratio
+         do k = 1, nlayer
+            do i = 1, ngrid
+               p_tmp = pplay(i,k)
+               call Psat_generic(zt(i,k),p_tmp,metallicity,psat_tmp,zqs(i,k)) ! calculates psat_tmp & zqs(i,k)
+               ! metallicity --- is not used here but necessary to call function Psat_generic
+               call Lcpdqsat_generic(zt(i,k),p_tmp,psat_tmp,zqs(i,k),dqsat(i,k),dlnpsat(i,k)) ! calculates dqsat(i,k) & dlnpsat(i,k)
+               call Tsat_generic(p_tmp,Tsat(i,k)) ! calculates Tsat(i,k)
+
+            enddo
          enddo
 
-      elseif (precip_scheme_generic.ge.2) then
-      
-         do i = 1, ngrid
-            if (rneb(i,k).GT.0.0) then
-               zoliq(i) = ql(i,k)
-               zrho(i)  = pplay(i,k) / ( zt(i,k) * R )
-               zdz(i)   = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
-               zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
-               zfrac(i) = MAX(zfrac(i), 0.0)
-               zfrac(i) = MIN(zfrac(i), 1.0)
-               zneb(i)  = MAX(rneb(i,k), seuil_neb)
-            endif
+         ! get column / layer conversion factor
+         do k = 1, nlayer
+            do i = 1, ngrid
+               dmass(i,k)=(pplev(i,k)-pplev(i,k+1))/g ! mass per m² in each layer
+            enddo
          enddo
 
-         !recalculate liquid GCS particle radii
-         call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice)
-
-         SELECT CASE(precip_scheme_generic)
-
-         !precip scheme from Sundquist 78
-         CASE(2)
-
-            do n = 1, ninter
+         ! Vertical loop (from top to bottom)
+         ! We carry the rain with us and calculate that added by warm/cold precipitation
+         ! processes and that subtracted by evaporation at each level.
+         ! We go from a layer to the next layer below and make the rain fall
+         do k = nlayer, 1, -1
+
+            if (evap_prec_generic) then ! note no rneb dependence!
+               do i = 1, ngrid
+                  if (precip_rate(i) .GT.0.) then ! if rain from upper layers has fallen in the current layer box
+
+                     if(zt(i,k).gt.Tsat(i,k))then ! if temperature of the layer box is greater than Tsat
+                        ! treat the case where all liquid water should boil
+                        zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*dmass(i,k)/RLVTT/ptimestep,precip_rate(i)) ! on évapore
+                        precip_rate(i)=MAX(precip_rate(i)-zqev,0.) ! we withdraw from precip_rate the evaporated part
+                        d_q(i,k)=zqev/dmass(i,k)*ptimestep ! quantité évaporée
+                        d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
+                     else
+                        zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*dmass(i,k)/(ptimestep*(1.d0+dqsat(i,k)))
+                        !max evaporation to reach saturation implictly accounting for temperature reduction
+                        zqevt= MAX (0.0, evap_coeff_generic*2.0e-5*(1.0-q(i,k)/zqs(i,k))    & !default was 2.e-5
+                              *sqrt(precip_rate(i))*dmass(i,k)/pplay(i,k)*zt(i,k)*R) ! BC modif here, is it R or r/(mu/1000) ?
+                        zqev  = MIN (zqev, zqevt)
+                        zqev  = MAX (zqev, 0.0) ! a priori inutile d'après les précédentes lignes
+                        precip_rate_tmp(i)= precip_rate(i) - zqev
+                        precip_rate_tmp(i)= max(precip_rate_tmp(i),0.0)
+
+                        d_q(i,k) = - (precip_rate_tmp(i)-precip_rate(i))/dmass(i,k)*ptimestep
+                        d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
+                        precip_rate(i)  = precip_rate_tmp(i)
+                     end if
+                  endif ! of if (precip_rate(i) .GT.0.)
+               enddo
+            endif ! of if (evap_prec_generic)
+
+            zoliq(1:ngrid) = 0.0
+
+
+            if(precip_scheme_generic.eq.1)then
+
+               do i = 1, ngrid
+                  t_tmp = zt(i,k)
+                  if (t_tmp .ge. T_h2O_ice_liq) then
+                     lconvert=rainthreshold_generic
+                  elseif (t_tmp .gt. t_crit) then
+                     lconvert=rainthreshold_generic*(1.- t_crit/t_tmp)
+                     lconvert=MAX(0.0,lconvert)              
+                  else
+                     lconvert=0.
+                  endif
+
+                  ! lconvert is between 0. and rainthreshold_generic
+
+
+                  if (ql(i,k).gt.1.e-9) then
+                     zneb(i)  = MAX(rneb(i,k), seuil_neb) ! in mesoscale rneb = 0 or 1
+                     if ((ql(i,k)/zneb(i)).gt.lconvert)then ! precipitate!
+                        d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
+                        precip_rate(i)   = precip_rate(i) - d_ql(i,k)*dmass(i,k)/ptimestep
+                     endif
+                  endif
+               enddo
+
+            elseif (precip_scheme_generic.ge.2) then
+            
                do i = 1, ngrid
                   if (rneb(i,k).GT.0.0) then
-                     ! this is the ONLY place where zneb, precip_timescale_generic and cloud_sat_generic are used
-
-                     zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale_generic)) * zoliq(i)      &
-                              * (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat_generic)**2)) * zfrac(i)
-                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
-                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
-                              *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
-                     ztot(i)  = zchau(i) + zfroi(i)
-
-                     if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
-                        ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
-                        zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
-                     endif
+                     zoliq(i) = ql(i,k)
+                     zrho(i)  = pplay(i,k) / ( zt(i,k) * R )
+                     zdz(i)   = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
+                     zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
+                     zfrac(i) = MAX(zfrac(i), 0.0)
+                     zfrac(i) = MIN(zfrac(i), 1.0)
+                     zneb(i)  = MAX(rneb(i,k), seuil_neb)
+                  endif
                enddo
-            enddo            
-
-         !precip scheme modified from Sundquist 78 (in q**3)
-         CASE(3)         
-            do n = 1, ninter
+
+               !recalculate liquid GCS particle radii
+               call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice)
+
+               SELECT CASE(precip_scheme_generic)
+
+               !precip scheme from Sundquist 78
+               CASE(2)
+
+                  do n = 1, ninter
+                     do i = 1, ngrid
+                        if (rneb(i,k).GT.0.0) then
+                           ! this is the ONLY place where zneb, precip_timescale_generic and cloud_sat_generic are used
+
+                           zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale_generic)) * zoliq(i)      &
+                                    * (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat_generic)**2)) * zfrac(i)
+                           zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+                           zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
+                                    *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
+                           ztot(i)  = zchau(i) + zfroi(i)
+
+                           if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
+                              ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
+                              zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
+                           endif
+                     enddo
+                  enddo            
+
+               !precip scheme modified from Sundquist 78 (in q**3)
+               CASE(3)   
+                  do n = 1, ninter
+                     do i = 1, ngrid
+                        if (rneb(i,k).GT.0.0) then
+                        ! this is the ONLY place where zneb, precip_timescale_generic and cloud_sat_generic are used
+
+                           zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale_generic*cloud_sat_generic**2)) * (zoliq(i)/zneb(i))**3  
+                           zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+                           zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
+                                       *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
+                           ztot(i)  = zchau(i) + zfroi(i)
+
+                           if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
+                              ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
+                              zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
+                           endif
+                     enddo
+                  enddo            
+
+               !precip scheme modified from Boucher 95
+               CASE(4)
+                  do n = 1, ninter
+                     do i = 1, ngrid
+                        if (rneb(i,k).GT.0.0) then
+                        ! this is the ONLY place where zneb and c1 are used
+                           zchau(i) = ptimestep/FLOAT(ninter) *c1* zrho(i) &
+                                       *(zoliq(i)/zneb(i))**2*reffh2oliq(i,k)*zneb(i)* zfrac(i) 
+                           zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
+                           zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
+                                       *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
+                           ztot(i)  = zchau(i) + zfroi(i)
+
+                           if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
+                           ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
+                           zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
+                        endif
+                     enddo
+                  enddo            
+
+               END SELECT ! precip_scheme_generic 
+
+               ! Change in cloud density and surface GCS values
                do i = 1, ngrid
                   if (rneb(i,k).GT.0.0) then
-                  ! this is the ONLY place where zneb, precip_timescale_generic and cloud_sat_generic are used
-
-                     zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale_generic*cloud_sat_generic**2)) * (zoliq(i)/zneb(i))**3  
-                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
-                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
-                                 *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
-                     ztot(i)  = zchau(i) + zfroi(i)
-
-                     if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
-                        ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
-                        zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
-                     endif
-               enddo
-            enddo            
-
-         !precip scheme modified from Boucher 95
-         CASE(4)
-            do n = 1, ninter
-               do i = 1, ngrid
-                  if (rneb(i,k).GT.0.0) then
-                  ! this is the ONLY place where zneb and c1 are used
-                     zchau(i) = ptimestep/FLOAT(ninter) *c1* zrho(i) &
-                                 *(zoliq(i)/zneb(i))**2*reffh2oliq(i,k)*zneb(i)* zfrac(i) 
-                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
-                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
-                                 *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
-                     ztot(i)  = zchau(i) + zfroi(i)
-
-                     if (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
-                     ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
-                     zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
+                     d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
+                     precip_rate(i)   = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmass(i,k)/ptimestep
                   endif
                enddo
-            enddo            
-
-         END SELECT ! precip_scheme_generic 
-
-         ! Change in cloud density and surface GCS values
+
+            endif ! if precip_scheme_generic=1
+
+         enddo ! of do k = nlayer, 1, -1
+
+         ! Rain or snow on the ground
          do i = 1, ngrid
-            if (rneb(i,k).GT.0.0) then
-               d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
-               precip_rate(i)   = precip_rate(i)+ MAX(ql(i,k)-zoliq(i),0.0)*dmass(i,k)/ptimestep
+            if(precip_rate(i).lt.0.0)then
+               print*,'Droplets of negative rain are falling...'
+               call abort
             endif
+            if (t(i,1) .LT. T_h2O_ice_liq) then
+               dqssnow_generic(i,i_ice_generic) = precip_rate(i)
+               dqsrain_generic(i,i_ice_generic) = 0.0
+            else
+               dqssnow_generic(i,i_ice_generic) = 0.0
+               dqsrain_generic(i,i_ice_generic) = precip_rate(i) ! liquid water = ice for now
+            endif
+
+            ! For now we force :
+            dqsrain_generic(i,i_ice_generic) = precip_rate(i)
+            dqssnow_generic(i,i_ice_generic) = 0.0
          enddo
 
-      endif ! if precip_scheme_generic=1
-
-   enddo ! of do k = nlayer, 1, -1
-
-   ! Rain or snow on the ground
-   do i = 1, ngrid
-      if(precip_rate(i).lt.0.0)then
-         print*,'Droplets of negative rain are falling...'
-         call abort
-      endif
-      if (t(i,1) .LT. T_h2O_ice_liq) then
-         dqssnow_generic(i) = precip_rate(i)
-         dqsrain_generic(i) = 0.0
-      else
-         dqssnow_generic(i) = 0.0
-         dqsrain_generic(i) = precip_rate(i) ! liquid water = ice for now
-      endif
-   enddo
-
-   ! now subroutine -----> GCM variables
-   if (evap_prec_generic) then
-      dqrain_generic(1:ngrid,1:nlayer,i_vap_generic)=d_q(1:ngrid,1:nlayer)/ptimestep
-      d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep
-      do i=1,ngrid
-         reevap_precip(i)=0.
-         do k=1,nlayer
-            reevap_precip(i)=reevap_precip(i)+dqrain_generic(i,k,i_vap_generic)*dmass(i,k)
-         enddo
-      enddo
-   else
-      dqrain_generic(1:ngrid,1:nlayer,i_vap_generic)=0.0
-      d_t(1:ngrid,1:nlayer)=0.0
-   endif
-   dqrain_generic(1:ngrid,1:nlayer,i_ice_generic) = d_ql(1:ngrid,1:nlayer)/ptimestep
+         ! now subroutine -----> GCM variables
+         if (evap_prec_generic) then
+            dq_rain_generic_vap(1:ngrid,1:nlayer,i_vap_generic)=d_q(1:ngrid,1:nlayer)/ptimestep
+            d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep
+            do i=1,ngrid
+               reevap_precip(i,i_vap_generic)=0.
+               do k=1,nlayer
+                  reevap_precip(i,i_vap_generic)=reevap_precip(i,i_vap_generic)+dq_rain_generic_vap(i,k,i_vap_generic)*dmass(i,k)
+               enddo
+            enddo
+         else
+            dq_rain_generic_vap(1:ngrid,1:nlayer,i_vap_generic)=0.0
+            d_t(1:ngrid,1:nlayer)=0.0
+         endif
+         dq_rain_generic_cld(1:ngrid,1:nlayer,i_ice_generic) = d_ql(1:ngrid,1:nlayer)/ptimestep
+
+      end if ! if(one_call_ice_vap_generic)
+
+   end do ! do iq=1,nq loop on tracers 
 
 end subroutine rain_generic