mode_compute_updraft.f90

!               3.2 Critical mixed fraction for KK+KKL flux level (ZKIC_F2) and
!                   for bottom of cloudy part (ZKIC), then a mean for the cloudy part
!                   (put also in ZKIC)
!
!                   computation by estimating unknown  
!                   T^mix r_c^mix and r_i^mix from enthalpy^mix and r_w^mix
!                   We determine the zero crossing of the linear curve
!                   evaluating the derivative using ZMIXF=0.1
                
ZKIC_INIT=0.1  ! starting value for critical mixed fraction for CLoudy Part
ZMIXTHL(:) = 300.
ZMIXRT(:) = 0.1



!  Compute thetaV of environment at the bottom of cloudy part
!    and cons then non cons. var. of mixture at the bottom of cloudy part

!   JKLIM computed to avoid KKL(KK-KKL) being < KKL*KKB
JKLIM=KKL*MAX(KKL*(KK-KKL),KKL*KKB)
DO JIJ=IIJB,IIJE
  IF(OTEST(JIJ) .AND. PPART_DRY(JIJ)>0.5) THEN
    ZDZ=ZDZ_STOP(JIJ)-0.5*(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))
    ZTHV(JIJ)= PTHVM(JIJ,KK)+ZCOEFF_PLUS_HALF(JIJ)*ZDZ
    ZMIXTHL(JIJ) = ZKIC_INIT * &
               (PTHLM(JIJ,KK)+ZDZ*(PTHLM(JIJ,KK+KKL)-PTHLM(JIJ,KK))/PDZZ(JIJ,KK+KKL)) + &
               (1. - ZKIC_INIT)*PTHL_UP(JIJ)
    ZMIXRT(JIJ)  = ZKIC_INIT * &
               (PRTM(JIJ,KK)+ZDZ*(PRTM(JIJ,KK+KKL)-PRTM(JIJ,KK))/PDZZ(JIJ,KK+KKL)) +   &
               (1. - ZKIC_INIT)*PRT_UP(JIJ)
  ELSEIF(OTEST(JIJ)) THEN
    ZDZ=0.5*(PZZ(JIJ,KK+KKL)-PZZ(JIJ,KK))-ZDZ_STOP(JIJ)
    ZTHV(JIJ)= PTHVM(JIJ,KK)-ZCOEFF_MINUS_HALF(JIJ)*ZDZ
    ZMIXTHL(JIJ) = ZKIC_INIT * &
               (PTHLM(JIJ,KK)-ZDZ*(PTHLM(JIJ,KK)-PTHLM(JIJ,JKLIM))/PDZZ(JIJ,KK)) + &
               (1. - ZKIC_INIT)*PTHL_UP(JIJ)
    ZMIXRT(JIJ)  = ZKIC_INIT * &
               (PRTM(JIJ,KK)-ZDZ*(PRTM(JIJ,KK)-PRTM(JIJ,JKLIM))/PDZZ(JIJ,KK)) + &
               (1. - ZKIC_INIT)*PRT_UP(JIJ)
  END IF
ENDDO
CALL TH_R_FROM_THL_RT(CST,NEBN,D%NIJT,NEBN%CFRAC_ICE_SHALLOW_MF,ZFRAC_ICE,&
             ZPRE,ZMIXTHL,ZMIXRT,&
             ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
             ZRSATW_ED, ZRSATI_ED,OOCEAN=.FALSE.,&
             PBUF=ZBUF, KB=D%NIJB, KE=D%NIJE)
!$mnh_expand_array(JIJ=IIJB:IIJE)
ZTHVMIX(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))

!  Compute cons then non cons. var. of mixture at the flux level KK+KKL  with initial ZKIC
ZMIXTHL(:) = ZKIC_INIT * 0.5*(PTHLM(:,KK)+PTHLM(:,KK+KKL))+&
                       & (1. - ZKIC_INIT)*PTHL_UP(:)
ZMIXRT(:)  = ZKIC_INIT * 0.5*(PRTM(:,KK)+PRTM(:,KK+KKL))+&
                       & (1. - ZKIC_INIT)*PRT_UP(:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
CALL TH_R_FROM_THL_RT(CST,NEBN,D%NIJT,NEBN%CFRAC_ICE_SHALLOW_MF,ZFRAC_ICE,&
             PPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
             ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
             ZRSATW_ED, ZRSATI_ED,OOCEAN=.FALSE.,&
             PBUF=ZBUF, KB=D%NIJB, KE=D%NIJE)
!$mnh_expand_array(JIJ=IIJB:IIJE)
ZTHVMIX_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
!$mnh_end_expand_array(JIJ=IIJB:IIJE)

!Computation of mean ZKIC over the cloudy part
DO JIJ=IIJB,IIJE
  IF (OTEST(JIJ)) THEN
    ! Compute ZKIC at the bottom of cloudy part
    ! Thetav_up at bottom is equal to Thetav_up at flux level KK
    IF (ABS(PTHV_UP(JIJ)-ZTHVMIX(JIJ))<1.E-10) THEN
      ZKIC(JIJ)=1.
    ELSE
      ZKIC(JIJ) = MAX(0.,PTHV_UP(JIJ)-ZTHV(JIJ))*ZKIC_INIT /  &  
                 (PTHV_UP(JIJ)-ZTHVMIX(JIJ))
    END IF
    ! Compute ZKIC_F2 at flux level KK+KKL
    IF (ABS(ZTHV_UP_F2(JIJ)-ZTHVMIX_F2(JIJ))<1.E-10) THEN
      ZKIC_F2(JIJ)=1.
    ELSE
      ZKIC_F2(JIJ) = MAX(0.,ZTHV_UP_F2(JIJ)-ZTHV_PLUS_HALF(JIJ))*ZKIC_INIT /  &  
                 (ZTHV_UP_F2(JIJ)-ZTHVMIX_F2(JIJ))
    END IF
    !Mean ZKIC over the cloudy part
    ZKIC(JIJ)=MAX(MIN(0.5*(ZKIC(JIJ)+ZKIC_F2(JIJ)),1.),0.)
  END IF
END DO

!               3.3 Integration of PDF
!                   According to Kain and Fritsch (1990), we replace delta Mt
!                   in eq. (7) and (8) using eq. (5). Here we compute the ratio
!                   of integrals without computing delta Me

!Constant PDF
!For this PDF, eq. (5) is delta Me=0.5*delta Mt
DO JIJ=IIJB,IIJE
  IF(OTEST(JIJ)) THEN
    ZEPSI(JIJ) = ZKIC(JIJ)**2. !integration multiplied by 2
    ZDELTA(JIJ) = (1.-ZKIC(JIJ))**2. !idem
  ENDIF
ENDDO

!Triangular PDF
!Calculus must be verified before activating this part, but in this state,
!results on ARM case are almost identical
!For this PDF, eq. (5) is also delta Me=0.5*delta Mt
!WHERE(OTEST(:))
!  !Integration multiplied by 2
!  WHERE(ZKIC<0.5)
!    ZEPSI(:)=8.*ZKIC(:)**3/3.
!    ZDELTA(:)=1.-4.*ZKIC(:)**2+8.*ZKIC(:)**3/3.
!  ELSEWHERE
!    ZEPSI(:)=5./3.-4*ZKIC(:)**2+8.*ZKIC(:)**3/3.
!    ZDELTA(:)=8.*(1.-ZKIC(:))**3/3.
!  ENDWHERE
!ENDWHERE

!               3.4 Computation of PENTR and PDETR
DO JIJ=IIJB,IIJE
  IF(OTEST(JIJ)) THEN
    ZEPSI_CLOUD=MIN(ZDELTA(JIJ), ZEPSI(JIJ))
    PENTR_CLD(JIJ) = (1.-PPART_DRY(JIJ))*ZCOEFFMF_CLOUD*PRHODREF(JIJ)*ZEPSI_CLOUD
    PDETR_CLD(JIJ) = (1.-PPART_DRY(JIJ))*ZCOEFFMF_CLOUD*PRHODREF(JIJ)*ZDELTA(JIJ)
    PENTR(JIJ) = PENTR(JIJ)+PENTR_CLD(JIJ)
    PDETR(JIJ) = PDETR(JIJ)+PDETR_CLD(JIJ)
  ELSE
    PENTR_CLD(JIJ) = 0.
    PDETR_CLD(JIJ) = 0.
  ENDIF
ENDDO

END SUBROUTINE COMPUTE_ENTR_DETR
INCLUDE "th_r_from_thl_rt.func.h"
INCLUDE "compute_frac_ice.func.h"
END SUBROUTINE COMPUTE_UPDRAFT
END MODULE MODE_COMPUTE_UPDRAFT