Skip to content
Snippets Groups Projects
rain_ice.f90 145 KiB
Newer Older

 IF (OWARM) THEN ! rain_ice_warm
   IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'AUTO_BU_RRC')
   IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'AUTO_BU_RRR')
   IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'ACCR_BU_RRC')
   IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'ACCR_BU_RRR')
   IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'REVA_BU_RTH')
   IF (LBUDGET_RV) CALL BUDGET (PRVS(:,:,:)*PRHODJ(:,:,:),6,'REVA_BU_RRV')
   IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'REVA_BU_RRR')
 ENDIF

 !rain_ice_fast_rs
 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'RIM_BU_RTH')
 IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'RIM_BU_RRC')
 IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:),10,'RIM_BU_RRS')
 IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'RIM_BU_RRG')
 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'ACC_BU_RTH')
 IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'ACC_BU_RRR')
 IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:),10,'ACC_BU_RRS')
 IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'ACC_BU_RRG')
 IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:),10,'CMEL_BU_RRS')
 IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'CMEL_BU_RRG')

 !rain_ice_fast_rg
 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'CFRZ_BU_RTH')
 IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'CFRZ_BU_RRR')
 IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:),9,'CFRZ_BU_RRI')
 IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'CFRZ_BU_RRG')
 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'WETG_BU_RTH')
 IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'WETG_BU_RRC')
 IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'WETG_BU_RRR')
 IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:),9,'WETG_BU_RRI')
 IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:),10,'WETG_BU_RRS')
 IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'WETG_BU_RRG')
 IF (LBUDGET_RH) CALL BUDGET (PRHS(:,:,:)*PRHODJ(:,:,:),12,'WETG_BU_RRH')
 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'DRYG_BU_RTH')
 IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'DRYG_BU_RRC')
 IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'DRYG_BU_RRR')
 IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:),9,'DRYG_BU_RRI')
 IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:),10,'DRYG_BU_RRS')
 IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'DRYG_BU_RRG')
 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'GMLT_BU_RTH')
 IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'GMLT_BU_RRR')
 IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'GMLT_BU_RRG')

 IF(KRR==7) THEN ! rain_ice_fast_rh
   IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'WETH_BU_RTH')
   IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'WETH_BU_RRC')
   IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'WETH_BU_RRR')
   IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:),9,'WETH_BU_RRI')
   IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:),10,'WETH_BU_RRS')
   IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'WETH_BU_RRG')
   IF (LBUDGET_RH) CALL BUDGET (PRHS(:,:,:)*PRHODJ(:,:,:),12,'WETH_BU_RRH')
   IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'HMLT_BU_RTH')
   IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8,'HMLT_BU_RRR')
   IF (LBUDGET_RH) CALL BUDGET (PRHS(:,:,:)*PRHODJ(:,:,:),12,'HMLT_BU_RRH')
 ENDIF

 !rain_ice_fast_ri
 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'IMLT_BU_RTH')
 IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'IMLT_BU_RRC')
 IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:),9,'IMLT_BU_RRI')
 IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:)*PRHODJ(:,:,:),4,'BERFI_BU_RTH')
 IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7,'BERFI_BU_RRC')
 IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:),9,'BERFI_BU_RRI')
!
END IF
!
!-------------------------------------------------------------------------------
!
!*       8.     COMPUTE THE SEDIMENTATION (RS) SOURCE
!               -------------------------------------
!
!*       8.1    time splitting loop initialization
!
ZTSPLITR= PTSTEP / FLOAT(KSPLITR)
!
!
IF (HSEDIM == 'STAT') THEN
  CALL RAIN_ICE_SEDIMENTATION_STAT
ELSEIF (HSEDIM == 'SPLI') THEN
  CALL RAIN_ICE_SEDIMENTATION_SPLIT
ELSE
  WRITE(*,*) ' STOP'
  WRITE(*,*) ' NO SEDIMENTATION SCHEME FOR HSEDIM=',HSEDIM
  CALL PRINT_MSG(NVERB_FATAL,'GEN','RAIN_ICE','')  
!sedimentation of rain fraction
CALL RAINFR_VERT(PRAINFR, PRRS(:,:,:)*PTSTEP)

!
!
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
!
CONTAINS
!
!
!-------------------------------------------------------------------------------
!
!
  SUBROUTINE RAIN_ICE_SEDIMENTATION_SPLIT
!
!*      0. DECLARATIONS
!          ------------
!
IMPLICIT NONE
!
!*       0.2  declaration of local variables
!
!
INTEGER , DIMENSION(SIZE(GSEDIMC)) :: IC1,IC2,IC3 ! Used to replace the COUNT
INTEGER , DIMENSION(SIZE(GSEDIMR)) :: IR1,IR2,IR3 ! Used to replace the COUNT
INTEGER , DIMENSION(SIZE(GSEDIMS)) :: IS1,IS2,IS3 ! Used to replace the COUNT
INTEGER , DIMENSION(SIZE(GSEDIMI)) :: II1,II2,II3 ! Used to replace the COUNT
INTEGER , DIMENSION(SIZE(GSEDIMG)) :: IG1,IG2,IG3 ! Used to replace the COUNT
INTEGER , DIMENSION(SIZE(GSEDIMH)) :: IH1,IH2,IH3 ! Used to replace the COUNT
INTEGER   :: ILENALLOCC,ILENALLOCR,ILENALLOCI,ILENALLOCS,ILENALLOCG,ILENALLOCH
INTEGER   :: ILISTLENC,ILISTLENR,ILISTLENI,ILISTLENS,ILISTLENG,ILISTLENH
INTEGER, ALLOCATABLE :: ILISTR(:),ILISTC(:),ILISTI(:),ILISTS(:),ILISTG(:),ILISTH(:)
! Optimization for NEC
!INTEGER, SAVE :: IOLDALLOCC = SIZE(PEXNREF,1)*SIZE(PEXNREF,2)*SIZE(PEXNREF,3)/10
!INTEGER, SAVE :: IOLDALLOCR = SIZE(PEXNREF,1)*SIZE(PEXNREF,2)*SIZE(PEXNREF,3)/10
!INTEGER, SAVE :: IOLDALLOCI = SIZE(PEXNREF,1)*SIZE(PEXNREF,2)*SIZE(PEXNREF,3)/10
!INTEGER, SAVE :: IOLDALLOCS = SIZE(PEXNREF,1)*SIZE(PEXNREF,2)*SIZE(PEXNREF,3)/10
!INTEGER, SAVE :: IOLDALLOCG = SIZE(PEXNREF,1)*SIZE(PEXNREF,2)*SIZE(PEXNREF,3)/10
!INTEGER, SAVE :: IOLDALLOCH = SIZE(PEXNREF,1)*SIZE(PEXNREF,2)*SIZE(PEXNREF,3)/10
INTEGER, SAVE :: IOLDALLOCC = 6000
INTEGER, SAVE :: IOLDALLOCR = 6000
INTEGER, SAVE :: IOLDALLOCI = 6000
INTEGER, SAVE :: IOLDALLOCS = 6000
INTEGER, SAVE :: IOLDALLOCG = 6000
INTEGER, SAVE :: IOLDALLOCH = 6000
!
REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3)) :: ZCONC3D !  droplet condensation
!-------------------------------------------------------------------------------
!
!
!        O. Initialization of for sedimentation
!
IF (OSEDIC) PINPRC (:,:) = 0.
IF (LDEPOSC) PINDEP (:,:) = 0.
PINPRR (:,:) = 0.
PINPRR3D (:,:,:) = 0.
PINPRS (:,:) = 0.
PINPRG (:,:) = 0.
IF ( KRR == 7 ) PINPRH (:,:) = 0.
IF (PRESENT(PFPR)) PFPR(:,:,:,:) = 0.
!
!*       1. Parameters for cloud sedimentation
!
   IF (OSEDIC) THEN
    ZRAY(:,:,:)   = 0.
    ZLBC(:,:,:)   = XLBC(1)
    ZFSEDC(:,:,:) = XFSEDC(1)
    ZCONC3D(:,:,:)= XCONC_LAND
    ZCONC_TMP(:,:)= XCONC_LAND
    IF (PRESENT(PSEA)) THEN 
      ZCONC_TMP(:,:)=PSEA(:,:)*XCONC_SEA+(1.-PSEA(:,:))*XCONC_LAND
      DO JK=IKTB,IKTE
        ZLBC(:,:,JK)   = PSEA(:,:)*XLBC(2)+(1.-PSEA(:,:))*XLBC(1)
        ZFSEDC(:,:,JK) = (PSEA(:,:)*XFSEDC(2)+(1.-PSEA(:,:))*XFSEDC(1))
        ZFSEDC(:,:,JK) = MAX(MIN(XFSEDC(1),XFSEDC(2)),ZFSEDC(:,:,JK))
        ZCONC3D(:,:,JK)= (1.-PTOWN(:,:))*ZCONC_TMP(:,:)+PTOWN(:,:)*XCONC_URBAN
        ZRAY(:,:,JK)   = 0.5*((1.-PSEA(:,:))*GAMMA(XNUC+1.0/XALPHAC)/(GAMMA(XNUC)) + &
                PSEA(:,:)*GAMMA(XNUC2+1.0/XALPHAC2)/(GAMMA(XNUC2)))
      END DO
    ELSE
        ZCONC3D(:,:,:) = XCONC_LAND
        ZRAY(:,:,:)  = 0.5*(GAMMA(XNUC+1.0/XALPHAC)/(GAMMA(XNUC)))
    END IF
    ZRAY(:,:,:)      = MAX(1.,ZRAY(:,:,:))
    ZLBC(:,:,:)      = MAX(MIN(XLBC(1),XLBC(2)),ZLBC(:,:,:))
   ENDIF
!
!*       2.    compute the fluxes
!
!  optimization by looking for locations where
!  the precipitating fields are larger than a minimal value only !!!
!  For optimization we consider each variable separately

ZRTMIN(:)    = XRTMIN(:) * ZINVTSTEP
IF (OSEDIC) GSEDIMC(:,:,:) = .FALSE.
GSEDIMR(:,:,:) = .FALSE.
GSEDIMI(:,:,:) = .FALSE.
GSEDIMS(:,:,:) = .FALSE.
GSEDIMG(:,:,:) = .FALSE.
IF ( KRR == 7 ) GSEDIMH(:,:,:) = .FALSE.
!
ILENALLOCR = 0
IF (OSEDIC) ILENALLOCC = 0
ILENALLOCI = 0
ILENALLOCS = 0
ILENALLOCG = 0
IF ( KRR == 7 ) ILENALLOCH = 0
!
! ZPiS = Specie i source creating during the current time step
! PRiS = Source of the previous time step
!
IF (OSEDIC) THEN
  ZPRCS(:,:,:) = 0.0
  ZPRCS(:,:,:) = PRCS(:,:,:)-PRCT(:,:,:)* ZINVTSTEP
  PRCS(:,:,:)  = PRCT(:,:,:)* ZINVTSTEP
END IF
ZPRRS(:,:,:) = 0.0
ZPRSS(:,:,:) = 0.0
ZPRGS(:,:,:) = 0.0
IF ( KRR == 7 ) ZPRHS(:,:,:) = 0.0
!
ZPRRS(:,:,:) = PRRS(:,:,:)-PRRT(:,:,:)* ZINVTSTEP
ZPRSS(:,:,:) = PRSS(:,:,:)-PRST(:,:,:)* ZINVTSTEP
ZPRGS(:,:,:) = PRGS(:,:,:)-PRGT(:,:,:)* ZINVTSTEP
IF ( KRR == 7 ) ZPRHS(:,:,:) = PRHS(:,:,:)-PRHT(:,:,:)* ZINVTSTEP
PRRS(:,:,:)  = PRRT(:,:,:)* ZINVTSTEP
PRSS(:,:,:)  = PRST(:,:,:)* ZINVTSTEP
PRGS(:,:,:)  = PRGT(:,:,:)* ZINVTSTEP
IF ( KRR == 7 ) PRHS(:,:,:)  = PRHT(:,:,:)* ZINVTSTEP
!
! PRiS = Source of the previous time step + source created during the subtime
! step
!
DO JN = 1 , KSPLITR
  IF( JN==1 ) THEN
   IF (OSEDIC) PRCS(:,:,:) = PRCS(:,:,:) + ZPRCS(:,:,:)/KSPLITR
   PRRS(:,:,:) = PRRS(:,:,:) + ZPRRS(:,:,:)/KSPLITR
   PRSS(:,:,:) = PRSS(:,:,:) + ZPRSS(:,:,:)/KSPLITR
   PRGS(:,:,:) = PRGS(:,:,:) + ZPRGS(:,:,:)/KSPLITR
   IF ( KRR == 7 ) PRHS(:,:,:) = PRHS(:,:,:) + ZPRHS(:,:,:)/KSPLITR
   DO JK = IKTB , IKTE
     ZW(:,:,JK) =ZTSPLITR/(PRHODREF(:,:,JK)* PDZZ(:,:,JK))
   END DO
 ELSE
   IF (OSEDIC) PRCS(:,:,:) = PRCS(:,:,:) + ZPRCS(:,:,:)*ZTSPLITR
   PRRS(:,:,:) = PRRS(:,:,:) + ZPRRS(:,:,:)*ZTSPLITR
   PRSS(:,:,:) = PRSS(:,:,:) + ZPRSS(:,:,:)*ZTSPLITR
   PRGS(:,:,:) = PRGS(:,:,:) + ZPRGS(:,:,:)*ZTSPLITR
   IF ( KRR == 7 ) PRHS(:,:,:) = PRHS(:,:,:) + ZPRHS(:,:,:)*ZTSPLITR
 END IF
 !
 IF (OSEDIC) GSEDIMC(IIB:IIE,IJB:IJE,IKTB:IKTE) =                &
                  PRCS(IIB:IIE,IJB:IJE,IKTB:IKTE)>ZRTMIN(2)
 GSEDIMR(IIB:IIE,IJB:IJE,IKTB:IKTE) =                            &
                  PRRS(IIB:IIE,IJB:IJE,IKTB:IKTE)>ZRTMIN(3)
 GSEDIMI(IIB:IIE,IJB:IJE,IKTB:IKTE) =                            &
                  PRIS(IIB:IIE,IJB:IJE,IKTB:IKTE)>ZRTMIN(4)
 GSEDIMS(IIB:IIE,IJB:IJE,IKTB:IKTE) =                            &
                  PRSS(IIB:IIE,IJB:IJE,IKTB:IKTE)>ZRTMIN(5)
 GSEDIMG(IIB:IIE,IJB:IJE,IKTB:IKTE) =                            &
                  PRGS(IIB:IIE,IJB:IJE,IKTB:IKTE)>ZRTMIN(6)
 IF ( KRR == 7 ) GSEDIMH(IIB:IIE,IJB:IJE,IKTB:IKTE) =            &
                  PRHS(IIB:IIE,IJB:IJE,IKTB:IKTE)>ZRTMIN(7)
!
 IF (OSEDIC) ISEDIMC = COUNTJV( GSEDIMC(:,:,:),IC1(:),IC2(:),IC3(:))
 ISEDIMR = COUNTJV( GSEDIMR(:,:,:),IR1(:),IR2(:),IR3(:))
 ISEDIMI = COUNTJV( GSEDIMI(:,:,:),II1(:),II2(:),II3(:))
 ISEDIMS = COUNTJV( GSEDIMS(:,:,:),IS1(:),IS2(:),IS3(:))
 ISEDIMG = COUNTJV( GSEDIMG(:,:,:),IG1(:),IG2(:),IG3(:))
 IF ( KRR == 7 ) ISEDIMH = COUNTJV( GSEDIMH(:,:,:),IH1(:),IH2(:),IH3(:))
!
!*       2.1   for cloud
!
 IF (OSEDIC) THEN
  ZWSED(:,:,:) = 0.
  IF( JN==1 ) PRCS(:,:,:) = PRCS(:,:,:) * PTSTEP
  IF( ISEDIMC >= 1 ) THEN
    IF ( ISEDIMC .GT. ILENALLOCC ) THEN
      IF ( ILENALLOCC .GT. 0 ) THEN
        DEALLOCATE (ZRCS, ZRHODREFC, ILISTC,ZWLBDC,ZCONC,ZRCT,  &
                    ZZT,ZPRES,ZRAY1D,ZFSEDC1D,ZWLBDA,ZCC )
      END IF
      ILENALLOCC = MAX (IOLDALLOCC, 2*ISEDIMC )
      IOLDALLOCC = ILENALLOCC
      ALLOCATE(ZRCS(ILENALLOCC), ZRHODREFC(ILENALLOCC), ILISTC(ILENALLOCC), &
        ZWLBDC(ILENALLOCC), ZCONC(ILENALLOCC), ZRCT(ILENALLOCC), ZZT(ILENALLOCC), &
        ZPRES(ILENALLOCC), ZRAY1D(ILENALLOCC), ZFSEDC1D(ILENALLOCC), &
        ZWLBDA(ILENALLOCC), ZCC(ILENALLOCC)  )
    END IF
!
    DO JL=1,ISEDIMC
      ZRCS(JL) = PRCS(IC1(JL),IC2(JL),IC3(JL))
      ZRHODREFC(JL) =  PRHODREF(IC1(JL),IC2(JL),IC3(JL))
      ZWLBDC(JL) = ZLBC(IC1(JL),IC2(JL),IC3(JL))
      ZCONC(JL) = ZCONC3D(IC1(JL),IC2(JL),IC3(JL))
      ZRCT(JL) = PRCT(IC1(JL),IC2(JL),IC3(JL))
      ZZT(JL) = PTHT(IC1(JL),IC2(JL),IC3(JL))
      ZPRES(JL) = PPABST(IC1(JL),IC2(JL),IC3(JL))
      ZRAY1D(JL) = ZRAY(IC1(JL),IC2(JL),IC3(JL))
      ZFSEDC1D(JL) = ZFSEDC(IC1(JL),IC2(JL),IC3(JL))
    END DO
!
    ILISTLENC = 0
    DO JL=1,ISEDIMC
     IF( ZRCS(JL) .GT. ZRTMIN(2) ) THEN
       ILISTLENC = ILISTLENC + 1
       ILISTC(ILISTLENC) = JL
     END IF
    END DO
       DO JJ = 1, ILISTLENC
          JL = ILISTC(JJ)
          IF (ZRCS(JL) .GT. ZRTMIN(2) .AND. ZRCT(JL) .GT. XRTMIN(2)) THEN
            ZWLBDC(JL) = ZWLBDC(JL) * ZCONC(JL) / (ZRHODREFC(JL) * ZRCT(JL))
            ZWLBDC(JL) = ZWLBDC(JL)**XLBEXC
            ZRAY1D(JL) = ZRAY1D(JL) / ZWLBDC(JL) !! ZRAY : mean diameter=M(1)/2
            ZZT(JL)    = ZZT(JL) * (ZPRES(JL)/XP00)**(XRD/XCPD)
            ZWLBDA(JL) = 6.6E-8*(101325./ZPRES(JL))*(ZZT(JL)/293.15)
            ZCC(JL)    = XCC*(1.+1.26*ZWLBDA(JL)/ZRAY1D(JL)) !! XCC modified for cloud
            ZWSED (IC1(JL),IC2(JL),IC3(JL))= ZRHODREFC(JL)**(-XCEXVT +1 ) *   &
              ZWLBDC(JL)**(-XDC)*ZCC(JL)*ZFSEDC1D(JL) * ZRCS(JL)
          END IF
       END DO
  END IF
       DO JK = IKTB , IKTE
         PRCS(:,:,JK) = PRCS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
       END DO
       IF (PRESENT(PFPR)) THEN
         DO JK = IKTB , IKTE
           PFPR(:,:,JK,2)=ZWSED(:,:,JK)
         ENDDO
       ENDIF
      PINPRC(:,:) = PINPRC(:,:) + ZWSED(:,:,IKB) / XRHOLW / KSPLITR
      IF( JN==KSPLITR ) THEN
        PRCS(:,:,:) = PRCS(:,:,:) * ZINVTSTEP
      END IF
 END IF
!
!*       2.2   for rain
!
  IF( JN==1 ) PRRS(:,:,:) = PRRS(:,:,:) * PTSTEP
  ZWSED(:,:,:) = 0.
  IF( ISEDIMR >= 1 ) THEN
    IF ( ISEDIMR .GT. ILENALLOCR ) THEN
      IF ( ILENALLOCR .GT. 0 ) THEN
        DEALLOCATE (ZRRS, ZRHODREFR, ILISTR)
      END IF
      ILENALLOCR = MAX (IOLDALLOCR, 2*ISEDIMR )
      IOLDALLOCR = ILENALLOCR
      ALLOCATE(ZRRS(ILENALLOCR), ZRHODREFR(ILENALLOCR), ILISTR(ILENALLOCR))
    END IF
!
    DO JL=1,ISEDIMR
      ZRRS(JL) = PRRS(IR1(JL),IR2(JL),IR3(JL))
      ZRHODREFR(JL) =  PRHODREF(IR1(JL),IR2(JL),IR3(JL))
    END DO
!
    ILISTLENR = 0
    DO JL=1,ISEDIMR
     IF( ZRRS(JL) .GT. ZRTMIN(3) ) THEN
       ILISTLENR = ILISTLENR + 1
       ILISTR(ILISTLENR) = JL
     END IF
    END DO
       DO JJ = 1, ILISTLENR
          JL = ILISTR(JJ)
           ZWSED (IR1(JL),IR2(JL),IR3(JL))= XFSEDR  * ZRRS(JL)**XEXSEDR *   &
                                        ZRHODREFR(JL)**(XEXSEDR-XCEXVT)
       END DO
  END IF
       DO JK = IKTB , IKTE
         PRRS(:,:,JK) = PRRS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
       END DO
       IF (PRESENT(PFPR)) THEN
         DO JK = IKTB , IKTE
           PFPR(:,:,JK,3)=ZWSED(:,:,JK)
         ENDDO
       ENDIF
       PINPRR(:,:) = PINPRR(:,:) + ZWSED(:,:,IKB)/XRHOLW/KSPLITR
       PINPRR3D(:,:,:) = PINPRR3D(:,:,:) + ZWSED(:,:,1:IKT)/XRHOLW/KSPLITR
      IF( JN==KSPLITR ) THEN
        PRRS(:,:,:) = PRRS(:,:,:) * ZINVTSTEP
      END IF
!
!*       2.3   for pristine ice
!
  IF( JN==1 ) PRIS(:,:,:) = PRIS(:,:,:) * PTSTEP
  ZWSED(:,:,:) = 0.
  IF( ISEDIMI >= 1 ) THEN
    IF ( ISEDIMI .GT. ILENALLOCI ) THEN
      IF ( ILENALLOCI .GT. 0 ) THEN
        DEALLOCATE (ZRIS, ZRHODREFI, ILISTI)
      END IF
      ILENALLOCI = MAX (IOLDALLOCI, 2*ISEDIMI )
      IOLDALLOCI = ILENALLOCI
      ALLOCATE(ZRIS(ILENALLOCI), ZRHODREFI(ILENALLOCI), ILISTI(ILENALLOCI))
    END IF
!
    DO JL=1,ISEDIMI
      ZRIS(JL) = PRIS(II1(JL),II2(JL),II3(JL))
      ZRHODREFI(JL) =  PRHODREF(II1(JL),II2(JL),II3(JL))
    END DO
!
    ILISTLENI = 0
    DO JL=1,ISEDIMI
     IF( ZRIS(JL) .GT.  MAX(ZRTMIN(4),1.0E-7 )) THEN ! limitation of the McF&H formula
       ILISTLENI = ILISTLENI + 1
       ILISTI(ILISTLENI) = JL
     END IF
    END DO
       DO JJ = 1, ILISTLENI
          JL = ILISTI(JJ)
              ZWSED (II1(JL),II2(JL),II3(JL))= XFSEDI * ZRIS(JL) *  &
                               ZRHODREFI(JL)**(1.0-XCEXVT) * & !    McF&H
                               MAX( 0.05E6,-0.15319E6-0.021454E6* &
                               ALOG(ZRHODREFI(JL)*ZRIS(JL)) )**XEXCSEDI
       END DO
  END IF
       DO JK = IKTB , IKTE
         PRIS(:,:,JK) = PRIS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
       END DO
       IF (PRESENT(PFPR)) THEN
         DO JK = IKTB , IKTE
           PFPR(:,:,JK,4)=ZWSED(:,:,JK)
         ENDDO
       ENDIF
      IF( JN==KSPLITR ) THEN
        PRIS(:,:,:) = PRIS(:,:,:) * ZINVTSTEP
      END IF
!
!*       2.4   for aggregates/snow
!
  IF( JN==1 ) PRSS(:,:,:) = PRSS(:,:,:) * PTSTEP
  ZWSED(:,:,:) = 0.
  IF( ISEDIMS >= 1 ) THEN
    IF ( ISEDIMS .GT. ILENALLOCS ) THEN
      IF ( ILENALLOCS .GT. 0 ) THEN
        DEALLOCATE (ZRSS, ZRHODREFS, ILISTS)
      END IF
      ILENALLOCS = MAX (IOLDALLOCS, 2*ISEDIMS )
      IOLDALLOCS = ILENALLOCS
      ALLOCATE(ZRSS(ILENALLOCS), ZRHODREFS(ILENALLOCS), ILISTS(ILENALLOCS))
    END IF
!
    DO JL=1,ISEDIMS
      ZRSS(JL) = PRSS(IS1(JL),IS2(JL),IS3(JL))
      ZRHODREFS(JL) =  PRHODREF(IS1(JL),IS2(JL),IS3(JL))
    END DO
!
    ILISTLENS = 0
    DO JL=1,ISEDIMS
     IF( ZRSS(JL) .GT. ZRTMIN(5) ) THEN
       ILISTLENS = ILISTLENS + 1
       ILISTS(ILISTLENS) = JL
     END IF
    END DO
       DO JJ = 1, ILISTLENS
          JL = ILISTS(JJ)
             ZWSED (IS1(JL),IS2(JL),IS3(JL))= XFSEDS * ZRSS(JL)**XEXSEDS *  &
                                        ZRHODREFS(JL)**(XEXSEDS-XCEXVT)
       END DO
  END IF
       DO JK = IKTB , IKTE
         PRSS(:,:,JK) = PRSS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
       END DO
       IF (PRESENT(PFPR)) THEN
         DO JK = IKTB , IKTE
           PFPR(:,:,JK,5)=ZWSED(:,:,JK)
         ENDDO
       ENDIF
      PINPRS(:,:) = PINPRS(:,:) + ZWSED(:,:,IKB)/XRHOLW/KSPLITR
      IF( JN==KSPLITR ) THEN
        PRSS(:,:,:) = PRSS(:,:,:) * ZINVTSTEP
      END IF
!
!*       2.5   for graupeln
!
  ZWSED(:,:,:) = 0.
  IF( JN==1 ) PRGS(:,:,:) = PRGS(:,:,:) * PTSTEP
  IF( ISEDIMG >= 1 ) THEN
    IF ( ISEDIMG .GT. ILENALLOCG ) THEN
      IF ( ILENALLOCG .GT. 0 ) THEN
        DEALLOCATE (ZRGS, ZRHODREFG, ILISTG)
      END IF
      ILENALLOCG = MAX (IOLDALLOCG, 2*ISEDIMG )
      IOLDALLOCG = ILENALLOCG
      ALLOCATE(ZRGS(ILENALLOCG), ZRHODREFG(ILENALLOCG), ILISTG(ILENALLOCG))
    END IF
!
    DO JL=1,ISEDIMG
      ZRGS(JL) = PRGS(IG1(JL),IG2(JL),IG3(JL))
      ZRHODREFG(JL) =  PRHODREF(IG1(JL),IG2(JL),IG3(JL))
    END DO
!
    ILISTLENG = 0
    DO JL=1,ISEDIMG
     IF( ZRGS(JL) .GT. ZRTMIN(6) ) THEN
       ILISTLENG = ILISTLENG + 1
       ILISTG(ILISTLENG) = JL
     END IF
    END DO
       DO JJ = 1, ILISTLENG
          JL = ILISTG(JJ)
             ZWSED (IG1(JL),IG2(JL),IG3(JL))= XFSEDG  * ZRGS(JL)**XEXSEDG *   &
                                        ZRHODREFG(JL)**(XEXSEDG-XCEXVT)
       END DO
END IF
       DO JK = IKTB , IKTE
         PRGS(:,:,JK) = PRGS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
       END DO
       IF (PRESENT(PFPR)) THEN
         DO JK = IKTB , IKTE
           PFPR(:,:,JK,6)=ZWSED(:,:,JK)
         ENDDO
       ENDIF
       PINPRG(:,:) = PINPRG(:,:) + ZWSED(:,:,IKB)/XRHOLW/KSPLITR
      IF( JN==KSPLITR ) THEN
        PRGS(:,:,:) = PRGS(:,:,:) * ZINVTSTEP
      END IF
!
!*       2.6   for hail
!
 IF ( KRR == 7 ) THEN
  IF( JN==1 ) PRHS(:,:,:) = PRHS(:,:,:) * PTSTEP
  ZWSED(:,:,:) = 0.
  IF( ISEDIMH >= 1 ) THEN
    IF ( ISEDIMH .GT. ILENALLOCH ) THEN
      IF ( ILENALLOCH .GT. 0 ) THEN
        DEALLOCATE (ZRHS, ZRHODREFH, ILISTH)
      END IF
      ILENALLOCH = MAX (IOLDALLOCH, 2*ISEDIMH )
      IOLDALLOCH = ILENALLOCH
      ALLOCATE(ZRHS(ILENALLOCH), ZRHODREFH(ILENALLOCH), ILISTH(ILENALLOCH))
    END IF
!
    DO JL=1,ISEDIMH
      ZRHS(JL) = PRHS(IH1(JL),IH2(JL),IH3(JL))
      ZRHODREFH(JL) =  PRHODREF(IH1(JL),IH2(JL),IH3(JL))
    END DO
!
    ILISTLENH = 0
    DO JL=1,ISEDIMH
     IF( ZRHS(JL) .GT. ZRTMIN(7) ) THEN
       ILISTLENH = ILISTLENH + 1
       ILISTH(ILISTLENH) = JL
     END IF
    END DO
       DO JJ = 1, ILISTLENH
          JL = ILISTH(JJ)
             ZWSED (IH1(JL),IH2(JL),IH3(JL))= XFSEDH  * ZRHS(JL)**XEXSEDH *   &
                                        ZRHODREFH(JL)**(XEXSEDH-XCEXVT)
       END DO
  END IF
       DO JK = IKTB , IKTE
         PRHS(:,:,JK) = PRHS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
       END DO
       IF (PRESENT(PFPR)) THEN
         DO JK = IKTB , IKTE
           PFPR(:,:,JK,7)=ZWSED(:,:,JK)
         ENDDO
       ENDIF
       PINPRH(:,:) = PINPRH(:,:) + ZWSED(:,:,IKB)/XRHOLW/KSPLITR
        PRHS(:,:,:) = PRHS(:,:,:) * ZINVTSTEP
      END IF
 END IF
!
END DO
!
IF (OSEDIC) THEN
   IF (ILENALLOCC .GT. 0) DEALLOCATE (ZRCS, ZRHODREFC,  &
  ILISTC,ZWLBDC,ZCONC,ZRCT, ZZT,ZPRES,ZRAY1D,ZFSEDC1D, ZWLBDA,ZCC)
END IF
IF (ILENALLOCR .GT. 0 ) DEALLOCATE(ZRHODREFR,ZRRS,ILISTR)
IF (ILENALLOCI .GT. 0 ) DEALLOCATE(ZRHODREFI,ZRIS,ILISTI)
IF (ILENALLOCS .GT. 0 ) DEALLOCATE(ZRHODREFS,ZRSS,ILISTS)
IF (ILENALLOCG .GT. 0 ) DEALLOCATE(ZRHODREFG,ZRGS,ILISTG)
IF (KRR == 7 .AND. (ILENALLOCH .GT. 0 )) DEALLOCATE(ZRHODREFH,ZRHS,ILISTH)
!
!*       2.3     budget storage
!
IF (LBUDGET_RC .AND. OSEDIC) &
                CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7 ,'SEDI_BU_RRC')
IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:)*PRHODJ(:,:,:),8 ,'SEDI_BU_RRR')
IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:)*PRHODJ(:,:,:),9 ,'SEDI_BU_RRI')
IF (LBUDGET_RS) CALL BUDGET (PRSS(:,:,:)*PRHODJ(:,:,:),10,'SEDI_BU_RRS')
IF (LBUDGET_RG) CALL BUDGET (PRGS(:,:,:)*PRHODJ(:,:,:),11,'SEDI_BU_RRG')
IF ( KRR == 7 .AND. LBUDGET_RH) &
                CALL BUDGET (PRHS(:,:,:)*PRHODJ(:,:,:),12,'SEDI_BU_RRH')
!*       2.4  DROPLET DEPOSITION AT THE 1ST LEVEL ABOVE GROUND
!
IF (LDEPOSC) THEN
  GDEP(:,:) = .FALSE.
  GDEP(IIB:IIE,IJB:IJE) =    PRCS(IIB:IIE,IJB:IJE,IKB) >0 
     PRCS(:,:,IKB) = PRCS(:,:,IKB) - XVDEPOSC * PRCT(:,:,IKB) / PDZZ(:,:,IKB)
     PINPRC(:,:) = PINPRC(:,:) + XVDEPOSC * PRCT(:,:,IKB) * PRHODREF(:,:,IKB) /XRHOLW 
     PINDEP(:,:) = XVDEPOSC * PRCT(:,:,IKB) * PRHODREF(:,:,IKB) /XRHOLW 
  END WHERE
END IF
!
!*       2.5     budget storage
!
IF ( LBUDGET_RC .AND. LDEPOSC ) &
   CALL BUDGET (PRCS(:,:,:)*PRHODJ(:,:,:),7 ,'DEPO_BU_RRC')
  END SUBROUTINE RAIN_ICE_SEDIMENTATION_SPLIT
!
!-------------------------------------------------------------------------------
!
 SUBROUTINE RAIN_ICE_SEDIMENTATION_STAT
!
!*      0. DECLARATIONS
!          ------------
!
IMPLICIT NONE
!
!*       0.2  declaration of local variables
!
!

REAL :: ZP1,ZP2,ZH,ZZWLBDA,ZZWLBDC,ZZCC
REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2)) :: ZQP
INTEGER :: JI,JJ,JK
INTEGER :: JCOUNT, JL
INTEGER, DIMENSION(SIZE(PRHODREF,1)*SIZE(PRHODREF,2)) :: I1, I2
REAL, DIMENSION(SIZE(PRHODREF,1),SIZE(PRHODREF,2),SIZE(PRHODREF,3)) :: ZCONC3D !  droplet condensation
!-------------------------------------------------------------------------------
!
!
!*       1. Parameters for cloud sedimentation
!
  IF (OSEDIC) THEN
    ZRAY(:,:,:)   = 0.
    ZLBC(:,:,:)   = XLBC(1)
    ZFSEDC(:,:,:) = XFSEDC(1)
    ZCONC3D(:,:,:)= XCONC_LAND
    ZCONC_TMP(:,:)= XCONC_LAND
    IF (PRESENT(PSEA)) THEN
      ZCONC_TMP(:,:)=PSEA(:,:)*XCONC_SEA+(1.-PSEA(:,:))*XCONC_LAND
      DO JK=IKTB,IKTE
        ZLBC(:,:,JK)   = PSEA(:,:)*XLBC(2)+(1.-PSEA(:,:))*XLBC(1)
        ZFSEDC(:,:,JK) = (PSEA(:,:)*XFSEDC(2)+(1.-PSEA(:,:))*XFSEDC(1))
        ZFSEDC(:,:,JK) = MAX(MIN(XFSEDC(1),XFSEDC(2)),ZFSEDC(:,:,JK))
        ZCONC3D(:,:,JK)= (1.-PTOWN(:,:))*ZCONC_TMP(:,:)+PTOWN(:,:)*XCONC_URBAN
        ZRAY(:,:,JK)   = 0.5*((1.-PSEA(:,:))*GAMMA(XNUC+1.0/XALPHAC)/(GAMMA(XNUC)) + &
                PSEA(:,:)*GAMMA(XNUC2+1.0/XALPHAC2)/(GAMMA(XNUC2)))
      END DO
    ELSE
        ZCONC3D(:,:,:) = XCONC_LAND
        ZRAY(:,:,:)  = 0.5*(GAMMA(XNUC+1.0/XALPHAC)/(GAMMA(XNUC)))
    END IF
    ZRAY(:,:,:)      = MAX(1.,ZRAY(:,:,:))
    ZLBC(:,:,:)      = MAX(MIN(XLBC(1),XLBC(2)),ZLBC(:,:,:))
  ENDIF
  IF (LDEPOSC) PINDEP (:,:) = 0.
!
!*       2.    compute the fluxes
!

ZRTMIN(:)    = XRTMIN(:) * ZINVTSTEP
!
IF (OSEDIC) THEN
  ZPRCS(:,:,:) = 0.0
  ZPRCS(:,:,:) = PRCS(:,:,:)-PRCT(:,:,:)* ZINVTSTEP
  PRCS(:,:,:)  = PRCT(:,:,:)* ZINVTSTEP
END IF
ZPRRS(:,:,:) = 0.0
ZPRSS(:,:,:) = 0.0
ZPRGS(:,:,:) = 0.0
IF ( KRR == 7 ) ZPRHS(:,:,:) = 0.0
!
ZPRRS(:,:,:) = PRRS(:,:,:)-PRRT(:,:,:)* ZINVTSTEP
ZPRSS(:,:,:) = PRSS(:,:,:)-PRST(:,:,:)* ZINVTSTEP
ZPRGS(:,:,:) = PRGS(:,:,:)-PRGT(:,:,:)* ZINVTSTEP
IF ( KRR == 7 ) ZPRHS(:,:,:) = PRHS(:,:,:)-PRHT(:,:,:)* ZINVTSTEP
PRRS(:,:,:)  = PRRT(:,:,:)* ZINVTSTEP
PRSS(:,:,:)  = PRST(:,:,:)* ZINVTSTEP
PRGS(:,:,:)  = PRGT(:,:,:)* ZINVTSTEP
IF ( KRR == 7 ) PRHS(:,:,:)  = PRHT(:,:,:)* ZINVTSTEP
!
IF (OSEDIC) PRCS(:,:,:) = PRCS(:,:,:) + ZPRCS(:,:,:)
PRRS(:,:,:) = PRRS(:,:,:) + ZPRRS(:,:,:)
PRSS(:,:,:) = PRSS(:,:,:) + ZPRSS(:,:,:)
PRGS(:,:,:) = PRGS(:,:,:) + ZPRGS(:,:,:)
IF ( KRR == 7 ) PRHS(:,:,:) = PRHS(:,:,:) + ZPRHS(:,:,:)
DO JK = IKTB , IKTE
  ZW(:,:,JK) =PTSTEP/(PRHODREF(:,:,JK)* PDZZ(:,:,JK) )
END DO

!
!*       2.1   for cloud
!
 IF (OSEDIC) THEN
     PRCS(:,:,:) = PRCS(:,:,:) * PTSTEP
     ZWSED(:,:,:) = 0.
     ZWSEDW1(:,:,:) = 0.
     ZWSEDW2(:,:,:) = 0.

! calculation of P1, P2 and sedimentation flux
     DO JK = IKE , IKB, -1*KKL
       !estimation of q' taking into account incomming ZWSED
       ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)

       JCOUNT=COUNTJV2((PRCS(:,:,JK) > ZRTMIN(2) .AND. PRCT(:,:,JK) > ZRTMIN(2)) .OR. &
                       (ZQP(:,:) > ZRTMIN(2)),I1(:),I2(:))
       DO JL=1, JCOUNT
         JI=I1(JL)
         JJ=I2(JL)
         !calculation of w
         ! mars 2009 : ajout d'un test
         !IF ( PRCS(JI,JJ,JK) > ZRTMIN(2) ) THEN
         IF(PRCS(JI,JJ,JK) > ZRTMIN(2) .AND. PRCT(JI,JJ,JK) > ZRTMIN(2)) THEN
           ZZWLBDA=6.6E-8*(101325./PPABST(JI,JJ,JK))*(PTHT(JI,JJ,JK)/293.15)
           ZZWLBDC=(ZLBC(JI,JJ,JK)*ZCONC3D(JI,JJ,JK)  &
                &/(PRHODREF(JI,JJ,JK)*PRCT(JI,JJ,JK)))**XLBEXC
           ZZCC=XCC*(1.+1.26*ZZWLBDA*ZZWLBDC/ZRAY(JI,JJ,JK)) !! ZCC  : Fall speed
           ZWSEDW1 (JI,JJ,JK)=PRHODREF(JI,JJ,JK)**(-XCEXVT ) *   &
             &  ZZWLBDC**(-XDC)*ZZCC*ZFSEDC(JI,JJ,JK)
         ENDIF
         IF ( ZQP(JI,JJ) > ZRTMIN(2) ) THEN
           ZZWLBDA=6.6E-8*(101325./PPABST(JI,JJ,JK))*(PTHT(JI,JJ,JK)/293.15)
           ZZWLBDC=(ZLBC(JI,JJ,JK)*ZCONC3D(JI,JJ,JK)  &
                &/(PRHODREF(JI,JJ,JK)*ZQP(JI,JJ)))**XLBEXC
           ZZCC=XCC*(1.+1.26*ZZWLBDA*ZZWLBDC/ZRAY(JI,JJ,JK)) !! ZCC  : Fall speed
           ZWSEDW2 (JI,JJ,JK)=PRHODREF(JI,JJ,JK)**(-XCEXVT ) *   &
             &  ZZWLBDC**(-XDC)*ZZCC*ZFSEDC(JI,JJ,JK)
         ENDIF
       ENDDO

       DO JJ = IJB, IJE
         DO JI = IIB, IIE
           ZH=PDZZ(JI,JJ,JK)
           ZP1 = MIN(1., ZWSEDW1(JI,JJ,JK) * PTSTEP / ZH)
           ! mars 2009 : correction : ZWSEDW1 =>  ZWSEDW2
           !IF (ZWSEDW1(JI,JJ,JK) /= 0.) THEN
           IF (ZWSEDW2(JI,JJ,JK) /= 0.) THEN
             ZP2 = MAX(0.,1 -  ZH &
           &  / (PTSTEP*ZWSEDW2(JI,JJ,JK)) )
           ELSE
             ZP2 = 0.
           ENDIF
           ZWSED (JI,JJ,JK)=ZP1*PRHODREF(JI,JJ,JK)*&
           &ZH*PRCS(JI,JJ,JK)&
           &* ZINVTSTEP+ ZP2 * ZWSED (JI,JJ,JK+KKL)
         ENDDO
       ENDDO

     DO JK = IKTB , IKTE
       PRCS(:,:,JK) = PRCS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
     END DO
     IF (PRESENT(PFPR)) THEN
       DO JK = IKTB , IKTE
         PFPR(:,:,JK,2)=ZWSED(:,:,JK)
       ENDDO
     ENDIF

     PINPRC(:,:) = ZWSED(:,:,IKB)/XRHOLW                        ! in m/s
     PRCS(:,:,:) = PRCS(:,:,:) * ZINVTSTEP
 ENDIF

!
!*       2.2   for rain
!

   PRRS(:,:,:) = PRRS(:,:,:) * PTSTEP
   ZWSED(:,:,:) = 0.
   ZWSEDW1(:,:,:) = 0.
   ZWSEDW2(:,:,:) = 0.

! calculation of ZP1, ZP2 and sedimentation flux
   DO JK = IKE , IKB, -1*KKL
     !estimation of q' taking into account incomming ZWSED
     ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)

     JCOUNT=COUNTJV2((PRRS(:,:,JK) > ZRTMIN(3)) .OR. &
                     (ZQP(:,:) > ZRTMIN(3)),I1(:),I2(:))
     DO JL=1, JCOUNT
       JI=I1(JL)
       JJ=I2(JL)
       !calculation of w
       IF ( PRRS(JI,JJ,JK) > ZRTMIN(3) ) THEN
         ZWSEDW1 (JI,JJ,JK)= XFSEDR *PRRS(JI,JJ,JK)**(XEXSEDR-1)* &
         PRHODREF(JI,JJ,JK)**(XEXSEDR-XCEXVT-1)
       ENDIF
       IF ( ZQP(JI,JJ) > ZRTMIN(3) ) THEN
         ZWSEDW2 (JI,JJ,JK)= XFSEDR *(ZQP(JI,JJ))**(XEXSEDR-1)* &
         PRHODREF(JI,JJ,JK)**(XEXSEDR-XCEXVT-1)
       ENDIF
     ENDDO
     DO JJ = IJB, IJE
       DO JI = IIB, IIE
         ZH=PDZZ(JI,JJ,JK)
         ZP1 = MIN(1., ZWSEDW1(JI,JJ,JK) * PTSTEP / ZH )
         IF (ZWSEDW2(JI,JJ,JK) /= 0.) THEN
           ZP2 = MAX(0.,1 -  ZH &
         & / (PTSTEP*ZWSEDW2(JI,JJ,JK)) )
         ELSE
           ZP2 = 0.
         ENDIF
         ZWSED (JI,JJ,JK)=ZP1*PRHODREF(JI,JJ,JK)*&
         &ZH*PRRS(JI,JJ,JK)&
         &* ZINVTSTEP+ ZP2 * ZWSED (JI,JJ,JK+KKL)
       ENDDO
     ENDDO

   DO JK = IKTB , IKTE
     PRRS(:,:,JK) = PRRS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
   ENDDO
   IF (PRESENT(PFPR)) THEN
     DO JK = IKTB , IKTE
       PFPR(:,:,JK,3)=ZWSED(:,:,JK)
     ENDDO
   ENDIF
   PINPRR(:,:) = ZWSED(:,:,IKB)/XRHOLW                        ! in m/s
   PINPRR3D(:,:,:) = ZWSED(:,:,1:IKT)/XRHOLW                        ! in m/s
   PRRS(:,:,:) = PRRS(:,:,:) * ZINVTSTEP

!
!*       2.3   for pristine ice
!

   PRIS(:,:,:) = PRIS(:,:,:) * PTSTEP
   ZWSED(:,:,:) = 0.
   ZWSEDW1(:,:,:) = 0.
   ZWSEDW2(:,:,:) = 0.
! calculation of ZP1, ZP2 and sedimentation flux
   DO JK = IKE , IKB, -1*KKL
     !estimation of q' taking into account incomming ZWSED
     ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)

     JCOUNT=COUNTJV2((PRIS(:,:,JK) > MAX(ZRTMIN(4),1.0E-7 )) .OR. &
                     (ZQP(:,:) > MAX(ZRTMIN(4),1.0E-7 )),I1(:),I2(:))
     DO JL=1, JCOUNT
       JI=I1(JL)
       JJ=I2(JL)
       !calculation of w
       IF ( PRIS(JI,JJ,JK) > MAX(ZRTMIN(4),1.0E-7 ) ) THEN
         ZWSEDW1 (JI,JJ,JK)= XFSEDI *  &
         &  PRHODREF(JI,JJ,JK)**(XCEXVT) * & !    McF&H
         &  MAX( 0.05E6,-0.15319E6-0.021454E6* &
         &  ALOG(PRHODREF(JI,JJ,JK)*PRIS(JI,JJ,JK)) )**XEXCSEDI
       ENDIF
       IF ( ZQP(JI,JJ) > MAX(ZRTMIN(4),1.0E-7 ) ) THEN
         ZWSEDW2 (JI,JJ,JK)= XFSEDI *  &
         &  PRHODREF(JI,JJ,JK)**(XCEXVT) * & !    McF&H
         &  MAX( 0.05E6,-0.15319E6-0.021454E6* &
         &  ALOG(PRHODREF(JI,JJ,JK)*ZQP(JI,JJ)) )**XEXCSEDI
       ENDIF
     ENDDO
     DO JJ = IJB, IJE
       DO JI = IIB, IIE
         ZH=PDZZ(JI,JJ,JK)
         ZP1 = MIN(1., ZWSEDW1(JI,JJ,JK) * PTSTEP / ZH )
         IF (ZWSEDW2(JI,JJ,JK) /= 0.) THEN
           ZP2 = MAX(0.,1 - ZH  &
           &  / (PTSTEP*ZWSEDW2(JI,JJ,JK)) )
         ELSE
           ZP2 = 0.
         ENDIF
         ZWSED (JI,JJ,JK)=ZP1*PRHODREF(JI,JJ,JK)*&
         &ZH*PRIS(JI,JJ,JK)&
         &* ZINVTSTEP+ ZP2 * ZWSED (JI,JJ,JK+KKL)
       ENDDO
     ENDDO

   DO JK = IKTB , IKTE
     PRIS(:,:,JK) = PRIS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
   ENDDO
   IF (PRESENT(PFPR)) THEN
     DO JK = IKTB , IKTE
       PFPR(:,:,JK,4)=ZWSED(:,:,JK)
     ENDDO
   ENDIF

   PRIS(:,:,:) = PRIS(:,:,:) * ZINVTSTEP


!
!*       2.4   for aggregates/snow
!

   PRSS(:,:,:) = PRSS(:,:,:) * PTSTEP
   ZWSED(:,:,:) = 0.
   ZWSEDW1(:,:,:) = 0.
   ZWSEDW2(:,:,:) = 0.

! calculation of ZP1, ZP2 and sedimentation flux
   DO JK = IKE , IKB, -1*KKL
     !estimation of q' taking into account incomming ZWSED
     ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)

     JCOUNT=COUNTJV2((PRSS(:,:,JK) > ZRTMIN(5)) .OR. &
                     (ZQP(:,:) > ZRTMIN(5)),I1(:),I2(:))
     DO JL=1, JCOUNT
       JI=I1(JL)
       JJ=I2(JL)
       !calculation of w
       IF (PRSS(JI,JJ,JK) > ZRTMIN(5) ) THEN
         ZWSEDW1(JI,JJ,JK)=XFSEDS*(PRSS(JI,JJ,JK))**(XEXSEDS-1)*&
         PRHODREF(JI,JJ,JK)**(XEXSEDS-XCEXVT-1)
       ENDIF
       IF ( ZQP(JI,JJ) > ZRTMIN(5) ) THEN
         ZWSEDW2(JI,JJ,JK)=XFSEDS*(ZQP(JI,JJ))**(XEXSEDS-1)*&
         PRHODREF(JI,JJ,JK)**(XEXSEDS-XCEXVT-1)
       ENDIF
     ENDDO
     DO JJ = IJB, IJE
       DO JI = IIB, IIE
         ZH=PDZZ(JI,JJ,JK)
         ZP1 = MIN(1., ZWSEDW1(JI,JJ,JK) * PTSTEP / ZH )
         IF (ZWSEDW2(JI,JJ,JK) /= 0.) THEN
          / (PTSTEP*ZWSEDW2(JI,JJ,JK)) )
         ELSE
           ZP2 = 0.
         ENDIF
         ZWSED (JI,JJ,JK)=ZP1*PRHODREF(JI,JJ,JK)*&
         &ZH*PRSS(JI,JJ,JK)&
         &* ZINVTSTEP+ ZP2 * ZWSED (JI,JJ,JK+KKL)
       ENDDO
     ENDDO

   DO JK = IKTB , IKTE
     PRSS(:,:,JK) = PRSS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
   ENDDO
   IF (PRESENT(PFPR)) THEN
     DO JK = IKTB , IKTE
       PFPR(:,:,JK,5)=ZWSED(:,:,JK)
     ENDDO
   ENDIF

   PINPRS(:,:) = ZWSED(:,:,IKB)/XRHOLW                        ! in m/s
   PRSS(:,:,:) = PRSS(:,:,:) * ZINVTSTEP


!
!*       2.5   for graupeln
!

   PRGS(:,:,:) = PRGS(:,:,:) * PTSTEP
   ZWSED(:,:,:) = 0.
   ZWSEDW1(:,:,:) = 0.
   ZWSEDW2(:,:,:) = 0.

! calculation of ZP1, ZP2 and sedimentation flux
   DO JK = IKE,  IKB, -1*KKL
     !estimation of q' taking into account incomming ZWSED
     ZQP(:,:)=ZWSED(:,:,JK+KKL)*ZW(:,:,JK)

     JCOUNT=COUNTJV2((PRGS(:,:,JK) > ZRTMIN(6)) .OR. &
                     (ZQP(:,:) > ZRTMIN(6)),I1(:),I2(:))
     DO JL=1, JCOUNT
       JI=I1(JL)
       JJ=I2(JL)
       !calculation of w
       IF ( PRGS(JI,JJ,JK) > ZRTMIN(6) ) THEN
         ZWSEDW1 (JI,JJ,JK)= XFSEDG*(PRGS(JI,JJ,JK))**(XEXSEDG-1) * &
                                  PRHODREF(JI,JJ,JK)**(XEXSEDG-XCEXVT-1)
       ENDIF
       IF ( ZQP(JI,JJ) > ZRTMIN(6) ) THEN
         ZWSEDW2 (JI,JJ,JK)= XFSEDG*(ZQP(JI,JJ))**(XEXSEDG-1) * &
                                  PRHODREF(JI,JJ,JK)**(XEXSEDG-XCEXVT-1)
       ENDIF
     ENDDO
     DO JJ = IJB, IJE
       DO JI = IIB, IIE
         ZH=PDZZ(JI,JJ,JK)
         ZP1 = MIN(1., ZWSEDW1(JI,JJ,JK) * PTSTEP / ZH )
         IF (ZWSEDW2(JI,JJ,JK) /= 0.) THEN
           ZP2 = MAX(0.,1 - ZH &
         & / (PTSTEP*ZWSEDW2(JI,JJ,JK)) )
         ELSE
           ZP2 = 0.
         ENDIF
         ZWSED (JI,JJ,JK)=ZP1*PRHODREF(JI,JJ,JK)*&
         &ZH*PRGS(JI,JJ,JK)&
         &* ZINVTSTEP+ ZP2 * ZWSED (JI,JJ,JK+KKL)
       ENDDO
     ENDDO

   DO JK = IKTB , IKTE
         PRGS(:,:,JK) = PRGS(:,:,JK) + ZW(:,:,JK)*(ZWSED(:,:,JK+KKL)-ZWSED(:,:,JK))
   ENDDO
   IF (PRESENT(PFPR)) THEN
     DO JK = IKTB , IKTE
       PFPR(:,:,JK,6)=ZWSED(:,:,JK)
     ENDDO
   ENDIF

   PINPRG(:,:) = ZWSED(:,:,IKB)/XRHOLW                        ! in m/s
   PRGS(:,:,:) = PRGS(:,:,:) * ZINVTSTEP

!
!*       2.6   for hail