phys/module_bl_mfshconvpbl.F

   1      MODULE MODULE_BL_MFSHCONVPBL
   2
   3
   4
   5
   6 USE MODULE_MODEL_CONSTANTS
   7
   8 REAL,PARAMETER :: XG      = 9.80665
   9
  10 REAL,PARAMETER :: XP00= 1.E5               ! Reference pressure
  11 !
  12 !
  13 REAL,PARAMETER :: XMD= 28.9644E-3
  14 REAL,PARAMETER :: XMV= 18.0153E-3            ! Molar mass of dry air and molar mass of vapor
  15 REAL,PARAMETER :: XRD=R_D
  16 REAL,PARAMETER :: XRV=R_V        ! Gaz constant for dry air, gaz constant for vapor
  17 REAL,PARAMETER :: XCPD=7.* XRD /2.
  18 REAL,PARAMETER :: XCPV=4.* XRV          ! Cpd (dry air), Cpv (vapor)
  19 REAL,PARAMETER :: XCL= 4.218E+3
  20 REAL,PARAMETER :: XCI= 2.106E+3      ! Cl (liquid), Ci (ice)
  21 REAL,PARAMETER :: XTT= 273.16        ! Triple point temperature
  22 REAL,PARAMETER :: XLVTT=2.5008E+6              ! Vaporization heat constant
  23 REAL,PARAMETER :: XLSTT=2.8345E+6   ! Sublimation heat constant
  24                                !  temperature
  25 REAL,PARAMETER :: XGAMW = (XCL - XCPV) / XRV! Constants for saturation vapor
  26 REAL,PARAMETER :: XBETAW= (XLVTT/XRV) + (XGAMW * XTT)
  27 !The use of intrinsics in an initialization expressions is a F2003 feature.
  28 !For backward compatibility, hard coded Log(644.11) & Log(XTT) here
  29 !REAL,PARAMETER :: XALPW= LOG(611.14) + (XBETAW /XTT) + (XGAMW *LOG(XTT))
  30 REAL,PARAMETER :: LOG_611_14 = 6.415326
  31 REAL,PARAMETER :: LOG_XTT    = 5.610058
  32 REAL,PARAMETER :: XALPW= LOG_611_14 + (XBETAW /XTT) + (XGAMW *LOG_XTT)
  33                                 !  pressure  function
  34 REAL,PARAMETER :: XGAMI  = (XCI - XCPV) / XRV
  35 REAL,PARAMETER :: XBETAI = (XLSTT/XRV) + (XGAMI * XTT)
  36 !REAL,PARAMETER :: XALPI  = LOG(611.14) + (XBETAI /XTT) + (XGAMI *LOG(XTT))
  37 REAL,PARAMETER :: XALPI  = LOG_611_14 + (XBETAI /XTT) + (XGAMI *LOG_XTT)
  38 REAL,PARAMETER :: XLINI = 0.32
  39
  40
  41 REAL, PARAMETER :: XALP_PERT   = 0.3  ! coefficient for the perturbation of
  42                    ! theta_l and r_t at the first level of
  43                    ! the updraft
  44 REAL, PARAMETER ::XABUO       = 1.   ! coefficient of the buoyancy term in the w_up equation
  45 REAL, PARAMETER ::XBENTR      = 1.   ! coefficient of the entrainment term in thew_up equation
  46 REAL, PARAMETER ::XBDETR      = 0.   ! coefficient of the detrainment term in thew_up equation
  47
  48
  49 REAL, PARAMETER ::XCMF        = 0.065! coefficient for the mass flux at the firstlevel 0.065
  50                                 ! of the updraft (closure) XCMF        = 0.065
  51 REAL, PARAMETER ::XENTR_DRY   = 0.55 ! coefficient for entrainment in dry part XENTR_DRY   = 0.55
  52 REAL, PARAMETER ::XDETR_DRY   = 10.  ! coefficient for detrainment in dry part XDETR_DRY   = 10.
  53 REAL, PARAMETER ::XDETR_LUP   = 1.0   !  coefficient for detrainment in dry part XDETR_LUP   = 1.
  54
  55 REAL, PARAMETER ::XENTR_MF    = 0.035! entrainment constant (m/Pa) = 0.2 (m)  XENTR_MF    = 0.035
  56 REAL, PARAMETER ::XCRAD_MF    = 50.  ! cloud radius in cloudy part
  57 REAL, PARAMETER ::XKCF_MF     = 2.75 ! coefficient for cloud fraction
  58 REAL, PARAMETER ::XKRC_MF     = 1.   ! coefficient for convective rc
  59 REAL, PARAMETER ::XTAUSIGMF   = 600.
  60 REAL, PARAMETER ::XPRES_UV    = 0.5  ! coefficient for pressure term in wind mixing
  61 !
  62 REAL, PARAMETER ::XFRAC_UP_MAX= 0.33 ! maximum Updraft fraction
  63
  64
  65 !
  66
  67      CONTAINS
  68
  69
  70       SUBROUTINE MFSHCONVPBL (DT,STEPBL,HT,DZ           &
  71                     ,RHO,PMID,PINT,TH,EXNER             &
  72                     ,QV, QC, U, V                       &
  73                     ,HFX, QFX, TKE                      &
  74                     ,RUBLTEN,RVBLTEN,RTHBLTEN           &
  75                     ,RQVBLTEN,RQCBLTEN                  &
  76                     ,IDS,IDE,JDS,JDE,KDS,KDE            &
  77                     ,IMS,IME,JMS,JME,KMS,KME            &
  78                     ,ITS,ITE,JTS,JTE,KTS,KTE,KRR        &
  79                     ,MASSFLUX_EDKF, ENTR_EDKF, DETR_EDKF &
  80                     ,THL_UP, THV_UP, RT_UP, RV_UP       &
  81                     ,RC_UP, U_UP, V_UP, FRAC_UP, RC_MF  &
  82                     ,WTHV,PLM_BL89 )
  83
  84       IMPLICIT NONE
  85 !
  86 !----------------------------------------------------------------------
  87       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
  88      &                     ,IMS,IME,JMS,JME,KMS,KME                    &
  89      &                     ,ITS,ITE,JTS,JTE,KTS,KTE
  90 !
  91       INTEGER,INTENT(IN) :: KRR
  92
  93       INTEGER,INTENT(IN) :: STEPBL
  94
  95       REAL,INTENT(IN) :: DT
  96 !
  97       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT, HFX, QFX
  98 !
  99       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ          &
 100      &                                                     ,EXNER       &
 101      &                                                     ,PMID,PINT   &
 102      &                                                     ,QV,QC,RHO   &
 103      &                                                     ,TH,U,V,TKE
 104 !
 105       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) :: RQCBLTEN,RQVBLTEN       &
 106      &                                        ,RTHBLTEN                                &
 107      &                                        ,RUBLTEN,RVBLTEN
 108
 109       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),OPTIONAL,INTENT(OUT) ::              &
 110      &                                         MASSFLUX_EDKF, ENTR_EDKF, DETR_EDKF &
 111      &                                        ,THL_UP, THV_UP, RT_UP, RV_UP        &
 112      &                                        ,RC_UP, U_UP, V_UP, FRAC_UP
 113
 114       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),OPTIONAL,INTENT(INOUT) ::            &
 115      &                                         RC_MF
 116
 117       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) ::   WTHV
 118       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) ::   PLM_BL89
 119
 120 !
 121 !Local declaration
 122
 123 INTEGER   :: KRRL      ! number of liquid water var.
 124 INTEGER   :: KRRI      ! number of ice water var.
 125 LOGICAL   :: OMIXUV    ! True if mixing of momentum
 126
 127 REAL   :: PIMPL_MF     ! degre of implicitness
 128 REAL   ::  PTSTEP   ! Dynamical timestep
 129 REAL   ::  PTSTEP_MET! Timestep for meteorological variables
 130
 131 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)   :: PZZ       !  Height at the flux point
 132
 133 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)   :: PZZM      !  Height at the mass point
 134
 135 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)   :: PDZZ      ! depth between mass levels
 136
 137 REAL, DIMENSION(ITS:ITE,JTS:JTE)   ::  PSFTH,PSFRV
 138                                             ! normal surface fluxes of theta,rv
 139 !
 140 !    prognostic variables at t- deltat
 141 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PPABSM      ! Pressure at mass point
 142 !REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PPABSF    ! Pressure at flux point
 143 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PEXNM       ! Exner function at t-dt
 144 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PRHODREF    ! dry density of the
 145                                                      ! reference state
 146 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PRHODJ    ! dry density of the
 147
 148 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PTKEM       ! TKE
 149 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE) ::  PUM,PVM     ! momentum
 150
 151 !   thermodynamical variables which are transformed in conservative var.
 152 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     ::  PTHM       ! pot. temp. = PTHLM in turb.f90
 153 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE,KRR) ::  PRM        ! water species
 154
 155 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     ::  PRUS,PRVS,PRTHS
 156 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE,KRR) ::  PRRS
 157
 158 ! For diagnostic output
 159 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PEMF, PENTR, PDETR
 160 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PTHL_UP, PRT_UP, PRV_UP, PRC_UP
 161 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PU_UP, PV_UP, PTHV_UP, PFRAC_UP
 162 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PRC_MF
 163 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: WTHV_MF
 164 REAL, DIMENSION(ITS:ITE,JTS:JTE,KTS:KTE)     :: PLM_MF
 165
 166 INTEGER :: I,J,K  ! loop variables
 167
 168 ! Transform WRF Variable to input of mass flux scheme
 169
 170 DO J=JTS,JTE
 171    DO K=KTS,KTE
 172       DO I=ITS,ITE
 173         IF (K==KTS) PZZ(I,J,K)=0.
 174         PEMF(I,J,K)=0.
 175         PENTR(I,J,K)=0.
 176         PDETR(I,J,K)=0.
 177         PTHL_UP(I,J,K)=0.
 178         PTHV_UP(I,J,K)=0.
 179         PRT_UP(I,J,K)=0.
 180         PRV_UP(I,J,K)=0.
 181         PRC_UP(I,J,K)=0.
 182         PU_UP(I,J,K)=0.
 183         PV_UP(I,J,K)=0.
 184         PFRAC_UP(I,J,K)=0.
 185         WTHV_MF(I,J,K)=0.
 186         PTHM(I,J,K)=TH(I,K,J)
 187         PTKEM(I,J,K)=TKE(I,K,J)
 188         PRM(I,J,K,1)=QV(I,K,J)-RC_MF(I,K,J)
 189         PRM(I,J,K,2)=RC_MF(I,K,J)
 190         PUM(I,J,K)=U(I,K,J)
 191         PVM(I,J,K)=V(I,K,J)
 192         PRHODREF(I,J,K)=RHO(I,K,J)/(1.+QV(I,K,J))
 193         PEXNM(I,J,K)=EXNER(I,K,J)
 194         PPABSM(I,J,K)=PMID(I,K,J)
 195         IF (K/=KTE) THEN
 196             PZZ(I,J,K+1)=PZZ(I,J,K)+DZ(I,K,J)
 197             PZZM(I,J,K)=0.5*(PZZ(I,J,K+1)+PZZ(I,J,K)) ! z at mass point
 198         ELSE
 199             PZZM(I,J,K)=PZZ(I,J,K)+0.5*DZ(I,K-1,J) ! z at mass point
 200         ENDIF
 201         IF (K==KTS) THEN
 202            PDZZ(I,J,K)=2*(PZZM(I,J,K))
 203         ELSE
 204            PDZZ(I,J,K)=PZZM(I,J,K)-PZZM(I,J,K-1)
 205         ENDIF
 206
 207         PRHODJ(I,J,K)=PRHODREF(I,J,K)*DZ(I,K,J)
 208
 209
 210       ENDDO
 211    ENDDO
 212 ENDDO
 213
 214 ! fill the kte+1 level
 215 PTHM(:,:,KTE)=PTHM(:,:,KTE-1)
 216 PTKEM(:,:,KTE)=PTKEM(:,:,KTE-1)
 217 PRM(:,:,KTE,1)=PRM(:,:,KTE-1,1)
 218 PRM(:,:,KTE,2)=PRM(:,:,KTE-1,2)
 219 PUM(:,:,KTE)=PUM(:,:,KTE-1)
 220 PVM(:,:,KTE)=PVM(:,:,KTE-1)
 221 PRHODREF(:,:,KTE)=PRHODREF(:,:,KTE-1)
 222 PEXNM(:,:,KTE)=PEXNM(:,:,KTE-1)
 223 PPABSM(:,:,KTE)=PPABSM(:,:,KTE-1)
 224 PRHODJ(:,:,KTE)=PRHODJ(:,:,KTE-1)
 225
 226 PSFTH(:,:)=HFX(ITS:ITE,JTS:JTE)/(PRHODREF(:,:,KTS)*XCPD)
 227 PSFRV(:,:)=QFX(ITS:ITE,JTS:JTE)/(PRHODREF(:,:,KTS))
 228
 229
 230
 231 ! Assign some variables
 232
 233 OMIXUV=.FALSE.
 234 KRRL=1 !Qc is managed
 235 KRRI=0 !Qi not
 236 PIMPL_MF=0.
 237 PTSTEP=DT*STEPBL
 238 PTSTEP_MET=PTSTEP
 239
 240       CALL MFSHCONVPBL_CORE(KRR,KRRL,KRRI,                            &
 241                 OMIXUV,                                               &
 242                 PIMPL_MF,PTSTEP,PTSTEP_MET,                           &
 243                 PDZZ, PZZ,                                            &
 244                 PRHODJ, PRHODREF,                                     &
 245                 PPABSM, PEXNM,                                        &
 246                 PSFTH,PSFRV,                                          &
 247                 PTHM,PRM,PUM,PVM,PTKEM,                               &
 248                 PRTHS,PRRS,PRUS,PRVS,PEMF, PENTR, PDETR,              &
 249                 PTHL_UP, PRT_UP, PRV_UP, PRC_UP,                      &
 250                 PU_UP, PV_UP, PTHV_UP, PFRAC_UP, PRC_MF, WTHV_MF,PLM_MF  )
 251
 252
 253 DO J=JTS,JTE
 254    DO K=KTS,KTE
 255       DO I=ITS,ITE
 256          RQCBLTEN(I,K,J)=PRRS(I,J,K,2)
 257          RQVBLTEN(I,K,J)=PRRS(I,J,K,1)
 258          RTHBLTEN(I,K,J)=PRTHS(I,J,K)
 259          RUBLTEN(I,K,J)=PRUS(I,J,K)
 260          RVBLTEN(I,K,J)=PRVS(I,J,K)
 261          WTHV(I,K,J)=WTHV_MF(I,J,K)
 262          PLM_BL89(I,K,J)=PLM_MF(I,J,K)
 263       ENDDO
 264    ENDDO
 265 ENDDO
 266
 267       IF ( PRESENT(MASSFLUX_EDKF) ) THEN
 268          DO J=JTS,JTE
 269             DO K=KTS,KTE
 270                DO I=ITS,ITE
 271                   MASSFLUX_EDKF(I,K,J)=PEMF(I,J,K)
 272                   ENTR_EDKF(I,K,J)=PENTR(I,J,K)
 273                   DETR_EDKF(I,K,J)=PDETR(I,J,K)
 274                   THL_UP(I,K,J)=PTHL_UP(I,J,K)
 275                   THV_UP(I,K,J)=PTHV_UP(I,J,K)
 276                   RT_UP(I,K,J)=PRT_UP(I,J,K)
 277                   RV_UP(I,K,J)=PRV_UP(I,J,K)
 278                   RC_UP(I,K,J)=PRC_UP(I,J,K)
 279                   U_UP(I,K,J)=PU_UP(I,J,K)
 280                   V_UP(I,K,J)=PV_UP(I,J,K)
 281                   FRAC_UP(I,K,J)=PFRAC_UP(I,J,K)
 282                   RC_MF(I,K,J)=PRC_MF(I,J,K)
 283                ENDDO
 284             ENDDO
 285          ENDDO
 286       ENDIF
 287
 288
 289 END SUBROUTINE MFSHCONVPBL
 290
 291 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 292 !   WRAPPER from WRF to MASS FLUX SCHEME
 293 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 294
 295       SUBROUTINE MFSHCONVPBL_CORE(KRR,KRRL,KRRI,                            &
 296                 OMIXUV,                       &
 297                 PIMPL_MF,PTSTEP,PTSTEP_MET,                 &
 298                 PDZZ, PZZ,                                            &
 299                 PRHODJ, PRHODREF,                                     &
 300                 PPABSM, PEXNM,                                        &
 301                 PSFTH,PSFRV,                                          &
 302                 PTHM,PRM,PUM,PVM,PTKEM,                               &
 303                 PRTHS,PRRS,PRUS,PRVS,  PEMF, PENTR, PDETR,            &
 304                 PTHL_UP, PRT_UP, PRV_UP, PRC_UP,                      &
 305                 PU_UP, PV_UP, PTHV_UP, PFRAC_UP,  PRC_MF,             &
 306                 PFLXZTHVMF,PLM       )
 307 !!
 308 !!****  *MFSHCONVPBL_CORE* - Interfacing routine
 309 !!
 310
 311 !! --------------------------------------------------------------------------
 312 !
 313
 314
 315 IMPLICIT NONE
 316
 317 INTEGER,                INTENT(IN)   :: KRR          ! number of moist var.
 318 INTEGER,                INTENT(IN)   :: KRRL         ! number of liquid water var.
 319 INTEGER,                INTENT(IN)   :: KRRI         ! number of ice water var.
 320 LOGICAL,                INTENT(IN)   :: OMIXUV    ! True if mixing of momentum
 321
 322 REAL,                   INTENT(IN)   :: PIMPL_MF     ! degre of implicitness
 323 REAL,                 INTENT(IN)     ::  PTSTEP   ! Dynamical timestep
 324 REAL,                 INTENT(IN)     ::  PTSTEP_MET! Timestep for meteorological variables
 325 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PZZ         ! Height of flux point
 326 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PDZZ        ! Metric coefficients
 327 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PRHODJ      ! dry density * Grid size
 328 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PRHODREF    ! dry density of the
 329                                                      ! reference state
 330 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PPABSM      ! Pressure at time t-1
 331 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PEXNM       ! Exner function at t-dt
 332
 333 REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTH,PSFRV ! normal surface fluxes of theta and Rv
 334 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PTHM        ! Theta at t-dt
 335 REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM         ! water var. at t-dt
 336 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PUM,PVM     ! wind components at t-dt
 337 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKEM       ! tke at t-dt
 338
 339
 340 REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PRUS,PRVS,PRTHS ! Meso-NH sources
 341 REAL, DIMENSION(:,:,:,:), INTENT(OUT) ::  PRRS
 342
 343 REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PEMF, PENTR, PDETR
 344 REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PTHL_UP, PRT_UP, PRV_UP, PRC_UP
 345 REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PU_UP, PV_UP, PTHV_UP, PFRAC_UP
 346 REAL, DIMENSION(:,:,:), INTENT(INOUT) ::  PRC_MF
 347 REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PFLXZTHVMF
 348 REAL, DIMENSION(:,:,:), INTENT(OUT) ::  PLM
 349
 350
 351 !
 352 !                     0.2  Declaration of local variables
 353 !
 354 REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2),SIZE(PTHM,3)) ::     &
 355           ZEXN,ZCPH,                                      &
 356           PRV,PRL,PTH,                                    &
 357           ZTM,                                            &  ! Temperature at t-dt
 358           ZLVOCPEXN,                                      &  !
 359           ZCF_MF,                                         &
 360           ZLSOCPEXN,                                      &  !
 361           ZAMOIST,                                        &  !
 362           ZATHETA,                                        &  !
 363           ZTHLM,                                          &  !
 364           ZRTM,                                           &  !
 365           ZTHVM,ZTHVREF,ZUMM,ZVMM,                        &  !
 366           ZRI_UP,ZW_UP,                                   & !
 367           ZEMF_O_RHODREF,                                 & ! entrainment/detrainment
 368           ZTHLDT,ZRTDT,ZUDT,ZVDT,                         & ! tendencies
 369           ZFLXZTHMF,ZFLXZRMF,ZFLXZUMF,ZFLXZVMF   ! fluxes
 370
 371
 372 INTEGER,DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2))  :: IKLCL,IKETL,IKCTL
 373
 374 INTEGER :: IKU, IKB, IKE
 375 INTEGER :: JI,JJ,JK,JSV                          ! Loop counters
 376 INTEGER :: IRESP    ! error code
 377
 378
 379 !------------------------------------------------------------------------
 380
 381 !!! 1. Initialisation
 382
 383 ! Internal Domain
 384 IKU=SIZE(PTHM,3)
 385 IKB=1              ! Modif WRF JP
 386 IKE=IKU-1
 387
 388
 389 ! number of scalar var
 390
 391   ZUMM=PUM   !Modif WRF JP
 392   ZVMM=PVM   !Modif WRF JP
 393
 394 ! Thermodynamics functions
 395
 396   CALL COMPUTE_FUNCTION_THERMO_MF( KRR,KRRL,KRRI,                   &
 397                                    PTHM,PRM,PEXNM,PPABSM,           &
 398                                    ZTM,ZLVOCPEXN,ZLSOCPEXN,         &
 399                                    ZAMOIST,ZATHETA                  )
 400
 401
 402
 403 ! Conservative variables at t-dt
 404
 405   CALL THL_RT_FROM_TH_R_MF( KRR,KRRL,KRRI,                         &
 406                                   PTHM, PRM, ZLVOCPEXN, ZLSOCPEXN,       &
 407                                   ZTHLM, ZRTM                            )
 408
 409
 410
 411 ! Virtual potential temperature at t-dt
 412
 413 ZTHVM(:,:,:) = PTHM(:,:,:)*((1.+XRV / XRD *PRM(:,:,:,1))/(1.+ZRTM(:,:,:)))
 414
 415
 416        ZTHVREF=XG/ZTHVM
 417        CALL BL89(PZZ,PDZZ,ZTHVREF,ZTHLM,KRR, &
 418                  PRM,PTKEM,PLM)
 419
 420 !!! 2. Compute updraft
 421
 422
 423       CALL UPDRAFT_SOPE (KRR,KRRL,KRRI,OMIXUV,                     &
 424                          PZZ,PDZZ,PSFTH,PSFRV,PPABSM,PRHODREF,     &
 425                          PTKEM,PTHM,PRM,ZTHLM,ZRTM,ZUMM,ZVMM,      &
 426                          PTHL_UP,PRT_UP,PRV_UP,PU_UP,PV_UP,        &
 427                          PRC_UP,ZRI_UP,PTHV_UP,ZW_UP,PFRAC_UP,PEMF,&
 428                          PDETR,PENTR,IKLCL,IKETL,IKCTL )
 429
 430 !!! 3. Compute fluxes of conservative variables and their divergence = tendency
 431
 432        ZEMF_O_RHODREF=PEMF/PRHODREF
 433        CALL MF_TURB(OMIXUV, PIMPL_MF, PTSTEP,PTSTEP_MET,       &
 434                 PDZZ, PRHODJ, ZTHLM,ZTHVM,ZRTM,ZUMM,ZVMM,           &
 435                 ZTHLDT,ZRTDT,ZUDT,ZVDT,                            &
 436                 ZEMF_O_RHODREF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP,&
 437                 ZFLXZTHMF,PFLXZTHVMF,ZFLXZRMF,ZFLXZUMF,ZFLXZVMF   )
 438
 439
 440 !!! 5. Compute diagnostic convective cloud fraction and content
 441 ! ! !  ONLY liquid cloud implemented (yet)
 442
 443       CALL COMPUTE_MF_CLOUD(KRRL,ZTHLM,PRC_UP,PFRAC_UP,PDZZ,IKLCL,  &
 444                             PRC_MF,ZCF_MF                      )
 445
 446
 447 !!! 6. Compute tendency terms for pronostic variables
 448
 449 ZEXN(:,:,:)=(PPABSM(:,:,:)/XP00) ** (XRD/XCPD)
 450 !
 451
 452
 453 PRV(:,:,:)=PRM(:,:,:,1)-PRC_MF(:,:,:)
 454 PRL(:,:,:)=PRC_MF(:,:,:)
 455       !       2.1 Cph
 456 ZCPH(:,:,:)=XCPD+ XCPV * PRV(:,:,:)+ XCL * PRL(:,:,:)
 457
 458 PTH(:,:,:)=(ZTHLM(:,:,:)+ZTHLDT(:,:,:))+(XLVTT/(ZCPH*ZEXN(:,:,:))*PRL(:,:,:))
 459
 460 PRTHS(:,:,:)  = ZTHLDT(:,:,:)
 461 PRTHS(:,:,:)  = (PTH(:,:,:)-PTHM(:,:,:))/PTSTEP_MET
 462 PRRS(:,:,:,2) = (PRC_MF-PRM(:,:,:,2))/PTSTEP_MET
 463 PRRS(:,:,:,1) = ZRTDT(:,:,:)-PRRS(:,:,:,2)
 464
 465 PRTHS(:,:,:)  = ZTHLDT(:,:,:)
 466 PRRS(:,:,:,1) = ZRTDT(:,:,:)
 467 PRRS(:,:,:,2) = 0
 468 PRUS(:,:,:)   = ZUDT(:,:,:)
 469 PRVS(:,:,:)   = ZVDT(:,:,:)
 470
 471
 472 END SUBROUTINE MFSHCONVPBL_CORE
 473
 474
 475 !     ###################################################################
 476       SUBROUTINE COMPUTE_BL89_ML(PDZZ2D, &
 477              PTKEM2D,PG_O_THVREF2D,PVPT,KK,OUPORDN,PLWORK)
 478 !     ###################################################################
 479 !!
 480 !!     COMPUTE_BL89_ML routine to:
 481 !!
 482 !-------------------------------------------------------------------------------
 483 !
 484 !*       0.    DECLARATIONS
 485 !
 486 !              ------------
 487 !USE MODD_CST
 488 !
 489 !
 490 IMPLICIT NONE
 491 !
 492 !          0.1 arguments
 493 !
 494 REAL, DIMENSION(:,:),   INTENT(IN)  :: PDZZ2D
 495 REAL, DIMENSION(:,:),   INTENT(IN)  :: PTKEM2D
 496 REAL, DIMENSION(:,:),   INTENT(IN)  :: PG_O_THVREF2D
 497 REAL, DIMENSION(:,:),   INTENT(IN)  :: PVPT
 498 INTEGER,                INTENT(IN)  :: KK
 499 LOGICAL,                INTENT(IN)  :: OUPORDN
 500 REAL, DIMENSION(:),     INTENT(OUT) :: PLWORK
 501
 502 !          0.2 Local variable
 503 !
 504 REAL, DIMENSION(SIZE(PTKEM2D,1)) :: ZLWORK1,ZLWORK2 ! Temporary mixing length
 505 REAL, DIMENSION(SIZE(PTKEM2D,1)) :: ZINTE,ZPOTE     ! TKE and potential energy
 506                                                     ! between 2 levels
 507 INTEGER :: IKB,IKE
 508 !
 509 REAL, DIMENSION(SIZE(PTKEM2D,1),SIZE(PTKEM2D,2)) :: ZDELTVPT,ZHLVPT
 510                       !Virtual Potential Temp at Half level and DeltaThv between
 511                       !2 levels
 512 REAL, DIMENSION(SIZE(PTKEM2D,1)) :: ZTH! Potential Temp
 513
 514 INTEGER :: IIJU,IKU             !Internal Domain
 515 INTEGER :: J1D                  !horizontal loop counter
 516 INTEGER :: JKK                  !loop counters
 517 INTEGER :: JRR                  !moist loop counter
 518 INTEGER :: JIJK                 !loop counters
 519 REAL    :: ZTEST,ZTEST0,ZTESTM  !test for vectorization
 520 !-------------------------------------------------------------------------------------
 521 !
 522 !*       1.    INITIALISATION
 523 !              --------------
 524 IIJU=SIZE(PTKEM2D,1)
 525 !
 526 IKB = 1 !Modif WRF JP
 527 IKE = SIZE(PTKEM2D,2)-1 !Modif WRF JP
 528 IKU = SIZE(PTKEM2D,2)
 529
 530 ZDELTVPT(:,2:IKU)=PVPT(:,2:IKU)-PVPT(:,1:IKU-1)
 531 ZDELTVPT(:,1)=0.
 532
 533 ! to prevent division by zero
 534 WHERE (ABS(ZDELTVPT(:,:))<1.E-10)
 535   ZDELTVPT(:,:)=1.E-10
 536 END WHERE
 537 !
 538 ZHLVPT(:,2:IKU)= 0.5 * ( PVPT(:,2:IKU)+PVPT(:,1:IKU-1) )
 539 ZHLVPT(:,1)    = PVPT(:,1)
 540
 541 !
 542 !
 543 !
 544 !*       2.    CALCULATION OF THE UPWARD MIXING LENGTH
 545 !              ---------------------------------------
 546 !
 547
 548 IF (OUPORDN.EQV..TRUE.) THEN
 549  ZINTE(:)=PTKEM2D(:,KK)
 550  PLWORK=0.
 551  ZLWORK1=0.
 552  ZLWORK2=0.
 553  ZTESTM=1.
 554  ZTH(:)=PVPT(:,KK)
 555  DO JKK=KK+1,IKE
 556     IF(ZTESTM > 0.) THEN
 557       ZTESTM=0
 558       DO J1D=1,IIJU
 559         ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
 560         ZPOTE(J1D) = ZTEST0*(PG_O_THVREF2D(J1D,KK)      *      &
 561             (ZHLVPT(J1D,JKK) - ZTH(J1D)))  * PDZZ2D(J1D,JKK) !particle keeps its temperature
 562         ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
 563         ZTESTM=ZTESTM+ZTEST0
 564         ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
 565         !ZLWORK2 jump of the last reached level
 566         ZLWORK2(J1D)=        ( - PG_O_THVREF2D(J1D,KK) *                     &
 567             (  PVPT(J1D,JKK-1) - ZTH(J1D) )                              &
 568           + SQRT (ABS(                                                       &
 569             ( PG_O_THVREF2D(J1D,KK) * (PVPT(J1D,JKK-1) - ZTH(J1D)) )**2  &
 570             + 2. * ZINTE(J1D) * PG_O_THVREF2D(J1D,KK)                        &
 571                  * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) ))    ) /             &
 572         ( PG_O_THVREF2D(J1D,KK) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
 573       !
 574         PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+  &
 575                                     (1-ZTEST)*ZLWORK2(J1D))
 576         ZINTE(J1D) = ZINTE(J1D) - ZPOTE(J1D)
 577       END DO
 578     ENDIF
 579   END DO
 580 ENDIF
 581 !!
 582 !*       2.    CALCULATION OF THE DOWNWARD MIXING LENGTH
 583 !              ---------------------------------------
 584 !
 585
 586 IF (OUPORDN.EQV..FALSE.) THEN
 587  ZINTE(:)=PTKEM2D(:,KK)
 588  PLWORK=0.
 589  ZLWORK1=0.
 590  ZLWORK2=0.
 591  ZTESTM=1.
 592  ZTH(:)=PVPT(:,KK)
 593  DO JKK=KK,IKB,-1
 594     IF(ZTESTM > 0.) THEN
 595       ZTESTM=0
 596       DO J1D=1,IIJU
 597         ZTEST0=0.5+SIGN(0.5,ZINTE(J1D))
 598         ZPOTE(J1D) = -ZTEST0*(PG_O_THVREF2D(J1D,KK)      *      &
 599             (ZHLVPT(J1D,JKK) - ZTH(J1D)))  * PDZZ2D(J1D,JKK) !particle keeps its temperature
 600         ZTEST =0.5+SIGN(0.5,ZINTE(J1D)-ZPOTE(J1D))
 601         ZTESTM=ZTESTM+ZTEST0
 602         ZLWORK1(J1D)=PDZZ2D(J1D,JKK)
 603         ZLWORK2(J1D)=        ( + PG_O_THVREF2D(J1D,KK) *                     &
 604             (  PVPT(J1D,JKK) - ZTH(J1D) )                              &
 605           + SQRT (ABS(                                                       &
 606             ( PG_O_THVREF2D(J1D,KK) * (PVPT(J1D,JKK) - ZTH(J1D)) )**2  &
 607             + 2. * ZINTE(J1D) * PG_O_THVREF2D(J1D,KK)                        &
 608                  * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) ))    ) /             &
 609         ( PG_O_THVREF2D(J1D,KK) * ZDELTVPT(J1D,JKK) / PDZZ2D(J1D,JKK) )
 610       !
 611         PLWORK(J1D)=PLWORK(J1D)+ZTEST0*(ZTEST*ZLWORK1(J1D)+  &
 612                                     (1-ZTEST)*ZLWORK2(J1D))
 613         ZINTE(J1D) = ZINTE(J1D) - ZPOTE(J1D)
 614       END DO
 615     ENDIF
 616   END DO
 617 ENDIF
 618
 619 END SUBROUTINE COMPUTE_BL89_ML
 620 !
 621
 622 !
 623 !         #############################################################
 624           SUBROUTINE COMPUTE_ENTR_DETR(OTEST,OTESTLCL,&
 625                             HFRAC_ICE,PFRAC_ICE,KK,PPABSM,PZZ,PDZZ,&
 626                             PTHVM,PTHLM,PRTM,PW_UP2,&
 627                             PTHL_UP,PRT_UP,PLUP,&
 628                             PENTR,PDETR,PBUO_INTEG)
 629 !         #############################################################
 630
 631 !!
 632 !!***COMPUTE_ENTR_DETR* - calculates caracteristics of the updraft or downdraft
 633 !!                       using model of the EDMF scheme
 634 !!
 635 !!
 636 !! --------------------------------------------------------------------------
 637 !
 638
 639
 640
 641 IMPLICIT NONE
 642 !
 643 !
 644 !*                    1.1  Declaration of Arguments
 645 !
 646 !
 647
 648 LOGICAL,DIMENSION(:),INTENT(INOUT)  :: OTEST ! test to see if updraft is running
 649 LOGICAL,DIMENSION(:),INTENT(INOUT)  :: OTESTLCL !test of condensation
 650 CHARACTER*1                         :: HFRAC_ICE ! frac_ice can be compute using
 651                                               ! Temperature (T) or prescribed
 652                                               ! (Y)
 653 REAL, DIMENSION(:), INTENT(INOUT)   :: PFRAC_ICE ! if frac_ice is prescribed
 654 INTEGER,            INTENT(IN)      :: KK  ! level where E and D are computed
 655 !
 656 !    prognostic variables at t- deltat
 657 !
 658 REAL, DIMENSION(:,:),   INTENT(IN) ::  PPABSM      ! Pressure at time t-1
 659 REAL, DIMENSION(:,:),   INTENT(IN) ::  PZZ       !  Height at the flux point
 660 REAL, DIMENSION(:,:),   INTENT(IN) ::  PDZZ       !  metrics coefficient
 661 REAL, DIMENSION(:,:),   INTENT(IN) ::  PTHVM      ! ThetaV environment
 662
 663 !
 664 !   thermodynamical variables which are transformed in conservative var.
 665 !
 666 REAL, DIMENSION(:),   INTENT(IN)     ::  PTHLM     ! Thetal
 667 REAL, DIMENSION(:),   INTENT(IN)     ::  PRTM      ! total mixing ratio
 668 REAL, DIMENSION(:,:), INTENT(INOUT)  ::  PW_UP2    ! Vertical velocity^2
 669 REAL, DIMENSION(:),   INTENT(IN)     ::  PTHL_UP,PRT_UP  ! updraft properties
 670 REAL, DIMENSION(:),   INTENT(IN)     ::  PLUP      ! LUP compute from the ground
 671 REAL, DIMENSION(:),   INTENT(INOUT)  ::  PENTR     ! Mass flux entrainment of the updraft
 672 REAL, DIMENSION(:),   INTENT(INOUT)  ::  PDETR     ! Mass flux detrainment of the updraft
 673 REAL, DIMENSION(:),   INTENT(INOUT)  ::  PBUO_INTEG! Integral Buoyancy
 674 !
 675 !
 676 !                       1.2  Declaration of local variables
 677 !
 678 !
 679
 680 ! Variables for cloudy part
 681
 682 REAL, DIMENSION(SIZE(PTHLM))   :: ZKIC           ! fraction of env. mass in the muxtures
 683 REAL, DIMENSION(SIZE(PTHLM))   :: ZEPSI,ZDELTA   ! factor entrainment detrainment
 684 REAL, DIMENSION(SIZE(PTHLM))   :: ZEPSI_CLOUD    ! factor entrainment detrainment
 685 REAL, DIMENSION(SIZE(PTHLM))   :: ZCOEFFMF_CLOUD ! factor for compputing entr. detr.
 686
 687 REAL, DIMENSION(SIZE(PTHLM))   :: ZMIXTHL,ZMIXRT ! Thetal and rt in the mixtures
 688 !
 689 REAL, DIMENSION(SIZE(PTHLM))   :: ZTHMIX,ZTHVMIX ! Theta and Thetav  of mixtures
 690 REAL, DIMENSION(SIZE(PTHLM))   :: ZRVMIX,ZRCMIX,ZRIMIX ! mixing ratios in mixtures
 691
 692 REAL, DIMENSION(SIZE(PTHLM))   :: ZTHMIX_F2     ! Theta and Thetav  of mixtures
 693 REAL, DIMENSION(SIZE(PTHLM))   :: ZRVMIX_F2,ZRCMIX_F2,ZRIMIX_F2 ! mixing ratios in mixtures
 694
 695 REAL, DIMENSION(SIZE(PTHLM))   :: ZTHV_UP       ! thvup at mass point kk
 696
 697
 698 REAL, DIMENSION(SIZE(PTHLM))   :: ZTHVMIX_1,ZTHVMIX_2 ! Theta and Thetav  of mixtures
 699
 700
 701 ! Variables for dry part
 702
 703 REAL, DIMENSION(SIZE(PTHLM))   :: ZBUO_INTEG,&         ! Temporary integral Buoyancy
 704                                   ZDZ_HALF,&           ! half-DeltaZ between 2 flux points
 705                                   ZDZ_STOP,&           ! Exact Height of the LCL
 706                                   ZTHV_MINUS_HALF,&    ! Thv at flux point(kk)
 707                                   ZTHV_PLUS_HALF,&     ! Thv at flux point(kk+1)
 708                                   ZCOEFF_MINUS_HALF,&  ! Variation of Thv between mass points kk-1 and kk
 709                                   ZCOEFF_PLUS_HALF,&   ! Variation of Thv between mass points kk and kk+1
 710                                   ZCOTHVU_MINUS_HALF,& ! Variation of Thvup between flux point kk and mass point kk
 711                                   ZCOTHVU_PLUS_HALF,&  ! Variation of Thvup between mass point kk and flux point kk+1
 712                                   ZW2_HALF             ! w**2 at mass point KK
 713
 714 REAL, DIMENSION(SIZE(PTHLM))   :: ZCOPRE_MINUS_HALF,&  ! Variation of pressure between mass points kk-1 and kk
 715                                   ZCOPRE_PLUS_HALF,&   ! Variation of pressure between mass points kk and kk+1
 716                                   ZPRE_MINUS_HALF,&    ! pressure at flux point kk
 717                                   ZPRE_PLUS_HALF,&     ! pressure at flux point kk+1
 718                                   ZTHV_UP_F1,&         ! thv_up at flux point kk
 719                                   ZTHV_UP_F2           ! thv_up at flux point kk+1
 720 REAL, DIMENSION(SIZE(PTHLM))   :: ZCOEFF_QSAT,&        ! variation of Qsat at the transition between dry part and cloudy part
 721                                   ZRC_ORD,&            !
 722                                   ZPART_DRY            ! part of dry part at the transition level
 723 !
 724 REAL, DIMENSION(SIZE(PTHLM))   :: ZQVSAT ! QV at saturation
 725
 726 REAL, DIMENSION(SIZE(PTHLM))   :: PT ! temperature to compute saturation vapour mixing ratio
 727
 728 REAL, DIMENSION(SIZE(PTHVM,1),SIZE(PTHVM,2))   ::ZG_O_THVREF
 729 !
 730 LOGICAL, DIMENSION(SIZE(OTEST,1)) :: GTEST_LOCAL_LCL,& ! true if LCL found between flux point KK and mass point KK
 731                                      GTEST_LOCAL_LCL2  ! true if LCL found between mass point KK and flux point KK+1
 732 !
 733 REAL     :: ZRDORV       ! RD/RV
 734 REAL     :: ZRVORD       ! RV/RD
 735 INTEGER  :: ILON, ITEST, IKB         !
 736
 737
 738 !----------------------------------------------------------------------------------
 739
 740 !                1.3 Initialisation
 741 !                ------------------
 742   IKB=1 ! modif WRF JP
 743
 744
 745   ZRDORV   = XRD / XRV   !=0.622
 746   ZRVORD   = XRV / XRD   !=1.607
 747   ZG_O_THVREF=XG/PTHVM
 748
 749   ZCOEFF_QSAT=0.
 750   ZRC_ORD=0.
 751   ZPART_DRY=1.
 752   GTEST_LOCAL_LCL=.FALSE.
 753   ZDZ_HALF(:) = (PZZ(:,KK+1)-PZZ(:,KK))/2.
 754   ZDZ_STOP(:) = ZDZ_HALF(:)
 755
 756   ZKIC(:)=0.1  ! starting value for critical mixed fraction for CLoudy Part
 757
 758
 759 !                Computation of KIC
 760 !                ---------------------
 761
 762 !        2.1    Compute critical mixed fraction by estimating unknown
 763 !               T^mix r_c^mix and r_i^mix from thl^mix and r_t^mix
 764 !               We determine the zero crossing of the linear curve
 765 !               evaluating the derivative using ZMIXF=0.1.
 766 !               -----------------------------------------------------
 767
 768   ZMIXTHL(:) = ZKIC(:) * PTHLM(:)+(1. - ZKIC(:))*PTHL_UP(:)
 769   ZMIXRT(:)  = ZKIC(:) * PRTM(:)+(1. - ZKIC(:))*PRT_UP(:)
 770
 771 ! MIXTURE FOR CLOUDY PART
 772   !  Compute pressure at flux level KK and at flux Level KK+1
 773
 774
 775 IF (KK==IKB) THEN !MODIF WRF JP
 776   ZCOPRE_MINUS_HALF(:) = 0. !MODIF WRF JP
 777 ELSE!MODIF WRF JP
 778   ZCOPRE_MINUS_HALF(:) = ((PPABSM(:,KK)-PPABSM(:,KK-1))/PDZZ(:,KK))
 779 ENDIF!MODIF WRF JP
 780   ZCOPRE_PLUS_HALF(:) = ((PPABSM(:,KK+1)-PPABSM(:,KK))/PDZZ(:,KK+1))
 781
 782 IF (KK==IKB) THEN !MODIF WRF JP
 783   ZPRE_MINUS_HALF(:)= PPABSM(:,KK)
 784 ELSE!MODIF WRF JP
 785   ZPRE_MINUS_HALF(:)= ZCOPRE_MINUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-1))+PPABSM(:,KK-1)
 786 ENDIF!MODIF WRF JP
 787   ZPRE_PLUS_HALF(:) = ZCOPRE_PLUS_HALF*0.5*(PZZ(:,KK+1)-PZZ(:,KK))+PPABSM(:,KK)
 788
 789
 790   !  Compute non cons. var. of mixture at the mass level
 791   CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
 792                PPABSM(:,KK),ZMIXTHL,ZMIXRT,&
 793                ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX)
 794
 795
 796   ! Compute theta_v of mixture at mass level KK for KF90
 797   ZTHVMIX_1(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
 798
 799   !  Compute non cons. var. of mixture at the flux level KK+1
 800   CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
 801                ZPRE_PLUS_HALF,ZMIXTHL,ZMIXRT,&
 802                ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX)
 803
 804
 805   ! compute theta_v of mixture at the flux level KK+1 for KF90
 806   ZTHVMIX_2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
 807
 808
 809 !        2.1    Compute critical mixed fraction by estimating unknown
 810 !               T^mix r_c^mix and r_i^mix from thl^mix and r_t^mix
 811 !               We determine the zero crossing of the linear curve
 812 !               evaluating the derivative using ZMIXF=0.1.
 813 !               -----------------------------------------------------
 814
 815
 816 ! THV_UP FOR DRY PART
 817   ! Compute theta_v of updraft at flux level KK
 818   CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
 819                ZPRE_MINUS_HALF,PTHL_UP,PRT_UP,&
 820                ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX)
 821   ZTHV_UP_F1(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
 822
 823   ! Compute theta_v of updraft at mass level KK
 824   CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
 825                PPABSM(:,KK),PTHL_UP,PRT_UP,&
 826                ZTHMIX,ZRVMIX,ZRCMIX,ZRIMIX)
 827   ZTHV_UP(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
 828
 829   ! Compute theta_v of updraft at flux level KK+1
 830   CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,&
 831                ZPRE_PLUS_HALF,PTHL_UP,PRT_UP,&
 832                ZTHMIX_F2,ZRVMIX_F2,ZRCMIX_F2,ZRIMIX_F2)
 833   ZTHV_UP_F2(:) = ZTHMIX_F2(:)*(1.+ZRVORD*ZRVMIX_F2(:))/(1.+PRT_UP(:))
 834
 835
 836 !
 837 !*   2.2     Compute final values for entr. and detr.
 838 !            ----------------------------------------
 839 !
 840 ! Dry PART
 841
 842  ! Computation of integral entrainment and detrainment between flux level KK
 843  ! and mass level KK
 844
 845   WHERE ((ZRCMIX(:)>0.).AND.(.NOT.OTESTLCL))
 846 ! If rc is found between flux level KK and mass level KK
 847 ! a part of dry entrainment/detrainment is defined
 848 ! the exact height of LCL is also determined
 849      ZCOEFF_QSAT(:) = (ZRCMIX_F2(:) - ZRCMIX(:))/ ZDZ_HALF(:)
 850      ZRC_ORD(:) = ZRCMIX(:) - ZCOEFF_QSAT(:) * ZDZ_HALF(:)
 851      ZDZ_STOP = (- ZRC_ORD(:)/ZCOEFF_QSAT(:))
 852      ZPART_DRY(:) = ZDZ_STOP / (PZZ(:,KK+1)-PZZ(:,KK))
 853      GTEST_LOCAL_LCL(:)=.TRUE.
 854
 855   ENDWHERE
 856
 857 IF (KK==IKB) THEN !MODIF WRF JP
 858   ZCOEFF_MINUS_HALF = 0.
 859 ELSE!MODIF WRF JP
 860   ZCOEFF_MINUS_HALF = ((PTHVM(:,KK)-PTHVM(:,KK-1))/PDZZ(:,KK))
 861 ENDIF!MODIF WRF JP
 862
 863   ZCOEFF_PLUS_HALF  = ((PTHVM(:,KK+1)-PTHVM(:,KK))/PDZZ(:,KK+1))
 864
 865   ZCOTHVU_MINUS_HALF = (ZTHV_UP(:)-ZTHV_UP_F1(:))/ZDZ_HALF(:)
 866   ZCOTHVU_PLUS_HALF  = (ZTHV_UP_F2(:)-ZTHV_UP(:))/ZDZ_HALF(:)
 867
 868 IF (KK==IKB) THEN !MODIF WRF JP
 869   ZTHV_MINUS_HALF = PTHVM(:,KK)
 870 ELSE!MODIF WRF JP
 871   ZTHV_MINUS_HALF = ZCOEFF_MINUS_HALF*0.5*(PZZ(:,KK)-PZZ(:,KK-1))+PTHVM(:,KK-1)
 872 ENDIF!MODIF WRF JP
 873
 874   ZTHV_PLUS_HALF  = ZCOEFF_PLUS_HALF*0.5*(PZZ(:,KK+1)-PZZ(:,KK))+ ZTHV_MINUS_HALF  ! BUG JP 16082010
 875
 876   ! Integral Buoyancy between flux level KK and mass level KK
 877
 878   PBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_HALF(:)*&
 879               (0.5*( ZCOTHVU_MINUS_HALF - ZCOEFF_MINUS_HALF)*ZDZ_HALF(:) &
 880                 - ZTHV_MINUS_HALF + ZTHV_UP_F1(:) )
 881
 882
 883
 884   WHERE ((OTEST).AND.(.NOT.OTESTLCL))
 885      PENTR=0.
 886      PDETR=0.
 887
 888
 889      ZBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_STOP(:)*&
 890                 (0.5 * (  - ZCOEFF_MINUS_HALF)* ZDZ_STOP(:) &
 891                   - ZTHV_MINUS_HALF + ZTHV_UP_F1(:) )
 892
 893
 894      WHERE (ZBUO_INTEG(:)>=0.)
 895          PENTR = 0.5/(XABUO-XBENTR*XENTR_DRY)*&
 896                  LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/PW_UP2(:,KK))* &
 897                  ZBUO_INTEG)
 898          PDETR = 0.
 899
 900          ZW2_HALF = PW_UP2(:,KK) +  2*(XABUO-XBENTR*XENTR_DRY)*(ZBUO_INTEG)
 901      ELSEWHERE
 902          PENTR = 0.
 903          PDETR = 0.5/(XABUO)*&
 904                  LOG(1.+ (2.*(XABUO)/PW_UP2(:,KK))* &
 905                  MAX(0.,-ZBUO_INTEG))
 906
 907          ZW2_HALF = PW_UP2(:,KK) +  2*(XABUO)*(ZBUO_INTEG)
 908      ENDWHERE
 909
 910  ENDWHERE
 911
 912
 913  ZDZ_STOP(:) = ZDZ_HALF(:)
 914
 915 ! total Integral Buoyancy between flux level KK and flux level KK+1
 916
 917  PBUO_INTEG = PBUO_INTEG + ZG_O_THVREF(:,KK)*ZDZ_HALF(:)*&
 918                 (0.5*(ZCOTHVU_PLUS_HALF - ZCOEFF_PLUS_HALF)* ZDZ_HALF(:) - &
 919                 PTHVM(:,KK) + ZTHV_UP(:) )
 920
 921
 922  WHERE ((((ZRCMIX_F2(:)>0.).AND.(ZRCMIX(:)<=0.)).AND.(.NOT.OTESTLCL)).AND.(.NOT.GTEST_LOCAL_LCL(:)))
 923  ! If rc is found between mass level KK and flux level KK+1
 924  ! a part of dry entrainment is defined
 925  ! the exact height of LCL is also determined
 926      PT(:)=ZTHMIX_F2(:)*((PPABSM(:,KK+1)/XP00)**(XRD/XCPD))
 927      ZQVSAT(:)=EXP( XALPW - XBETAW/PT(:) - XGAMW*LOG(PT(:))  )
 928      ZQVSAT(:)=XRD/XRV*ZQVSAT(:)/PPABSM(:,KK+1)   &
 929                    / (1.+(XRD/XRV-1.)*ZQVSAT(:)/PPABSM(:,KK+1))
 930      ZCOEFF_QSAT(:) = (PRT_UP(:) - ZQVSAT(:) - &
 931                 ZRCMIX(:))/ (0.5* (PZZ(:,KK+2)-PZZ(:,KK+1)))
 932      ZRC_ORD(:) = ZRCMIX_F2(:) - ZCOEFF_QSAT(:) * ZDZ_HALF(:)
 933      ZDZ_STOP = (- ZRC_ORD(:)/ZCOEFF_QSAT(:))
 934      ZPART_DRY(:) = 0.5+ZDZ_STOP  / (PZZ(:,KK+1)-PZZ(:,KK))
 935      GTEST_LOCAL_LCL2(:)=.TRUE.
 936  ENDWHERE
 937
 938
 939  WHERE (((OTEST).AND.(.NOT.OTESTLCL)).AND.(.NOT.GTEST_LOCAL_LCL(:)))
 940
 941      ZBUO_INTEG = ZG_O_THVREF(:,KK)*ZDZ_STOP(:)*&
 942                 (0.5*( - ZCOEFF_PLUS_HALF)* ZDZ_STOP(:)&
 943                 - PTHVM(:,KK) + ZTHV_UP(:) )
 944
 945
 946
 947      WHERE (ZW2_HALF>0.)
 948         WHERE (ZBUO_INTEG(:)>=0.)
 949            PENTR = PENTR + 0.5/(XABUO-XBENTR*XENTR_DRY)* &
 950                 LOG(1.+ (2.*(XABUO-XBENTR*XENTR_DRY)/ZW2_HALF(:)) * ZBUO_INTEG)
 951
 952            PDETR = PDETR
 953
 954         ELSEWHERE
 955           PENTR = PENTR
 956
 957           PDETR = PDETR + 0.5/(XABUO)* &
 958                 LOG(1.+ (2.*(XABUO)/ZW2_HALF(:)) * &
 959                 MAX(-ZBUO_INTEG,0.))
 960         ENDWHERE
 961      ELSEWHERE
 962         ! if w**2<0 the updraft is stopped
 963            OTEST=.FALSE.
 964            PENTR = PENTR
 965            PDETR = PDETR
 966      ENDWHERE
 967
 968  ENDWHERE
 969  PENTR = XENTR_DRY*PENTR/(PZZ(:,KK+1)-PZZ(:,KK))
 970  PDETR = XDETR_DRY*PDETR/(PZZ(:,KK+1)-PZZ(:,KK))
 971  PDETR = MAX(ZPART_DRY(:)*XDETR_LUP/(PLUP-0.5*(PZZ(:,KK)+PZZ(:,KK+1))),PDETR)
 972
 973
 974
 975
 976 ! compute final value of critical mixed fraction using theta_v
 977 ! of mixture, grid-scale and updraft in cloud
 978
 979  WHERE ((OTEST).AND.(OTESTLCL))
 980      ZKIC(:) = MAX(0.,ZTHV_UP(:)-PTHVM(:,KK))*ZKIC(:) /  &
 981                  (ZTHV_UP(:)-ZTHVMIX_1(:)+1.E-10)
 982
 983      ZKIC(:) = MAX(0., MIN(1., ZKIC(:)))
 984
 985      ZEPSI(:) = ZKIC(:) **2.
 986      ZDELTA(:) = (1.-ZKIC(:))**2.
 987      ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
 988      ZCOEFFMF_CLOUD(:)=XENTR_MF * XG / XCRAD_MF
 989      PENTR(:) = ZCOEFFMF_CLOUD(:)*ZEPSI_CLOUD(:)
 990      PDETR(:) = ZCOEFFMF_CLOUD(:)*ZDELTA(:)
 991  ENDWHERE
 992
 993 ! compute final value of critical mixed fraction using theta_v
 994 ! of mixture, grid-scale and updraft in cloud
 995
 996
 997  WHERE (((OTEST).AND.(.NOT.(OTESTLCL))).AND.((GTEST_LOCAL_LCL(:).OR.GTEST_LOCAL_LCL2(:))))
 998      ZKIC(:) = MAX(0.,ZTHV_UP_F2(:)-ZTHV_PLUS_HALF)*ZKIC(:) /  &
 999                        (ZTHV_UP_F2(:)-ZTHVMIX_2(:)+1.E-10)
1000      ZKIC(:) = MAX(0., MIN(1., ZKIC(:)))
1001      ZEPSI(:) = ZKIC(:) **2.
1002      ZDELTA(:) = (1.-ZKIC(:))**2.
1003      ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
1004      ZCOEFFMF_CLOUD(:)=XENTR_MF * XG / XCRAD_MF
1005      PENTR(:) = PENTR+(1.-ZPART_DRY(:))*ZCOEFFMF_CLOUD(:)*ZEPSI_CLOUD(:)
1006      PDETR(:) = PDETR+(1.-ZPART_DRY(:))*ZCOEFFMF_CLOUD(:)*ZDELTA(:)
1007  ENDWHERE
1008
1009
1010
1011 END SUBROUTINE COMPUTE_ENTR_DETR
1012
1013
1014 !     #################################################################
1015       SUBROUTINE COMPUTE_FUNCTION_THERMO_MF( KRR,KRRL,KRRI,                  &
1016                                        PTH, PR, PEXN, PPABS,                 &
1017                                        PT,PLVOCPEXN,PLSOCPEXN,               &
1018                                        PAMOIST,PATHETA                       )
1019 !     #################################################################
1020 !
1021 !!
1022 !!****  *COMPUTE_FUNCTION_THERMO_MF* -
1023 !!
1024 !!
1025 !! --------------------------------------------------------------------------
1026 !
1027
1028 !
1029 IMPLICIT NONE
1030 !
1031 !
1032 !*      0.1  declarations of arguments
1033 !
1034 INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
1035 INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
1036 INTEGER,                INTENT(IN)   :: KRRI          ! number of ice water var.
1037
1038 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PTH      ! theta
1039 REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PR       ! water species
1040 REAL, DIMENSION(:,:,:)  , INTENT(IN) :: PPABS,PEXN    ! pressure, Exner funct.
1041
1042 REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PT      ! temperature
1043 REAL, DIMENSION(:,:,:), INTENT(OUT)   :: PLVOCPEXN,PLSOCPEXN     ! L/(cp*Pi)
1044
1045 REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PAMOIST,PATHETA
1046 !
1047 !-------------------------------------------------------------------------------
1048 !
1049 !*       0.2   Declarations of local variables
1050 !
1051 REAL                :: ZEPS         ! XMV / XMD
1052 REAL, DIMENSION(SIZE(PTH,1),SIZE(PTH,2),SIZE(PTH,3)) ::     &
1053           ZCP,                        &  ! Cp
1054           ZE,                         &  ! Saturation mixing ratio
1055           ZDEDT,                      &  ! Saturation mixing ratio derivative
1056           ZFRAC_ICE,                  &  ! Ice fraction
1057           ZAMOIST_W,                  &  ! Coefficients for s = f (Thetal,Rnp)
1058           ZATHETA_W,                  &  !
1059           ZAMOIST_I,                  &  !
1060           ZATHETA_I                      !
1061
1062 INTEGER             :: JRR
1063 !
1064 !-------------------------------------------------------------------------------
1065 !
1066   ZEPS = XMV / XMD
1067
1068   PLVOCPEXN(:,:,:) = 0.
1069   PLSOCPEXN(:,:,:) = 0.
1070   PAMOIST(:,:,:) = 0.
1071   PATHETA(:,:,:) = 0.
1072   ZFRAC_ICE(:,:,:) = 0.0
1073   ZAMOIST_W(:,:,:) = 0.0
1074   ZATHETA_W(:,:,:) = 0.0
1075   ZAMOIST_I(:,:,:) = 0.0
1076   ZATHETA_I(:,:,:) = 0.0
1077 !
1078 !*       Cph
1079 !
1080 ZCP=XCPD
1081
1082 IF (KRR > 0) ZCP(:,:,:) = ZCP(:,:,:) + XCPV * PR(:,:,:,1)
1083
1084 DO JRR = 2,1+KRRL  ! loop on the liquid components
1085    ZCP(:,:,:)  = ZCP(:,:,:) + XCL * PR(:,:,:,JRR)
1086 END DO
1087
1088 DO JRR = 2+KRRL,1+KRRL+KRRI ! loop on the solid components
1089   ZCP(:,:,:)  = ZCP(:,:,:)  + XCI * PR(:,:,:,JRR)
1090 END DO
1091
1092 !*      Temperature
1093 !
1094 PT(:,:,:) =  PTH(:,:,:) * PEXN(:,:,:)
1095 !
1096 !
1097 !! Liquid water
1098 !
1099 IF ( KRRL >= 1 ) THEN
1100 !
1101 !*       Lv/Cph
1102 !
1103   PLVOCPEXN(:,:,:) = (XLVTT + (XCPV-XCL) *  (PT(:,:,:)-XTT) ) / ZCP(:,:,:)
1104 !
1105 !*      Saturation vapor pressure with respect to water
1106 !
1107   ZE(:,:,:) =  EXP( XALPW - XBETAW/PT(:,:,:) - XGAMW*ALOG( PT(:,:,:) ) )
1108 !
1109 !*      Saturation  mixing ratio with respect to water
1110 !
1111   ZE(:,:,:) =  ZE(:,:,:) * ZEPS / ( PPABS(:,:,:) - ZE(:,:,:) )
1112 !
1113 !*      Compute the saturation mixing ratio derivative (rvs')
1114 !
1115   ZDEDT(:,:,:) = ( XBETAW / PT(:,:,:)  - XGAMW ) / PT(:,:,:)   &
1116                  * ZE(:,:,:) * ( 1. + ZE(:,:,:) / ZEPS )
1117 !
1118 !*      Compute Amoist
1119 !
1120   ZAMOIST_W(:,:,:)=  0.5 / ( 1.0 + ZDEDT(:,:,:) * PLVOCPEXN(:,:,:) )
1121 !
1122 !*      Compute Atheta
1123 !
1124   ZATHETA_W(:,:,:)= ZAMOIST_W(:,:,:) * PEXN(:,:,:) *                       &
1125         ( ( ZE(:,:,:) - PR(:,:,:,1) ) * PLVOCPEXN(:,:,:) /                  &
1126           ( 1. + ZDEDT(:,:,:) * PLVOCPEXN(:,:,:) )           *              &
1127           (                                                                 &
1128            ZE(:,:,:) * (1. + ZE(:,:,:)/ZEPS)                                &
1129                         * ( -2.*XBETAW/PT(:,:,:) + XGAMW ) / PT(:,:,:)**2   &
1130           +ZDEDT(:,:,:) * (1. + 2. * ZE(:,:,:)/ZEPS)                        &
1131                         * ( XBETAW/PT(:,:,:) - XGAMW ) / PT(:,:,:)          &
1132           )                                                                 &
1133          - ZDEDT(:,:,:)                                                     &
1134         )
1135
1136 !
1137 !! Solid water
1138 !
1139   IF ( KRRI >= 1 ) THEN
1140
1141 !
1142 !*       Fraction of ice
1143 !
1144     WHERE(PR(:,:,:,2)+PR(:,:,:,4)>0.0)
1145       ZFRAC_ICE(:,:,:) = PR(:,:,:,4) / ( PR(:,:,:,2)+PR(:,:,:,4) )
1146     END WHERE
1147 !
1148 !*       Ls/Cph
1149 !
1150     PLSOCPEXN(:,:,:) = (XLSTT + (XCPV-XCI) *  (PT(:,:,:)-XTT) ) / ZCP(:,:,:)
1151 !
1152 !*      Saturation vapor pressure with respect to ice
1153 !
1154     ZE(:,:,:) =  EXP( XALPI - XBETAI/PT(:,:,:) - XGAMI*ALOG( PT(:,:,:) ) )
1155 !
1156 !*      Saturation  mixing ratio with respect to ice
1157 !
1158     ZE(:,:,:) =  ZE(:,:,:) * ZEPS / ( PPABS(:,:,:) - ZE(:,:,:) )
1159 !
1160 !*      Compute the saturation mixing ratio derivative (rvs')
1161 !
1162     ZDEDT(:,:,:) = ( XBETAI / PT(:,:,:)  - XGAMI ) / PT(:,:,:)   &
1163                    * ZE(:,:,:) * ( 1. + ZE(:,:,:) / ZEPS )
1164 !
1165 !*      Compute Amoist
1166 !
1167     ZAMOIST_I(:,:,:)=  0.5 / ( 1.0 + ZDEDT(:,:,:) * PLSOCPEXN(:,:,:) )
1168 !
1169 !*      Compute Atheta
1170 !
1171     ZATHETA_I(:,:,:)= ZAMOIST_I(:,:,:) * PEXN(:,:,:) *                     &
1172         ( ( ZE(:,:,:) - PR(:,:,:,1) ) * PLSOCPEXN(:,:,:) /                  &
1173           ( 1. + ZDEDT(:,:,:) * PLSOCPEXN(:,:,:) )           *              &
1174           (                                                                 &
1175            ZE(:,:,:) * (1. + ZE(:,:,:)/ZEPS)                                &
1176                         * ( -2.*XBETAI/PT(:,:,:) + XGAMI ) / PT(:,:,:)**2   &
1177           +ZDEDT(:,:,:) * (1. + 2. * ZE(:,:,:)/ZEPS)                        &
1178                         * ( XBETAI/PT(:,:,:) - XGAMI ) / PT(:,:,:)          &
1179           )                                                                 &
1180          - ZDEDT(:,:,:)                                                     &
1181         )
1182
1183
1184   ENDIF
1185
1186   PAMOIST(:,:,:) = (1.0-ZFRAC_ICE(:,:,:))*ZAMOIST_W(:,:,:) &
1187                          +ZFRAC_ICE(:,:,:) *ZAMOIST_I(:,:,:)
1188   PATHETA(:,:,:) = (1.0-ZFRAC_ICE(:,:,:))*ZATHETA_W(:,:,:) &
1189                          +ZFRAC_ICE(:,:,:) *ZATHETA_I(:,:,:)
1190
1191 !
1192 !*      Lv/Cph/Exner and  Ls/Cph/Exner
1193 !
1194   PLVOCPEXN(:,:,:) = PLVOCPEXN(:,:,:) / PEXN(:,:,:)
1195   PLSOCPEXN(:,:,:) = PLSOCPEXN(:,:,:) / PEXN(:,:,:)
1196 !
1197 ENDIF
1198
1199 END SUBROUTINE COMPUTE_FUNCTION_THERMO_MF
1200 !
1201 !     #################################################################
1202       SUBROUTINE COMPUTE_UPDRAFT(OMIXUV,PZZ,PDZZ,KK,              &
1203                                  PSFTH,PSFRV,                     &
1204                                  PPABSM,PRHODREF,PUM,PVM, PTKEM,  &
1205                                  PTHM,PRVM,PRCM,PRIM,PTHLM,PRTM,  &
1206                                  PTHL_UP,PRT_UP,                  &
1207                                  PRV_UP,PRC_UP,PRI_UP,PTHV_UP,    &
1208                                  PW_UP,PU_UP, PV_UP,              &
1209                                  PFRAC_UP,PEMF,PDETR,PENTR,       &
1210                                  KKLCL,KKETL,KKCTL)
1211 !     #################################################################
1212 !!
1213 !!****  *COMPUTE_UPDRAFT* - calculates caracteristics of the updraft
1214 !!
1215 !! --------------------------------------------------------------------------
1216
1217
1218
1219 IMPLICIT NONE
1220
1221 !*                    1.1  Declaration of Arguments
1222 !
1223 !
1224 !
1225 LOGICAL,                INTENT(IN) :: OMIXUV    ! True if mixing of momentum
1226 REAL, DIMENSION(:,:), INTENT(IN)   :: PZZ       !  Height at the flux point
1227 REAL, DIMENSION(:,:), INTENT(IN)   :: PDZZ      !  Metrics coefficient
1228
1229 INTEGER,              INTENT(IN)   ::  KK
1230 REAL, DIMENSION(:),   INTENT(IN)   ::  PSFTH,PSFRV
1231 ! normal surface fluxes of theta,rv,(u,v) parallel to the orography
1232 !
1233 REAL, DIMENSION(:,:),   INTENT(IN) ::  PPABSM     ! Pressure at t-dt
1234 REAL, DIMENSION(:,:),   INTENT(IN) ::  PRHODREF   ! dry density of the
1235                                                   ! reference state
1236 REAL, DIMENSION(:,:),   INTENT(IN) ::  PUM        ! u mean wind
1237 REAL, DIMENSION(:,:),   INTENT(IN) ::  PVM        ! v mean wind
1238 REAL, DIMENSION(:,:),   INTENT(IN) ::  PTKEM      ! TKE at t-dt
1239 !
1240 REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTHM           ! liquid pot. temp. at t-dt
1241 REAL, DIMENSION(:,:),   INTENT(IN)   ::  PRVM,PRCM,PRIM ! vapor mixing ratio at t-dt
1242 REAL, DIMENSION(:,:),   INTENT(IN)   ::  PTHLM,PRTM     ! cons. var. at t-dt
1243
1244
1245 REAL, DIMENSION(:,:),   INTENT(OUT)  ::  PTHL_UP,PRT_UP   ! updraft properties
1246 REAL, DIMENSION(:,:),   INTENT(OUT)  ::  PU_UP, PV_UP     ! updraft wind components
1247 REAL, DIMENSION(:,:),   INTENT(OUT)  ::  PRV_UP,PRC_UP, & ! updraft rv, rc
1248                                          PRI_UP,PTHV_UP,& ! updraft ri, THv
1249                                          PW_UP,PFRAC_UP   ! updraft w, fraction
1250
1251
1252 REAL, DIMENSION(:,:),   INTENT(OUT)  ::  PEMF,PDETR,PENTR ! Mass_flux,
1253                                                           ! detrainment,entrainment
1254 INTEGER, DIMENSION(:),  INTENT(OUT)  ::  KKLCL,KKETL,KKCTL! LCL, ETL, CTL
1255 !                       1.2  Declaration of local variables
1256 !
1257 !
1258 ! Mean environment variables at t-dt at flux point
1259 REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) ::    &
1260                         ZTHM_F,ZRVM_F,ZRCM_F,ZRIM_F   ! Theta,rv,rl,ri of
1261                                                       ! updraft environnement
1262 REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) ::    &
1263                         ZRTM_F, ZTHLM_F, ZTKEM_F,&    ! rt, thetal,TKE,pressure,
1264                         ZUM_F,ZVM_F,ZRHO_F,      &    ! density,momentum
1265                         ZPRES_F,ZTHVM_F,ZTHVM,   &    ! interpolated at the flux point
1266                         ZG_O_THVREF,             &    ! g*ThetaV ref
1267                         ZW_UP2                        ! w**2  of the updraft
1268
1269
1270 REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) ::  &
1271                         ZTH_UP,                  &    ! updraft THETA
1272                         ZFRAC_ICE                     ! Ice fraction
1273
1274 REAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) ::  ZCOEF  ! diminution coefficient for too high clouds
1275
1276 REAL, DIMENSION(SIZE(PSFTH,1) )            ::  ZWTHVSURF  ! Surface w'thetav'
1277 CHARACTER(LEN=1)                           ::  YFRAC_ICE  ! Ice Fraction
1278                                                           ! precribed or computed
1279
1280 REAL  :: ZRDORV       ! RD/RV
1281 REAL  :: ZRVORD       ! RV/RD
1282
1283
1284 REAL, DIMENSION(SIZE(PTHM,1)) :: ZMIX1,ZMIX2,ZMIX3
1285 REAL, DIMENSION(SIZE(PTHM,1)) :: ZBUO_INTEG           ! Integrated Buoyancy
1286
1287 REAL, DIMENSION(SIZE(PTHM,1)) :: ZLUP         ! Upward Mixing length from the ground
1288 REAL, DIMENSION(SIZE(PTHM,1)) :: ZDEPTH       ! Deepness limit for cloud
1289
1290 INTEGER  :: IKB,IKE            ! index value for the Beginning and the End
1291                                ! of the physical domain for the mass points
1292 INTEGER  :: IKU           ! array size in k
1293 INTEGER  :: JK,JI,JSV          ! loop counters
1294
1295 LOGICAL, DIMENSION(SIZE(PTHM,1)) ::  GTEST,GTESTLCL,GTESTETL
1296                                ! Test if the ascent continue, if LCL or ETL is reached
1297 LOGICAL                          ::  GLMIX
1298                                ! To choose upward or downward mixing length
1299 LOGICAL, DIMENSION(SIZE(PTHM,1))              :: GWORK1
1300 LOGICAL, DIMENSION(SIZE(PTHM,1),SIZE(PTHM,2)) :: GWORK2
1301
1302 INTEGER  :: ITEST
1303 REAL  :: ZTMAX,ZRMAX  ! control value
1304 ! Thresholds for the  perturbation of
1305 ! theta_l and r_t at the first level of the updraft
1306 ZTMAX=2.0
1307 ZRMAX=1.E-3
1308 !------------------------------------------------------------------------
1309
1310 !                     INITIALISATION
1311
1312 !                 Local variables, internal domain
1313 ! Internal Domain
1314 IKU=SIZE(PTHM,2)
1315 IKB=1        ! modif WRF JP
1316 IKE=IKU-1      ! modif WRF JP
1317
1318 ! Initialisation of intersesting Level :LCL,ETL,CTL
1319 KKLCL(:)=IKE
1320 KKETL(:)=IKE
1321 KKCTL(:)=IKE
1322
1323 !
1324 ! Initialisation
1325 PRV_UP(:,:)=0.
1326 PRC_UP(:,:)=0.
1327
1328 PW_UP(:,:)=0.
1329 PEMF(:,:)=0.
1330 PDETR(:,:)=0.
1331 PENTR(:,:)=0.
1332 ZTH_UP(:,:)=0.
1333 PFRAC_UP(:,:)=0.
1334 PTHV_UP(:,:)=0.
1335
1336 !no ice cloud coded yet
1337 PRI_UP(:,:)=0.
1338 ZFRAC_ICE(:,:)=0.
1339 YFRAC_ICE='T'
1340
1341 ZBUO_INTEG=0.
1342 ! Initialisation of the constants
1343 ZRDORV   = XRD / XRV   !=0.622
1344 ZRVORD   = (XRV / XRD)
1345
1346
1347 ! Initialisation of environment variables at t-dt
1348
1349 ! variables at flux level
1350
1351
1352 ZTHM_F (:,IKB+1:IKU) = 0.5*(PTHM(:,IKB:IKU-1)+PTHM(:,IKB+1:IKU))
1353 ZTHLM_F(:,IKB+1:IKU) = 0.5*(PTHLM(:,IKB:IKU-1)+PTHLM(:,IKB+1:IKU))
1354 ZRTM_F (:,IKB+1:IKU) = 0.5*(PRTM(:,IKB:IKU-1)+PRTM(:,IKB+1:IKU))
1355 ZTKEM_F(:,IKB+1:IKU) = 0.5*(PTKEM(:,IKB:IKU-1)+PTKEM(:,IKB+1:IKU))
1356 ZPRES_F(:,IKB+1:IKU) = 0.5*(PPABSM(:,IKB:IKU-1)+PPABSM(:,IKB+1:IKU))
1357 ZRHO_F (:,IKB+1:IKU) = 0.5*(PRHODREF(:,IKB:IKU-1)+PRHODREF(:,IKB+1:IKU))
1358 ZRVM_F (:,IKB+1:IKU) = 0.5*(PRVM(:,IKB:IKU-1)+PRVM(:,IKB+1:IKU))
1359 ZUM_F  (:,IKB+1:IKU) = 0.5*(PUM(:,IKB:IKU-1)+PUM(:,IKB+1:IKU))
1360 ZVM_F  (:,IKB+1:IKU) = 0.5*(PVM(:,IKB:IKU-1)+PVM(:,IKB+1:IKU))
1361
1362
1363
1364 ZTHM_F (:,IKB) = PTHM(:,IKB)
1365 ZTHLM_F(:,IKB) = PTHLM(:,IKB)
1366 ZRTM_F (:,IKB) = PRTM (:,IKB)
1367 ZTKEM_F(:,IKB) = PTKEM(:,IKB)
1368 ZPRES_F(:,IKB) = PPABSM(:,IKB)
1369 ZRHO_F(:,IKB)  = PRHODREF(:,IKB)
1370 ZRVM_F(:,IKB)  = PRVM(:,IKB)
1371 ZUM_F(:,IKB)   = PUM(:,IKB)
1372 ZVM_F(:,IKB)   = PVM(:,IKB)
1373
1374
1375 ! thetav at mass and flux levels
1376 ZTHVM_F(:,:)=ZTHM_F(:,:)*((1.+ZRVORD*ZRVM_F(:,:))/(1.+ZRTM_F(:,:)))
1377 ZTHVM(:,:)=PTHM(:,:)*((1.+ZRVORD*PRVM(:,:))/(1.+PRTM(:,:)))
1378
1379
1380 !
1381 !          Initialisation of updraft characteristics
1382 PTHL_UP(:,:)=ZTHLM_F(:,:)
1383 PRT_UP(:,:)=ZRTM_F(:,:)
1384 ZW_UP2(:,:)=0.
1385 PTHV_UP(:,:)=ZTHVM_F(:,:)
1386 PU_UP(:,:)=ZUM_F(:,:)
1387 PV_UP(:,:)=ZVM_F(:,:)
1388
1389
1390 ! Computation or initialisation of updraft characteristics at the KK level
1391 ! thetal_up,rt_up,thetaV_up, w-squared,Buoyancy term and mass flux (PEMF)
1392
1393 ! 03/2009
1394 !PTHL_UP(:,KK)= ZTHLM_F(:,KK)+(PSFTH(:)/SQRT(ZTKEM_F(:,KK)))*XALP_PERT
1395 !PRT_UP(:,KK) = ZRTM_F(:,KK)+(PSFRV(:)/SQRT(ZTKEM_F(:,KK)))*XALP_PERT
1396 PTHL_UP(:,KK)= ZTHLM_F(:,KK)+MAX(0.,MIN(ZTMAX,(PSFTH(:)/SQRT(ZTKEM_F(:,KK)))*XALP_PERT))
1397 PRT_UP(:,KK) = ZRTM_F(:,KK)+MAX(0.,MIN(ZRMAX,(PSFRV(:)/SQRT(ZTKEM_F(:,KK)))*XALP_PERT))
1398
1399
1400
1401 ZW_UP2(:,KK) = MAX(0.0001,(2./3.)*ZTKEM_F(:,KK))
1402
1403
1404 ! Computation of non conservative variable for the KK level of the updraft
1405 ! (all or nothing ajustement)
1406 CALL TH_R_FROM_THL_RT_2D(YFRAC_ICE,ZFRAC_ICE(:,KK:KK),ZPRES_F(:,KK:KK), &
1407              PTHL_UP(:,KK:KK),PRT_UP(:,KK:KK),ZTH_UP(:,KK:KK), &
1408              PRV_UP(:,KK:KK),PRC_UP(:,KK:KK),PRI_UP(:,KK:KK))
1409
1410 ! compute updraft thevav and buoyancy term at KK level
1411 PTHV_UP(:,KK) = ZTH_UP(:,KK)*((1+ZRVORD*PRV_UP(:,KK))/(1+PRT_UP(:,KK)))
1412
1413 ! Closure assumption for mass flux at KK level
1414 ZG_O_THVREF=XG/ZTHVM_F
1415
1416 ! compute L_up
1417 GLMIX=.TRUE.
1418 ZTKEM_F(:,KK)=0.
1419
1420 CALL COMPUTE_BL89_ML(PDZZ,ZTKEM_F,ZG_O_THVREF,ZTHVM_F,KK,GLMIX,ZLUP)
1421 ZLUP(:)=MAX(ZLUP(:),1.E-10)
1422
1423 ! Compute Buoyancy flux at the ground
1424 ZWTHVSURF(:) = (ZTHVM_F(:,IKB)/ZTHM_F(:,IKB))*PSFTH(:)+     &
1425                 (0.61*ZTHM_F(:,IKB))*PSFRV(:)
1426
1427 ! Mass flux at KK level (updraft triggered if PSFTH>0.)
1428 WHERE (ZWTHVSURF(:)>0.)
1429   PEMF(:,KK) = XCMF * ZRHO_F(:,KK) * ((ZG_O_THVREF(:,KK))*ZWTHVSURF*ZLUP)**(1./3.)
1430   PFRAC_UP(:,KK)=MIN(PEMF(:,KK)/(SQRT(ZW_UP2(:,KK))*ZRHO_F(:,KK)),XFRAC_UP_MAX)
1431   ZW_UP2(:,KK)=(PEMF(:,KK)/(PFRAC_UP(:,KK)*ZRHO_F(:,KK)))**2
1432   GTEST(:)=.TRUE.
1433 ELSEWHERE
1434   PEMF(:,KK) =0.
1435   GTEST(:)=.FALSE.
1436 ENDWHERE
1437
1438
1439 !--------------------------------------------------------------------------
1440
1441 !                        3. Iteration
1442 !                           ---------
1443 !
1444 ! If GTEST = T the updraft starts from the KK level and stops when GTEST becomes F
1445 !
1446 !
1447 GTESTLCL(:)=.FALSE.
1448 GTESTETL(:)=.FALSE.
1449
1450 !       Loop on vertical level
1451
1452 DO JK=KK,IKE-1
1453 ! IF the updraft top is reached for all column, stop the loop on levels
1454 ITEST=COUNT(GTEST)
1455 IF (ITEST==0) CYCLE
1456 !       Computation of entrainment and detrainment with KF90
1457 !       parameterization in clouds and LR01 in subcloud layer
1458
1459
1460 ! to find the LCL (check if JK is LCL or not)
1461
1462   WHERE(GTEST)
1463   WHERE ((PRC_UP(:,JK)>0.).AND.(.NOT.(GTESTLCL)))
1464       KKLCL(:) = JK
1465       GTESTLCL(:)=.TRUE.
1466   ENDWHERE
1467   ENDWHERE
1468
1469
1470 ! COMPUTE PENTR and PDETR at mass level JK
1471
1472  CALL COMPUTE_ENTR_DETR(GTEST,GTESTLCL,YFRAC_ICE,ZFRAC_ICE(:,JK),JK,&
1473                        PPABSM(:,:),PZZ(:,:),PDZZ(:,:),ZTHVM(:,:),&
1474                        PTHLM(:,JK),PRTM(:,JK),ZW_UP2(:,:),   &
1475                        PTHL_UP(:,JK),PRT_UP(:,JK),ZLUP(:), &
1476                        PENTR(:,JK),PDETR(:,JK),ZBUO_INTEG)
1477
1478
1479   IF (JK==KK) THEN
1480        PDETR(:,JK)=0.
1481   ENDIF
1482
1483
1484 !       Computation of updraft characteristics at level JK+1
1485   WHERE(GTEST)
1486     ZMIX1(:)=0.5*(PZZ(:,JK+1)-PZZ(:,JK))*(PENTR(:,JK)-PDETR(:,JK))
1487     PEMF(:,JK+1)=PEMF(:,JK)*EXP(2*ZMIX1(:))
1488   ENDWHERE
1489
1490
1491
1492 ! stop the updraft if MF becomes negative
1493   WHERE (GTEST.AND.(PEMF(:,JK+1)<=0.))
1494     PEMF(:,JK+1)=0.
1495     GTEST(:)=.FALSE.
1496     KKCTL(:) = JK+1
1497   ENDWHERE
1498
1499
1500 ! If the updraft did not stop, compute cons updraft characteritics at jk+1
1501   WHERE(GTEST)
1502     ZMIX2(:) = (PZZ(:,JK+1)-PZZ(:,JK))*PENTR(:,JK) !&
1503     ZMIX3(:) = (PZZ(:,JK+1)-PZZ(:,JK))*PDETR(:,JK) !&
1504
1505     PTHL_UP(:,JK+1) = (PTHL_UP(:,JK)*(1.-0.5*ZMIX2(:)) + PTHLM(:,JK)*ZMIX2(:)) &
1506                           /(1.+0.5*ZMIX2(:))
1507     PRT_UP(:,JK+1) = (PRT_UP (:,JK)*(1.-0.5*ZMIX2(:)) + PRTM(:,JK)*ZMIX2(:))   &
1508                           /(1.+0.5*ZMIX2(:))
1509   ENDWHERE
1510
1511
1512   IF(OMIXUV) THEN
1513     WHERE(GTEST)
1514       PU_UP(:,JK+1) = (PU_UP (:,JK)*(1-0.5*ZMIX2(:)) + PUM(:,JK)*ZMIX2(:)+ &
1515                         0.5*XPRES_UV*(PZZ(:,JK+1)-PZZ(:,JK))*&
1516                         ((PUM(:,JK+1)-PUM(:,JK))/PDZZ(:,JK+1)+&
1517                          (PUM(:,JK)-PUM(:,JK-1))/PDZZ(:,JK))        )   &
1518                         /(1+0.5*ZMIX2(:))
1519       PV_UP(:,JK+1) = (PV_UP (:,JK)*(1-0.5*ZMIX2(:)) + PVM(:,JK)*ZMIX2(:)+ &
1520                         0.5*XPRES_UV*(PZZ(:,JK+1)-PZZ(:,JK))*&
1521                         ((PVM(:,JK+1)-PVM(:,JK))/PDZZ(:,JK+1)+&
1522                          (PVM(:,JK)-PVM(:,JK-1))/PDZZ(:,JK))    )   &
1523                         /(1+0.5*ZMIX2(:))
1524     ENDWHERE
1525   ENDIF
1526
1527 ! Compute non cons. var. at level JK+1
1528   CALL TH_R_FROM_THL_RT_2D(YFRAC_ICE,ZFRAC_ICE(:,JK+1:JK+1),ZPRES_F(:,JK+1:JK+1), &
1529           PTHL_UP(:,JK+1:JK+1),PRT_UP(:,JK+1:JK+1),ZTH_UP(:,JK+1:JK+1), &
1530           PRV_UP(:,JK+1:JK+1),PRC_UP(:,JK+1:JK+1),PRI_UP(:,JK+1:JK+1))
1531
1532
1533 ! Compute the updraft theta_v, buoyancy and w**2 for level JK+1
1534   WHERE(GTEST)
1535       PTHV_UP(:,JK+1) = ZTH_UP(:,JK+1)*((1+ZRVORD*PRV_UP(:,JK+1))/(1+PRT_UP(:,JK+1)))
1536    WHERE (.NOT.(GTESTLCL))
1537       WHERE (ZBUO_INTEG(:)>0.)
1538         ZW_UP2(:,JK+1)  = ZW_UP2(:,JK) + 2.*(XABUO-XBENTR*XENTR_DRY)* ZBUO_INTEG(:)
1539       ENDWHERE
1540       WHERE (ZBUO_INTEG(:)<=0.)
1541         ZW_UP2(:,JK+1)  = ZW_UP2(:,JK) + 2.*XABUO* ZBUO_INTEG(:)
1542       ENDWHERE
1543    ENDWHERE
1544    WHERE (GTESTLCL)
1545       ZW_UP2(:,JK+1)  = ZW_UP2(:,JK)*(1.-(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:)))&
1546             /(1.+(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:))) &
1547             +2.*(XABUO)*ZBUO_INTEG/(1.+(XBDETR*ZMIX3(:)+XBENTR*ZMIX2(:)))
1548    ENDWHERE
1549  ENDWHERE
1550
1551
1552 ! Test if the updraft has reach the ETL
1553   GTESTETL(:)=.FALSE.
1554   WHERE (GTEST.AND.(ZBUO_INTEG(:)<=0.))
1555       KKETL(:) = JK+1
1556       GTESTETL(:)=.TRUE.
1557   ENDWHERE
1558
1559 ! Test is we have reached the top of the updraft
1560   WHERE (GTEST.AND.((ZW_UP2(:,JK+1)<=0.).OR.(PEMF(:,JK+1)<=0.)))
1561       ZW_UP2(:,JK+1)=0.
1562       PEMF(:,JK+1)=0.
1563       GTEST(:)=.FALSE.
1564       PTHL_UP(:,JK+1)=ZTHLM_F(:,JK+1)
1565       PRT_UP(:,JK+1)=ZRTM_F(:,JK+1)
1566       PRC_UP(:,JK+1)=0.
1567       PRV_UP(:,JK+1)=0.
1568       PTHV_UP(:,JK+1)=ZTHVM_F(:,JK+1)
1569       PFRAC_UP(:,JK+1)=0.
1570       KKCTL(:)=JK+1
1571   ENDWHERE
1572
1573 ! compute frac_up at JK+1
1574   WHERE (GTEST)
1575       PFRAC_UP(:,JK+1)=PEMF(:,JK+1)/(SQRT(ZW_UP2(:,JK+1))*ZRHO_F(:,JK+1))
1576   ENDWHERE
1577
1578 ! Updraft fraction must be smaller than XFRAC_UP_MAX
1579   WHERE (GTEST)
1580       PFRAC_UP(:,JK+1)=MIN(XFRAC_UP_MAX,PFRAC_UP(:,JK+1))
1581   ENDWHERE
1582
1583
1584 ! When cloudy and non-buoyant, updraft fraction must decrease
1585
1586   WHERE ((GTEST.AND.GTESTETL).AND.GTESTLCL)
1587     PFRAC_UP(:,JK+1)=MIN(PFRAC_UP(:,JK+1),PFRAC_UP(:,JK))
1588   ENDWHERE
1589
1590 ! Mass flux is updated with the new updraft fraction
1591
1592   PEMF(:,JK+1)=PFRAC_UP(:,JK+1)*SQRT(ZW_UP2(:,JK+1))*ZRHO_F(:,JK+1)
1593
1594 ENDDO
1595
1596
1597 PW_UP(:,:)=SQRT(ZW_UP2(:,:))
1598
1599 PEMF(:,KK) =0.
1600 PEMF(:,IKU)=0. ! modif WRF JP
1601 PFRAC_UP(:,IKU)=0.
1602
1603
1604 DO JI=1,SIZE(PTHM,1)
1605    ZDEPTH(JI) = (PZZ(JI,KKCTL(JI)) -  PZZ(JI,KKLCL(JI)) )
1606 END DO
1607
1608 GWORK1(:)= (GTESTLCL(:) .AND. (ZDEPTH(:) > 3000.) )
1609 GWORK2(:,:) = SPREAD( GWORK1(:), DIM=2, NCOPIES=IKU )
1610 ZCOEF(:,:) = SPREAD( (1.-(ZDEPTH(:)-3000.)/1000.), DIM=2, NCOPIES=IKU)
1611 ZCOEF=MIN(MAX(ZCOEF,0.),1.)
1612
1613 WHERE (GWORK2)
1614    PEMF(:,:)     = PEMF(:,:)     * ZCOEF(:,:)
1615    PFRAC_UP(:,:) = PFRAC_UP(:,:) * ZCOEF(:,:)
1616 ENDWHERE
1617
1618
1619 END SUBROUTINE COMPUTE_UPDRAFT
1620
1621
1622
1623 !     #################################################################
1624       SUBROUTINE MF_TURB(OMIXUV, PIMPL, PTSTEP,                       &
1625                 PTSTEP_MET,  PDZZ,                          &
1626                 PRHODJ,                                               &
1627                 PTHLM,PTHVM,PRTM,PUM,PVM,                             &
1628                 PTHLDT,PRTDT,PUDT,PVDT,                               &
1629                 PEMF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP,              &
1630                 PFLXZTHMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF       )
1631
1632 !     #################################################################
1633 !
1634 !
1635 !!****  *MF_TURB* - computes the MF_turbulent source terms for the prognostic
1636 !!                  variables.
1637 !!
1638
1639 !! --------------------------------------------------------------------------
1640
1641 !
1642 !
1643 IMPLICIT NONE
1644 !
1645 !
1646 !*      0.1  declarations of arguments
1647 !
1648 !
1649 LOGICAL,                INTENT(IN)   :: OMIXUV      ! True if mixing of momentum
1650 REAL,                   INTENT(IN)   :: PIMPL       ! degree of implicitness
1651 REAL,                 INTENT(IN)     ::  PTSTEP   ! Dynamical timestep
1652 REAL,                 INTENT(IN)     ::  PTSTEP_MET! Timestep for meteorological variables
1653 !
1654 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PDZZ        ! metric coefficients
1655
1656 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PRHODJ      ! dry density * Grid size
1657
1658 !   Conservative var. at t-dt
1659 REAL, DIMENSION(:,:,:), INTENT(IN) ::  PTHLM        ! conservative pot. temp.
1660 REAL, DIMENSION(:,:,:), INTENT(IN) ::  PRTM         ! water var.  where
1661 !  Virtual potential temperature at t-dt
1662 REAL, DIMENSION(:,:,:), INTENT(IN) ::  PTHVM
1663 !  Momentum at t-dt
1664 REAL, DIMENSION(:,:,:), INTENT(IN) ::  PUM
1665 REAL, DIMENSION(:,:,:), INTENT(IN) ::  PVM
1666 !
1667 ! Tendencies of conservative variables
1668 REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PTHLDT
1669
1670 REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PRTDT
1671 ! Tendencies of momentum
1672 REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PUDT
1673 REAL, DIMENSION(:,:,:),   INTENT(OUT) ::  PVDT
1674 ! Tendencies of scalar variables
1675
1676
1677 ! Updraft characteritics
1678 REAL, DIMENSION(:,:,:), INTENT(IN)   ::  PEMF,PTHL_UP,PTHV_UP,PRT_UP,PU_UP,PV_UP
1679 ! Fluxes
1680 REAL, DIMENSION(:,:,:), INTENT(OUT)  ::  PFLXZTHMF,PFLXZTHVMF,PFLXZRMF,PFLXZUMF,PFLXZVMF
1681
1682 !
1683 !
1684 !
1685 !-------------------------------------------------------------------------------
1686 !
1687 !       0.2  declaration of local variables
1688 !
1689
1690 REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) :: ZVARS
1691
1692 !
1693 !----------------------------------------------------------------------------
1694 !
1695 !*      1.PRELIMINARIES
1696 !         -------------
1697 !
1698 !
1699
1700 !
1701 PFLXZTHMF = 0.
1702 PFLXZRMF = 0.
1703 PFLXZTHVMF = 0.
1704 PFLXZUMF = 0.
1705 PFLXZVMF = 0.
1706 PTHLDT = 0.
1707 PRTDT = 0.
1708 PUDT = 0.
1709 PVDT = 0.
1710
1711 !
1712 !----------------------------------------------------------------------------
1713 !
1714 !*      2. COMPUTE THE MEAN FLUX OF CONSERVATIVE VARIABLES at time t-dt
1715 !          (equation (3) of Soares et al)
1716 !          + THE MEAN FLUX OF THETA_V (buoyancy flux)
1717 !          -----------------------------------------------
1718 !   ( Resulting fluxes are in flux level (w-point) as PEMF and PTHL_UP )
1719 !
1720
1721 PFLXZTHMF(:,:,:) = PEMF(:,:,:)*(PTHL_UP(:,:,:)-MZM(PTHLM(:,:,:)))
1722
1723 PFLXZRMF(:,:,:) =  PEMF(:,:,:)*(PRT_UP(:,:,:)-MZM(PRTM(:,:,:)))
1724
1725 PFLXZTHVMF(:,:,:) = PEMF(:,:,:)*(PTHV_UP(:,:,:)-MZM(PTHVM(:,:,:)))
1726
1727 PFLXZTHVMF(:,:,:) = 9.81/PTHVM(:,:,:)* PEMF(:,:,:)*(PTHV_UP(:,:,:)-MZM(PTHVM(:,:,:))) !JP
1728
1729 IF (OMIXUV) THEN
1730   PFLXZUMF(:,:,:) =  PEMF(:,:,:)*(PU_UP(:,:,:)-MZM(PUM(:,:,:)))
1731   PFLXZVMF(:,:,:) =  PEMF(:,:,:)*(PV_UP(:,:,:)-MZM(PVM(:,:,:)))
1732 ENDIF
1733 !
1734 !
1735 !----------------------------------------------------------------------------
1736 !
1737 !*      3. COMPUTE TENDENCIES OF CONSERVATIVE VARIABLES (or treated as such...)
1738 !          (implicit formulation)
1739 !          --------------------------------------------
1740 !
1741
1742 !
1743 !
1744 ! 3.1 Compute the tendency for the conservative potential temperature
1745 !     (PDZZ and flux in w-point and PRHODJ is mass point, result in mass point)
1746 !
1747 CALL TRIDIAG_MASSFLUX(PTHLM,PFLXZTHMF,-PEMF,PTSTEP_MET,PIMPL,  &
1748                       PDZZ,PRHODJ,ZVARS )
1749 ! compute new flux
1750 PFLXZTHMF(:,:,:) = PEMF(:,:,:)*(PTHL_UP(:,:,:)-MZM(ZVARS(:,:,:)))
1751
1752 !!! compute THL tendency
1753 !
1754 PTHLDT(:,:,:)= (ZVARS(:,:,:)-PTHLM(:,:,:))/PTSTEP_MET
1755
1756 !
1757 ! 3.2 Compute the tendency for the conservative mixing ratio
1758 !
1759 CALL TRIDIAG_MASSFLUX(PRTM(:,:,:),PFLXZRMF,-PEMF,PTSTEP_MET,PIMPL,  &
1760                                  PDZZ,PRHODJ,ZVARS )
1761 ! compute new flux
1762 PFLXZRMF(:,:,:) =  PEMF(:,:,:)*(PRT_UP(:,:,:)-MZM(ZVARS(:,:,:)))
1763
1764 !!! compute RT tendency
1765 PRTDT(:,:,:) = (ZVARS(:,:,:)-PRTM(:,:,:))/PTSTEP_MET
1766 !
1767 ! 3.3 Compute the tendency for the (non conservative but treated as it)  mixing ratio
1768 !
1769 !CALL TRIDIAG_MASSFLUX(PTHVM(:,:,:),PFLXZTHVMF,-PEMF,PTSTEP,PIMPL,  &
1770 !                                PDZZ,PRHODJ,ZVARS )
1771 ! compute new flux
1772 !PFLXZTHVMF(:,:,:) =  PEMF(:,:,:)*(PTHV_UP(:,:,:)-MZM(ZVARS(:,:,:)))
1773
1774
1775
1776 IF (OMIXUV) THEN
1777 !
1778 ! 3.3 Compute the tendency for the (non conservative but treated as it) zonal momentum
1779 !     (PDZZ and flux in w-point and PRHODJ is mass point, result in mass point)
1780 !
1781
1782  CALL TRIDIAG_MASSFLUX(PUM,PFLXZUMF,-PEMF,PTSTEP,PIMPL,  &
1783                                  PDZZ,PRHODJ,ZVARS )
1784 ! compute new flux
1785 PFLXZUMF(:,:,:) = PEMF(:,:,:)*(PU_UP(:,:,:)-MZM(ZVARS(:,:,:)))
1786
1787 ! compute U tendency
1788   PUDT(:,:,:)= (ZVARS(:,:,:)-PUM(:,:,:))/PTSTEP_MET
1789
1790 !
1791 !
1792 ! 3.4 Compute the tendency for the (non conservative but treated as it for the time beiing)
1793 !                                  meridian momentum
1794 !     (PDZZ and flux in w-point and PRHODJ is mass point, result in mass point)
1795 !
1796   CALL TRIDIAG_MASSFLUX(PVM,PFLXZVMF,-PEMF,PTSTEP,PIMPL,  &
1797                                  PDZZ,PRHODJ,ZVARS )
1798 ! compute new flux
1799 PFLXZVMF(:,:,:) = PEMF(:,:,:)*(PV_UP(:,:,:)-MZM(ZVARS(:,:,:)))
1800
1801 ! compute V tendency
1802   PVDT(:,:,:)= (ZVARS(:,:,:)-PVM(:,:,:))/PTSTEP_MET
1803
1804 ENDIF
1805
1806
1807
1808 !
1809 END SUBROUTINE MF_TURB
1810
1811 FUNCTION MZM(PA)  RESULT(PMZM)
1812
1813
1814 !
1815 IMPLICIT NONE
1816 !
1817 !*       0.1   Declarations of argument and result
1818 !              ------------------------------------
1819 !
1820
1821 REAL, DIMENSION(:,:,:), INTENT(IN)                :: PA     ! variable at mass localization
1822 REAL, DIMENSION(SIZE(PA,1),SIZE(PA,2),SIZE(PA,3)) :: PMZM   ! result at flux localization
1823
1824 !        0.2   Declarations of local variables
1825 !              -------------------------------
1826 !
1827
1828 INTEGER :: JK             ! Loop index in z direction
1829 INTEGER :: IKU            ! upper bound in z direction of PA
1830 !
1831 IKU = SIZE(PA,3)
1832
1833 !DO JK=2,IKU MODIF WRF JP
1834 DO JK=2,IKU
1835    PMZM(:,:,JK)=0.5*(PA(:,:,JK)+PA(:,:,JK-1))
1836 ENDDO
1837
1838 PMZM(:,:,1)=PA(:,:,2)
1839
1840
1841 END FUNCTION MZM
1842
1843
1844 !
1845 !
1846 !     #################################################################
1847       SUBROUTINE TH_R_FROM_THL_RT_1D(HFRAC_ICE,PFRAC_ICE,PP,             &
1848                                   PTHL, PRT, PTH, PRV, PRL, PRI       )
1849 !     #################################################################
1850 !
1851 !
1852 !!****  *TH_R_FROM_THL_RT_1D* - computes the non-conservative variables
1853 !!                          from conservative variables
1854 !!
1855
1856 !!
1857 !! --------------------------------------------------------------------------
1858 !
1859
1860 !
1861 IMPLICIT NONE
1862 !
1863 !
1864 !*      0.1  declarations of arguments
1865 !
1866 CHARACTER*1           , INTENT(IN) :: HFRAC_ICE
1867 REAL, DIMENSION(:), INTENT(INOUT) :: PFRAC_ICE
1868 REAL, DIMENSION(:), INTENT(IN) :: PP          ! Pressure
1869 REAL, DIMENSION(:), INTENT(IN) :: PTHL    ! thetal (or thetal tendency) to
1870                                              ! transform into th (or tendency)
1871 REAL, DIMENSION(:),INTENT(IN)  :: PRT    ! Total mixing ratios (or tendency) to
1872                                              ! transform into rv,rc and ri
1873                                              ! (or its tendency)
1874 REAL, DIMENSION(:), INTENT(OUT):: PTH    ! th (or th_l tendency)
1875 REAL, DIMENSION(:), INTENT(OUT):: PRV    ! vapor mixing ratio (or its tendency)
1876 REAL, DIMENSION(:), INTENT(OUT):: PRL    ! vapor mixing ratio (or its tendency)
1877 REAL, DIMENSION(:), INTENT(OUT):: PRI    ! vapor mixing ratio (or its tendency)
1878
1879 !
1880 !-------------------------------------------------------------------------------
1881 !
1882 !       0.2  declaration of local variables
1883 !
1884 INTEGER                       :: II
1885 ! Loop control
1886 REAL                          :: ZCOEF
1887 REAL, DIMENSION(SIZE(PP,1))   :: ZEXN,ZFOES,ZQVSAT
1888 REAL, DIMENSION(SIZE(PTHL,1)) :: ZRVSAT,ZCPH,ZRLTEMP
1889 REAL, DIMENSION(SIZE(PTHL,1)) :: ZT,ZLOVCPEXN,ZLOSCPEXN
1890 !----------------------------------------------------------------------------
1891 !
1892 !*      1 Initialisation
1893 !         --------------
1894 !
1895 !
1896 ZCOEF = 0.8
1897 PRL(:)=0.
1898 PRI(:)=0.
1899 ZRLTEMP(:)=0.
1900 PRV(:)=PRT(:)
1901 PTH(:)=PTHL(:)
1902 ZEXN(:)=(PP(:)/XP00) ** (XRD/XCPD)
1903 !
1904 !       2 Iteration
1905 !         ---------
1906
1907     DO II=1,20
1908       ZT(:)=PTH(:)*ZEXN(:)
1909
1910       WHERE (ZT(:) > 273.15)
1911 ! warm cloud
1912       ZFOES(:) = EXP( XALPW - XBETAW/ZT(:) - XGAMW*LOG(ZT(:))  )
1913       ZQVSAT(:) = XRD/XRV*ZFOES(:)/PP(:)   &
1914                    / (1.+(XRD/XRV-1.)*ZFOES(:)/PP(:))
1915
1916       ZRVSAT(:) = (1-ZCOEF)*ZQVSAT(:)*(1+PRT(:))+(ZCOEF)*PRV(:)
1917 !      CALL COMPUTE_FRAC_ICE_1D(HFRAC_ICE,PFRAC_ICE,PP,ZT)
1918       PFRAC_ICE(:)=0.
1919       ZRLTEMP(:)=MAX(0.,PRV(:)-ZRVSAT(:))
1920       PRV(:)=PRV(:)-ZRLTEMP(:)
1921       PRL(:)=PRL(:)+PRI(:)+ZRLTEMP(:)
1922       PRI(:)    = PFRAC_ICE(:)    * (PRL(:))
1923       PRL(:)    = (1-PFRAC_ICE(:))* (PRT(:) - PRV(:))
1924       !       2.1 Cph
1925       ZCPH(:)=XCPD+ XCPV * PRV(:)+ XCL * PRL(:) + XCI * PRI(:)
1926
1927       !       2.2 L/Cp/EXN
1928       ZLOVCPEXN(:) = (XLVTT + (XCPV-XCL) * (ZT(:)-XTT)) &
1929                         /(ZCPH*ZEXN(:))
1930       ZLOSCPEXN(:) = (XLSTT + (XCPV-XCI) * (ZT(:)-XTT)) &
1931                         /(ZCPH*ZEXN(:))
1932       PTH(:)=PTHL(:)+ZLOVCPEXN*PRL(:)+ZLOSCPEXN(:)*PRI(:)
1933
1934       ELSEWHERE
1935              ! cold shallow cloud not yet coded
1936              ! probleme also for the highest level of fire case
1937         PRL(:)=0.
1938         PRI(:)=0.
1939         PRV(:)=PRT(:)
1940         PTH(:)=PTHL(:)
1941       ENDWHERE
1942     ENDDO
1943
1944 END SUBROUTINE TH_R_FROM_THL_RT_1D
1945
1946
1947 !     #################################################################
1948       SUBROUTINE TH_R_FROM_THL_RT_2D(HFRAC_ICE,PFRAC_ICE,PP,             &
1949                                   PTHL, PRT, PTH, PRV, PRL, PRI       )
1950 !     #################################################################
1951 !
1952 !
1953 !!****  *TH_R_FROM_THL_RT_2D* - computes the non-conservative variables
1954 !!                          from conservative variables
1955 !!
1956 !!
1957 !!
1958 !! --------------------------------------------------------------------------
1959
1960
1961 !
1962 IMPLICIT NONE
1963 !
1964 !
1965 !*      0.1  declarations of arguments
1966 !
1967 CHARACTER*1         , INTENT(IN) :: HFRAC_ICE
1968 REAL, DIMENSION(:,:), INTENT(INOUT) :: PFRAC_ICE
1969 REAL, DIMENSION(:,:), INTENT(IN) :: PP          ! Pressure
1970 REAL, DIMENSION(:,:), INTENT(IN) :: PTHL    ! thetal (or thetal tendency) to
1971                                              ! transform into th (or tendency)
1972 REAL, DIMENSION(:,:),INTENT(IN)  :: PRT    ! Total mixing ratios (or tendency) to
1973                                              ! transform into rv,rc and ri
1974                                              ! (or its tendency)
1975 REAL, DIMENSION(:,:), INTENT(OUT):: PTH    ! th (or th_l tendency)
1976 REAL, DIMENSION(:,:), INTENT(OUT):: PRV    ! vapor mixing ratio (or its tendency)
1977 REAL, DIMENSION(:,:), INTENT(OUT):: PRL    ! vapor mixing ratio (or its tendency)
1978 REAL, DIMENSION(:,:), INTENT(OUT):: PRI    ! vapor mixing ratio (or its tendency)
1979
1980 !
1981 !-------------------------------------------------------------------------------
1982 !
1983 !       0.2  declaration of local variables
1984 !
1985 INTEGER                :: II
1986 ! Loop control
1987 REAL                                         :: ZCOEF
1988 REAL, DIMENSION(SIZE(PP,1),SIZE(PP,2))       :: ZEXN,ZFOES,ZQVSAT
1989 REAL, DIMENSION(SIZE(PTHL,1),SIZE(PTHL,2))   :: ZRVSAT,ZCPH,ZRLTEMP
1990 REAL, DIMENSION(SIZE(PTHL,1),SIZE(PTHL,2))   :: ZT,ZLOVCPEXN,ZLOSCPEXN
1991 !----------------------------------------------------------------------------
1992 !
1993 !*      1 Initialisation
1994 !         --------------
1995 !
1996 !
1997 ZCOEF = 0.8
1998 PRL(:,:)=0.
1999 PRI(:,:)=0.
2000 ZRLTEMP(:,:)=0.
2001 PRV(:,:)=PRT(:,:)
2002 PTH(:,:)=PTHL(:,:)
2003 ZEXN(:,:)=(PP(:,:)/XP00) ** (XRD/XCPD)
2004 !
2005 !       2 Iteration
2006 !         ---------
2007
2008     DO II=1,20
2009       ZT(:,:)=PTH(:,:)*ZEXN(:,:)
2010       WHERE (ZT(:,:) > 273.15)
2011 ! warm cloud
2012
2013       ZFOES(:,:) = EXP( XALPW - XBETAW/ZT(:,:) - XGAMW*LOG(ZT(:,:))  )
2014       ZQVSAT(:,:) = XRD/XRV*ZFOES(:,:)/PP(:,:)   &
2015                    / (1.+(XRD/XRV-1.)*ZFOES(:,:)/PP(:,:))
2016       ZRVSAT(:,:) = (1-ZCOEF)*ZQVSAT(:,:)*(1+PRT(:,:))+(ZCOEF)*PRV(:,:)
2017
2018 !      CALL COMPUTE_FRAC_ICE_2D(HFRAC_ICE,PFRAC_ICE,PP,ZT)
2019       PFRAC_ICE(:,:) = 0.
2020       ZRLTEMP(:,:)=MAX(0.,PRV(:,:)-ZRVSAT(:,:))
2021       PRV(:,:)=PRV(:,:)-ZRLTEMP(:,:)
2022       PRL(:,:)=PRL(:,:)+PRI(:,:)+ZRLTEMP(:,:)
2023       PRI(:,:)    = PFRAC_ICE(:,:)    * (PRL(:,:))
2024       PRL(:,:)    = (1-PFRAC_ICE(:,:))* (PRT(:,:) - PRV(:,:))
2025       !       2.1 Cph
2026         ZCPH(:,:)=XCPD+ XCPV * PRV(:,:)+ XCL * PRL(:,:) + XCI * PRI(:,:)
2027
2028         !       2.2 L/Cp/EXN
2029         ZLOVCPEXN(:,:) = (XLVTT + (XCPV-XCL) * (ZT(:,:)-XTT)) &
2030                         /(ZCPH*ZEXN(:,:))
2031         ZLOSCPEXN(:,:) = (XLSTT + (XCPV-XCI) * (ZT(:,:)-XTT)) &
2032                         /(ZCPH*ZEXN(:,:))
2033         PTH(:,:)=PTHL(:,:)+ZLOVCPEXN*PRL(:,:)+ZLOSCPEXN(:,:)*PRI(:,:)
2034       ELSEWHERE
2035         ! cold shallow cloud not yet coded
2036         ! probleme also for the highest level of fire case
2037         PRL(:,:)=0.
2038         PRI(:,:)=0.
2039         PRV(:,:)=PRT(:,:)
2040         PTH(:,:)=PTHL(:,:)
2041        ENDWHERE
2042     ENDDO
2043
2044 END SUBROUTINE TH_R_FROM_THL_RT_2D
2045
2046
2047
2048 !     #################################################################
2049       SUBROUTINE THL_RT_FROM_TH_R_MF( KRR,KRRL,KRRI,                  &
2050                                        PTH, PR, PLVOCPEXN, PLSOCPEXN, &
2051                                        PTHL, PRT                      )
2052 !     #################################################################
2053 !
2054 !!
2055 !!****  *THL_RT_FROM_TH_R* - computes the conservative variables THL and RT
2056 !!                           from TH and the non precipitating water species
2057 !!
2058 !! --------------------------------------------------------------------------
2059 !
2060 !*      0. DECLARATIONS
2061 !          ------------
2062 !
2063 !USE MODD_CST
2064 !
2065 IMPLICIT NONE
2066 !
2067 !
2068 !*      0.1  declarations of arguments
2069 !
2070 INTEGER,                INTENT(IN)   :: KRR           ! number of moist var.
2071 INTEGER,                INTENT(IN)   :: KRRL          ! number of liquid water var.
2072 INTEGER,                INTENT(IN)   :: KRRI          ! number of ice water var.
2073
2074 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PTH      ! theta
2075 REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PR       ! water species
2076 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PLVOCPEXN, PLSOCPEXN      ! L/(cp*Pi)
2077
2078 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PTHL     ! th_l
2079 REAL, DIMENSION(:,:,:), INTENT(OUT)  :: PRT      ! total non precip. water
2080 !
2081 !-------------------------------------------------------------------------------
2082 !
2083 !       0.2  declaration of local variables
2084 !
2085
2086 !----------------------------------------------------------------------------
2087
2088 !----------------------------------------------------------------------------
2089 !
2090 !
2091
2092 IF ( KRRL >= 1 ) THEN
2093   IF ( KRRI >= 1 ) THEN
2094     ! Rnp
2095     PRT(:,:,:)  = PR(:,:,:,1)  + PR(:,:,:,2)  + PR(:,:,:,4)
2096     ! Theta_l
2097     PTHL(:,:,:)  = PTH(:,:,:)  - PLVOCPEXN(:,:,:) * PR(:,:,:,2) &
2098                                - PLSOCPEXN(:,:,:) * PR(:,:,:,4)
2099   ELSE
2100     ! Rnp
2101     PRT(:,:,:)  = PR(:,:,:,1)  + PR(:,:,:,2)
2102     ! Theta_l
2103     PTHL(:,:,:) = PTH(:,:,:)  - PLVOCPEXN(:,:,:) * PR(:,:,:,2)
2104   END IF
2105 ELSE
2106     ! Rnp = rv
2107 PRT(:,:,:)  = PR(:,:,:,1)
2108     ! Theta_l = Theta
2109 PTHL(:,:,:) = PTH(:,:,:)
2110 END IF
2111
2112 END SUBROUTINE THL_RT_FROM_TH_R_MF
2113
2114
2115 !      #################################################
2116        SUBROUTINE TRIDIAG_MASSFLUX(PVARM,PF,PDFDT,PTSTEP,PIMPL,  &
2117                                  PDZZ,PRHODJ,PVARP             )
2118 !      #################################################
2119 !
2120 !
2121 !!****   *TRIDIAG_MASSFLUX* - routine to solve a time implicit scheme
2122 !!
2123 !!
2124 !! ---------------------------------------------------------------------
2125 !
2126 !*       0. DECLARATIONS
2127 !
2128 !
2129 !
2130 IMPLICIT NONE
2131 !
2132 !
2133 !*       0.1 declarations of arguments
2134 !
2135 REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARM   ! variable at t-1      at mass point
2136 REAL, DIMENSION(:,:,:), INTENT(IN) :: PF      ! flux in dT/dt=-dF/dz at flux point
2137 REAL, DIMENSION(:,:,:), INTENT(IN) :: PDFDT   ! dF/dT                at flux point
2138 REAL,                   INTENT(IN) :: PTSTEP  ! Double time step
2139 REAL,                   INTENT(IN) :: PIMPL   ! implicit weight
2140 REAL, DIMENSION(:,:,:), INTENT(IN) :: PDZZ    ! Dz                   at flux point
2141 REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ  ! (dry rho)*J          at mass point
2142 !
2143 REAL, DIMENSION(:,:,:), INTENT(OUT):: PVARP   ! variable at t+1      at mass point
2144 !
2145 !
2146 !*       0.2 declarations of local variables
2147 !
2148 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZRHODJ_DFDT_O_DZ
2149 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZMZM_RHODJ
2150 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZA, ZB, ZC
2151 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZY ,ZGAM
2152                                          ! RHS of the equation, 3D work array
2153 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2))                :: ZBET
2154                                          ! 2D work array
2155 INTEGER                              :: JK            ! loop counter
2156 INTEGER                              :: IKB,IKE       ! inner vertical limits
2157 !
2158 ! ---------------------------------------------------------------------------
2159 !
2160 !*      1.  Preliminaries
2161 !           -------------
2162 !
2163 IKB=1
2164 IKE=SIZE(PVARM,3)-1
2165 !
2166 ZMZM_RHODJ  = MZM(PRHODJ)
2167 ZRHODJ_DFDT_O_DZ = ZMZM_RHODJ*PDFDT/PDZZ
2168 !
2169 ZA=0.
2170 ZB=0.
2171 ZC=0.
2172 ZY=0.
2173 !
2174 !
2175 !*      2.  COMPUTE THE RIGHT HAND SIDE
2176 !           ---------------------------
2177 !
2178 ZY(:,:,IKB) = PRHODJ(:,:,IKB)*PVARM(:,:,IKB)/PTSTEP             &
2179     - ZMZM_RHODJ(:,:,IKB+1) * PF(:,:,IKB+1)/PDZZ(:,:,IKB+1)     &
2180     + ZMZM_RHODJ(:,:,IKB  ) * PF(:,:,IKB  )/PDZZ(:,:,IKB  )     &
2181     + ZRHODJ_DFDT_O_DZ(:,:,IKB+1) * 0.5*PIMPL * PVARM(:,:,IKB+1)    &
2182     + ZRHODJ_DFDT_O_DZ(:,:,IKB+1) * 0.5*PIMPL * PVARM(:,:,IKB  )
2183 !
2184 DO JK=IKB+1,IKE-1
2185   ZY(:,:,JK) = PRHODJ(:,:,JK)*PVARM(:,:,JK)/PTSTEP          &
2186     - ZMZM_RHODJ(:,:,JK+1) * PF(:,:,JK+1)/PDZZ(:,:,JK+1)    &
2187     + ZMZM_RHODJ(:,:,JK  ) * PF(:,:,JK  )/PDZZ(:,:,JK  )    &
2188     + ZRHODJ_DFDT_O_DZ(:,:,JK+1) * 0.5*PIMPL * PVARM(:,:,JK+1)  &
2189     + ZRHODJ_DFDT_O_DZ(:,:,JK+1) * 0.5*PIMPL * PVARM(:,:,JK  )  &
2190     - ZRHODJ_DFDT_O_DZ(:,:,JK  ) * 0.5*PIMPL * PVARM(:,:,JK  )  &
2191     - ZRHODJ_DFDT_O_DZ(:,:,JK  ) * 0.5*PIMPL * PVARM(:,:,JK-1)
2192 END DO
2193 !
2194 ZY(:,:,IKE) = PRHODJ(:,:,IKE)*PVARM(:,:,IKE)/PTSTEP         &
2195     - ZMZM_RHODJ(:,:,IKE+1) * PF(:,:,IKE+1)/PDZZ(:,:,IKE+1) &
2196     + ZMZM_RHODJ(:,:,IKE  ) * PF(:,:,IKE  )/PDZZ(:,:,IKE  ) &
2197     - ZRHODJ_DFDT_O_DZ(:,:,IKE ) * 0.5*PIMPL * PVARM(:,:,IKE  ) &
2198     - ZRHODJ_DFDT_O_DZ(:,:,IKE ) * 0.5*PIMPL * PVARM(:,:,IKE-1)
2199 !
2200 !
2201 !*       3.  INVERSION OF THE TRIDIAGONAL SYSTEM
2202 !            -----------------------------------
2203 !
2204 IF ( PIMPL > 1.E-10 ) THEN
2205 !
2206 !*       3.1 arrays A, B, C
2207 !            --------------
2208 !
2209   ZB(:,:,IKB) =   PRHODJ(:,:,IKB)/PTSTEP                   &
2210                 + ZRHODJ_DFDT_O_DZ(:,:,IKB+1) * 0.5*PIMPL
2211   ZC(:,:,IKB) =   ZRHODJ_DFDT_O_DZ(:,:,IKB+1) * 0.5*PIMPL
2212
2213   DO JK=IKB+1,IKE-1
2214     ZA(:,:,JK) = - ZRHODJ_DFDT_O_DZ(:,:,JK  ) * 0.5*PIMPL
2215     ZB(:,:,JK) =   PRHODJ(:,:,JK)/PTSTEP                   &
2216                  + ZRHODJ_DFDT_O_DZ(:,:,JK+1) * 0.5*PIMPL &
2217                  - ZRHODJ_DFDT_O_DZ(:,:,JK  ) * 0.5*PIMPL
2218     ZC(:,:,JK) =   ZRHODJ_DFDT_O_DZ(:,:,JK+1) * 0.5*PIMPL
2219   END DO
2220
2221   ZA(:,:,IKE) = - ZRHODJ_DFDT_O_DZ(:,:,IKE  ) * 0.5*PIMPL
2222   ZB(:,:,IKE) =   PRHODJ(:,:,IKE)/PTSTEP                   &
2223                 - ZRHODJ_DFDT_O_DZ(:,:,IKE  ) * 0.5*PIMPL
2224 !
2225 !*       3.2 going up
2226 !            --------
2227 !
2228   ZBET(:,:) = ZB(:,:,IKB)  ! bet = b(ikb)
2229   PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)
2230
2231   !
2232   DO JK = IKB+1,IKE-1
2233     ZGAM(:,:,JK) = ZC(:,:,JK-1) / ZBET(:,:)
2234                                                     ! gam(k) = c(k-1) / bet
2235     ZBET(:,:)    = ZB(:,:,JK) - ZA(:,:,JK) * ZGAM(:,:,JK)
2236                                                     ! bet = b(k) - a(k)* gam(k)
2237     PVARP(:,:,JK)= ( ZY(:,:,JK) - ZA(:,:,JK) * PVARP(:,:,JK-1) ) / ZBET(:,:)
2238                                         ! res(k) = (y(k) -a(k)*res(k-1))/ bet
2239   END DO
2240   ! special treatment for the last level
2241   ZGAM(:,:,IKE) = ZC(:,:,IKE-1) / ZBET(:,:)
2242                                                     ! gam(k) = c(k-1) / bet
2243   ZBET(:,:)     = ZB(:,:,IKE) - ZA(:,:,IKE) * ZGAM(:,:,IKE)
2244                                                     ! bet = b(k) - a(k)* gam(k)
2245   PVARP(:,:,IKE)= ( ZY(:,:,IKE) - ZA(:,:,IKE) * PVARP(:,:,IKE-1) ) / ZBET(:,:)
2246                                        ! res(k) = (y(k) -a(k)*res(k-1))/ bet
2247 !
2248 !*       3.3 going down
2249 !            ----------
2250 !
2251   DO JK = IKE-1,IKB,-1
2252     PVARP(:,:,JK) = PVARP(:,:,JK) - ZGAM(:,:,JK+1) * PVARP(:,:,JK+1)
2253   END DO
2254 !
2255 !
2256 ELSE
2257 !!! EXPLICIT FORMULATION
2258 !
2259   PVARP(:,:,IKB:IKE) = ZY(:,:,IKB:IKE) * PTSTEP / PRHODJ(:,:,IKB:IKE)
2260 !
2261 END IF
2262 !
2263 !
2264 !*       4.  FILL THE UPPER AND LOWER EXTERNAL VALUES
2265 !            ----------------------------------------
2266 !
2267 !PVARP(:,:,IKB-1)=PVARP(:,:,IKB) MODIF WRF JP
2268 PVARP(:,:,IKE+1)=PVARP(:,:,IKE)
2269 !
2270 !-------------------------------------------------------------------------------
2271 !
2272 END SUBROUTINE TRIDIAG_MASSFLUX
2273
2274 !
2275 !     #################################################################
2276       SUBROUTINE UPDRAFT_SOPE(KRR,KRRL,KRRI,OMIXUV,                   &
2277                          PZZ,PDZZ,PSFTH,PSFRV,PPABSM,PRHODREF,        &
2278                          PTKEM,PTHM,PRM,PTHLM,PRTM,PUM,PVM,           &
2279                          PTHL_UP,PRT_UP,PRV_UP,PU_UP,PV_UP,           &
2280                          PRC_UP,PRI_UP,PTHV_UP,PW_UP,PFRAC_UP,PEMF,   &
2281                          PDETR,PENTR,KKLCL,KKETL,KKCTL )
2282 !     #################################################################
2283 !!
2284 !!****  *UPDRAFT_SOPE* - Interfacing routine
2285 !!
2286 !! --------------------------------------------------------------------------
2287 !
2288 !*      0. DECLARATIONS
2289 !          ------------
2290 !
2291
2292 IMPLICIT NONE
2293
2294 !*                    1.1  Declaration of Arguments
2295 !
2296 !
2297 INTEGER,                INTENT(IN)   :: KRR       ! number of moist var.
2298 INTEGER,                INTENT(IN)   :: KRRL      ! number of liquid water var.
2299 INTEGER,                INTENT(IN)   :: KRRI      ! number of ice water var.
2300 LOGICAL,                INTENT(IN)   :: OMIXUV    ! True if mixing of momentum
2301
2302 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PZZ       !  Height at the flux point
2303
2304 REAL, DIMENSION(:,:,:), INTENT(IN)   :: PDZZ      ! depth between mass levels
2305
2306 REAL, DIMENSION(:,:),   INTENT(IN)   ::  PSFTH,PSFRV
2307                                             ! normal surface fluxes of theta,rv
2308 !
2309 !    prognostic variables at t- deltat
2310 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PPABSM      ! Pressure at time t-1
2311 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PRHODREF    ! dry density of the
2312                                                      ! reference state
2313 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKEM       ! TKE
2314 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PUM,PVM     ! momentum
2315
2316 !
2317 !   thermodynamical variables which are transformed in conservative var.
2318 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PTHM       ! pot. temp. = PTHLM in turb.f90
2319 REAL, DIMENSION(:,:,:,:), INTENT(IN)   ::  PRM        ! water species
2320 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PTHLM,PRTM !cons. var.
2321 REAL, DIMENSION(:,:,:),   INTENT(OUT)  ::  PTHL_UP,PRT_UP  ! updraft properties
2322 REAL, DIMENSION(:,:,:),   INTENT(OUT)  ::  PRV_UP,PRC_UP,PRI_UP,&!Thl,Rt,Rv,Rc,Ri
2323                                            PW_UP,PFRAC_UP,PEMF, &!w,Updraft Fraction, Mass Flux
2324                                            PDETR,PENTR,PTHV_UP, &!entrainment, detrainment, ThV
2325                                            PU_UP, PV_UP          !updraft wind component
2326
2327 INTEGER, DIMENSION(:,:),  INTENT(OUT)  ::  KKLCL,KKETL,KKCTL     !index for LCL,ETL,CTL
2328 !
2329 !
2330 !
2331 !                       1.2  Declaration of local variables
2332
2333 INTEGER :: IKB
2334
2335 !
2336 ! Variables to transform 3D fields in 2D fields
2337 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZPABSM,ZRHODREF
2338
2339 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZZZ,ZDZZ
2340
2341 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZTHLM,ZRTM,&
2342                                                               ZTHM,ZTKEM,&
2343                                                               ZUM,ZVM
2344
2345 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZRVM,ZRCM,ZRIM
2346
2347 REAL, DIMENSION(SIZE(PSFTH,1)*SIZE(PSFTH,2))    ::  ZSFTH,ZSFRV
2348
2349 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3))  :: ZTHL_UP,ZRT_UP,&
2350                                       ZRV_UP,ZRC_UP, ZRI_UP, &
2351                                       ZW_UP,ZU_UP,ZV_UP,ZTHV_UP,           &
2352                                       ZFRAC_UP,ZEMF_UP,ZENTR_UP,ZDETR_UP
2353
2354 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2))    :: ZFRAC_GRID
2355 INTEGER, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2)) :: JKETL,JKCTL,JKLCL
2356
2357 INTEGER :: IIU,IJU,IKU! Limit of the grid
2358 INTEGER :: J1D        ! horizontal loop counter
2359 INTEGER :: JK,JKK,JSV ! loop counters
2360 INTEGER :: JRR        ! moist loop counter
2361 REAL    :: ZRVORD     ! Rv/Rd (=1/0.622 cf glossaire)
2362
2363 !------------------------------------------------------------------------
2364
2365 !                     2. INITIALISATION
2366
2367 !                     2.1 Local variables, internal domain
2368
2369 IIU=SIZE(PTKEM,1)
2370 IJU=SIZE(PTKEM,2)
2371 !
2372 IKB = 1              ! Modif WRF JP
2373 IKU = SIZE(PTKEM,3)
2374
2375 ZRVORD = XRV / XRD
2376
2377
2378 DO JK=IKB,IKU
2379   ZZZ    (:,JK) = RESHAPE(PZZ    (:,:,JK),(/ IIU*IJU /) )
2380   ZDZZ   (:,JK) = RESHAPE(PDZZ    (:,:,JK),(/ IIU*IJU /) )
2381   ZTHM   (:,JK) = RESHAPE(PTHM   (:,:,JK),(/ IIU*IJU /) )
2382   ZTKEM  (:,JK) = RESHAPE(PTKEM  (:,:,JK),(/ IIU*IJU /) )
2383   ZPABSM (:,JK) = RESHAPE(PPABSM (:,:,JK),(/ IIU*IJU /) )
2384   ZRHODREF(:,JK) = RESHAPE(PRHODREF(:,:,JK),(/ IIU*IJU /) )
2385   ZRVM   (:,JK) = RESHAPE(PRM    (:,:,JK,1),(/ IIU*IJU /) )
2386   ZTHLM  (:,JK) = RESHAPE(PTHLM   (:,:,JK),(/ IIU*IJU /) )
2387   ZRTM   (:,JK) = RESHAPE(PRTM    (:,:,JK),(/ IIU*IJU /) )
2388   ZUM    (:,JK) = RESHAPE(PUM    (:,:,JK),(/ IIU*IJU /) )
2389   ZVM    (:,JK) = RESHAPE(PVM    (:,:,JK),(/ IIU*IJU /) )
2390 END DO
2391
2392 IF (KRRL>1) THEN
2393 DO JK=1,IKU
2394   ZRCM   (:,JK) = RESHAPE(PRM    (:,:,JK,2),(/ IIU*IJU /) )
2395 END DO
2396 ELSE
2397   ZRCM   (:,:) =0.
2398 ENDIF
2399 IF (KRRI>1) THEN
2400 DO JK=1,IKU
2401   ZRIM   (:,JK) = RESHAPE(PRM    (:,:,JK,4),(/ IIU*IJU /) )
2402 END DO
2403 ELSE
2404   ZRIM   (:,:) =0.
2405 ENDIF
2406
2407 ZSFTH(:)=RESHAPE(PSFTH(:,:),(/ IIU*IJU /) )
2408 ZSFRV(:)=RESHAPE(PSFRV(:,:),(/ IIU*IJU /) )
2409
2410 !Updraft begins at level 1 (Modif WRF)
2411 JK=IKB
2412
2413 !        6.2 compute properties of the updraft
2414 CALL COMPUTE_UPDRAFT(OMIXUV,ZZZ,ZDZZ,JK,      &
2415                      ZSFTH,ZSFRV,ZPABSM,ZRHODREF,ZUM,ZVM,ZTKEM, &
2416                      ZTHM,ZRVM,ZRCM,ZRIM,ZTHLM,ZRTM,       &
2417                      ZTHL_UP,ZRT_UP,ZRV_UP,ZRC_UP,ZRI_UP,&
2418                      ZTHV_UP,ZW_UP,ZU_UP,ZV_UP, &
2419                      ZFRAC_UP,ZEMF_UP,&
2420                      ZDETR_UP,ZENTR_UP,&
2421                      JKLCL,JKETL,JKCTL)
2422
2423
2424 PTHL_UP(:,:,:)= RESHAPE(ZTHL_UP(:,:), (/ IIU,IJU,IKU /))
2425 PRT_UP(:,:,:)=RESHAPE(ZRT_UP(:,:), (/ IIU,IJU,IKU /) )
2426 PRV_UP(:,:,:)=RESHAPE(ZRV_UP(:,:), (/ IIU,IJU,IKU /) )
2427 PRC_UP(:,:,:)=RESHAPE(ZRC_UP(:,:), (/ IIU,IJU,IKU /) )
2428 PRI_UP(:,:,:)=RESHAPE(ZRI_UP(:,:), (/ IIU,IJU,IKU /) )
2429 PW_UP(:,:,:)=RESHAPE(ZW_UP(:,:), (/ IIU,IJU,IKU /) )
2430 PU_UP(:,:,:)=RESHAPE(ZU_UP(:,:), (/ IIU,IJU,IKU /) )
2431 PV_UP(:,:,:)=RESHAPE(ZV_UP(:,:), (/ IIU,IJU,IKU /) )
2432 PEMF(:,:,:)=RESHAPE(ZEMF_UP(:,:), (/ IIU,IJU,IKU /) )
2433 PDETR(:,:,:)=RESHAPE(ZDETR_UP(:,:), (/ IIU,IJU,IKU /) )
2434 PENTR(:,:,:)=RESHAPE(ZENTR_UP(:,:), (/ IIU,IJU,IKU /) )
2435 PTHV_UP(:,:,:)=RESHAPE(ZTHV_UP(:,:), (/ IIU,IJU,IKU /) )
2436 KKETL(:,:)=RESHAPE(JKETL(:),(/ IIU,IJU/) )
2437 KKCTL(:,:)=RESHAPE(JKCTL(:),(/ IIU,IJU/) )
2438 KKLCL(:,:)=RESHAPE(JKLCL(:),(/ IIU,IJU/) )
2439 PFRAC_UP(:,:,:)=RESHAPE(ZFRAC_UP(:,:),(/ IIU,IJU,IKU /) )
2440
2441 END SUBROUTINE UPDRAFT_SOPE
2442
2443 !     #################################################################
2444       SUBROUTINE COMPUTE_MF_CLOUD(KRRL, PTHLM,PRC_UP, PFRAC_UP, PDZZ, KKLCL,&
2445                                   PRC_MF,PCF_MF                     )
2446 !     #################################################################
2447 !!
2448 !!****  *COMPUTE_MF_CLOUD* -
2449 !!       compute diagnostic subgrid cumulus cloud caracteristics
2450 !!
2451 !!    PURPOSE
2452 !!    -------
2453 !!****  The purpose of this routine is to compute the cloud fraction and
2454 !!      the mean cloud content associated with clouds described by the
2455 !!      mass flux scheme
2456 !!
2457 !
2458 !!**  METHOD
2459 !!    ------
2460 !!
2461 !!    EXTERNAL
2462 !!    --------
2463 !!
2464 !!    IMPLICIT ARGUMENTS
2465 !!    ------------------
2466 !!
2467 !!     REFERENCE
2468 !!     ---------
2469 !!
2470 !!
2471 !!     AUTHOR
2472 !!     ------
2473 !! --------------------------------------------------------------------------
2474 !
2475 !*      0. DECLARATIONS
2476 !          ------------
2477 !
2478
2479
2480 IMPLICIT NONE
2481
2482 !*                    1.1  Declaration of Arguments
2483 !
2484 !
2485 !
2486 INTEGER,                  INTENT(IN)   ::  KRRL         ! number of liquid water var.
2487                                                         ! scheme
2488 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PTHLM ! updraft characteritics
2489
2490 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRC_UP ! updraft characteritics
2491 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PFRAC_UP          ! Updraft Fraction
2492
2493 REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PDZZ
2494 INTEGER, DIMENSION(:,:),  INTENT(IN)   ::  KKLCL          ! index of updraft condensation level
2495 REAL, DIMENSION(:,:,:),   INTENT(OUT)  ::  PRC_MF, PCF_MF ! cloud content and
2496                                                           ! cloud fraction for MF scheme
2497
2498 !
2499 !                       1.2  Declaration of local variables
2500 !
2501 REAL, DIMENSION(SIZE(PTHLM,1),SIZE(PTHLM,2),SIZE(PTHLM,3)) ::  ZFLXZ
2502 INTEGER  :: IKU, IKB, IKE, JI,JJ,JK
2503 !------------------------------------------------------------------------
2504
2505 !                     1. INITIALISATION
2506
2507 !                     2.1 Internal domain
2508
2509 ! Internal Domain
2510 IKU=SIZE(PTHLM,3)
2511 IKB=1
2512 IKE=IKU-1
2513
2514 PCF_MF = 0.
2515 PRC_MF = 0.
2516
2517
2518 !
2519 !                    Direct cloud scheme
2520 !        This scheme may be activated only if the selected updraft model
2521 !        gives the updraft fraction as an output
2522 !
2523
2524 ! attention, les variables de l'updraft sont en niveaux flux
2525 ! ils faut les passer aux niveaux masse pour calculer PRC_MF et PCF_MF
2526 DO JI=1,SIZE(PCF_MF,1)
2527 DO JJ=1,SIZE(PCF_MF,2)
2528
2529 DO JK=KKLCL(JI,JJ),IKE
2530
2531 PCF_MF(JI,JJ,JK ) = XKCF_MF *0.5* (       &
2532   &    PFRAC_UP(JI,JJ,JK) +  PFRAC_UP(JI,JJ,JK+1) )
2533 PRC_MF(JI,JJ,JK)  = 0.5* XKCF_MF * ( PFRAC_UP(JI,JJ,JK)*PRC_UP(JI,JJ,JK)  &
2534                          + PFRAC_UP(JI,JJ,JK+1)*PRC_UP(JI,JJ,JK+1) )
2535
2536
2537 END DO
2538 END DO
2539 END DO
2540
2541
2542 END SUBROUTINE COMPUTE_MF_CLOUD
2543
2544
2545 !     #################################################################
2546       SUBROUTINE mfshconvpblinit(massflux_EDKF, entr_EDKF, detr_EDKF &
2547                                      ,thl_up, thv_up, rt_up       &
2548                                      ,rv_up, rc_up, u_up, v_up    &
2549                                      ,frac_up,RESTART,ALLOWED_TO_READ,     &
2550      &                      IDS,IDE,JDS,JDE,KDS,KDE,                    &
2551      &                      IMS,IME,JMS,JME,KMS,KME,                    &
2552      &                      ITS,ITE,JTS,JTE,KTS,KTE                 )
2553 !     #################################################################
2554
2555       IMPLICIT NONE
2556
2557       LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
2558       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,                    &
2559      &                      IMS,IME,JMS,JME,KMS,KME,                    &
2560      &                      ITS,ITE,JTS,JTE,KTS,KTE
2561
2562       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT),OPTIONAL ::             &
2563      &                                        MASSFLUX_EDKF, ENTR_EDKF, DETR_EDKF &
2564      &                                       ,THL_UP, THV_UP, RT_UP, RV_UP        &
2565      &                                       ,RC_UP, U_UP, V_UP, FRAC_UP
2566
2567       INTEGER :: I,J,K,ITF,JTF,KTF
2568
2569 !---------------------------------------------
2570       JTF=MIN0(JTE,JDE-1)
2571       KTF=MIN0(KTE,KDE-1)
2572       ITF=MIN0(ITE,IDE-1)
2573
2574 !     IF(.NOT.RESTART)THEN
2575       IF( PRESENT (MASSFLUX_EDKF) ) THEN
2576         DO J=JTS,JTF
2577         DO K=KTS,KTF
2578         DO I=ITS,ITF
2579           MASSFLUX_EDKF(I,K,J)=0.
2580           ENTR_EDKF(I,K,J)=0.
2581           DETR_EDKF(I,K,J)=0.
2582           THL_UP(I,K,J)=0.
2583           THV_UP(I,K,J)=0.
2584           RT_UP(I,K,J)=0.
2585           RV_UP(I,K,J)=0.
2586           RC_UP(I,K,J)=0.
2587           U_UP(I,K,J)=0.
2588           V_UP(I,K,J)=0.
2589           FRAC_UP(I,K,J)=0.
2590         ENDDO
2591         ENDDO
2592         ENDDO
2593       ENDIF
2594
2595 END SUBROUTINE mfshconvpblinit
2596
2597
2598 !     #########################################################
2599       SUBROUTINE BL89(PZZ,PDZZ,PTHVREF,PTHLM,KRR, &
2600                  PRM,PTKEM,PLM)
2601 !     #########################################################
2602 !
2603 !!****  *BL89* -
2604 !!
2605 !!    PURPOSE
2606 !!    -------
2607 !!    This routine computes the mixing length from Bougeault-Lacarrere 89
2608 !!    formula.
2609 !!
2610 !!**  METHOD
2611 !!    ------
2612 !!
2613 !*       0.1   Declaration of arguments
2614 !              ------------------------
2615 !
2616 INTEGER,                INTENT(IN)   :: KRR
2617
2618 REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PZZ
2619 REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PDZZ
2620 REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PTHVREF
2621 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTKEM       ! TKE
2622 !   thermodynamical variables PTHLM=Theta at the begining
2623 REAL, DIMENSION(:,:,:),   INTENT(IN) ::  PTHLM     ! conservative pot. temp.
2624 REAL, DIMENSION(:,:,:,:), INTENT(IN) ::  PRM       ! water var.
2625 REAL, DIMENSION(:,:,:),   INTENT(OUT) :: PLM       ! Mixing length
2626
2627 !*       0.2   Declaration of local variables
2628 !              ------------------------------
2629 !
2630 INTEGER :: IKB,IKE
2631
2632 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZRTM
2633 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) ::ZVPT  ! Virtual Potential Temp at half levels
2634 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2)) ::  ZLWORK
2635 !           ! downwards then upwards vertical displacement,
2636 REAL, DIMENSION(SIZE(PTKEM,1)*SIZE(PTKEM,2),SIZE(PTKEM,3)) :: ZZZ,ZDZZ,&
2637                                                               ZG_O_THVREF, &
2638                                                               ZTHM,ZTKEM,ZLM,&
2639                                                               ZLMUP,ZLMDN,&
2640                                                               ZLMTEST
2641 !           ! input and output arrays packed according one horizontal coord.
2642 REAL, DIMENSION(SIZE(PRM,1)*SIZE(PRM,2),SIZE(PRM,3),SIZE(PRM,4)) :: ZRM
2643 !           ! input array packed according one horizontal coord.
2644 REAL, DIMENSION(SIZE(PRM,1)*SIZE(PRM,2),SIZE(PRM,3)) :: ZSUM ! to replace SUM function
2645
2646 REAL   :: ZPOTE,ZLWORK1,ZLWORK2
2647
2648 INTEGER :: IIU,IJU,IKU,IICE
2649 INTEGER :: J1D        ! horizontal loop counter
2650 INTEGER :: JK,JKK     ! loop counters
2651 INTEGER :: JRR        ! moist loop counter
2652 REAL    :: ZRVORD     ! Rv/Rd (=1/0.622 cf glossaire)
2653 LOGICAL :: GUPORDN    !Test for computation of upward or downward mixing length
2654 REAL    :: Z2SQRT2
2655 !-------------------------------------------------------------------------------
2656 !
2657 Z2SQRT2=2.*SQRT(2.)
2658 IIU=SIZE(PTKEM,1)
2659 IJU=SIZE(PTKEM,2)
2660 !
2661 IKB =  1
2662 IKE = SIZE(PTKEM,3)-1
2663 IKU = SIZE(PTKEM,3)
2664 ZRVORD = XRV / XRD
2665 !-------------------------------------------------------------------------------
2666 !
2667 !*       1.    pack the horizontal dimensions into one
2668 !              ---------------------------------------
2669 !
2670 DO JK=1,IKU
2671   ZZZ    (:,JK)   = RESHAPE(PZZ    (:,:,JK),(/ IIU*IJU /) )
2672   ZDZZ   (:,JK)   = RESHAPE(PDZZ   (:,:,JK),(/ IIU*IJU /) )
2673   ZTHM   (:,JK)   = RESHAPE(PTHLM  (:,:,JK),(/ IIU*IJU /) )
2674   ZTKEM  (:,JK)   = RESHAPE(PTKEM  (:,:,JK),(/ IIU*IJU /) )
2675   ZG_O_THVREF(:,JK)   = RESHAPE(XG/PTHVREF(:,:,JK),(/ IIU*IJU /) )
2676   DO JRR=1,KRR
2677     ZRM  (:,JK,JRR) = RESHAPE(PRM    (:,:,JK,JRR),(/ IIU*IJU /) )
2678   END DO
2679 END DO
2680 !-------------------------------------------------------------------------------
2681 !
2682 !*       2.    Virtual potential temperature on the model grid
2683 !              -----------------------------------------------
2684 !
2685 IF(KRR /= 0) THEN
2686   ZSUM(:,:) = 0.
2687   DO JRR=1,KRR
2688     ZSUM(:,:) = ZSUM(:,:)+ZRM(:,:,JRR)
2689   ENDDO
2690   ZVPT(:,1:)=ZTHM(:,:) * ( 1. + ZRVORD*ZRM(:,:,1) )  &
2691                            / ( 1. + ZSUM(:,:) )
2692 ELSE
2693   ZVPT(:,1:)=ZTHM(:,:)
2694 END IF
2695 !
2696 !-------------------------------------------------------------------------------
2697 !
2698 !*       3.  loop on model levels
2699 !            --------------------
2700 !
2701 DO JK=IKB,IKE
2702 !-------------------------------------------------------------------------------
2703 !
2704 !*       4.  mixing length for a downwards displacement
2705 !            ------------------------------------------
2706    GUPORDN=.FALSE.
2707    CALL COMPUTE_BL89_ML(ZDZZ,ZTKEM,ZG_O_THVREF,ZVPT,JK,GUPORDN,ZLWORK)
2708
2709
2710
2711 !-------------------------------------------------------------------------------
2712 !
2713 !*       5.  intermediate storage of the final mixing length
2714 !            -----------------------------------------------
2715 !
2716    ZLMDN(:,JK)=MIN(ZLWORK(:),0.5*(ZZZ(:,JK)+ZZZ(:,JK+1))-ZZZ(:,IKB))
2717 !
2718 !-------------------------------------------------------------------------------
2719 !
2720 !*       6.  mixing length for an upwards displacement
2721 !            -----------------------------------------
2722
2723    GUPORDN=.TRUE.
2724    CALL COMPUTE_BL89_ML(ZDZZ,ZTKEM,ZG_O_THVREF,ZVPT,JK,GUPORDN,ZLWORK)
2725
2726    ZLMUP(:,JK)=ZLWORK(:)
2727
2728 !-------------------------------------------------------------------------------
2729 !
2730   DO J1D=1,IIU*IJU
2731     ZLWORK1=MAX(ZLMDN(J1D,JK),1.E-10)
2732     ZLWORK2=MAX(ZLMUP(J1D,JK),1.E-10)
2733     ZPOTE = ZLWORK1 / ZLWORK2
2734     ZLWORK2=1.d0 + ZPOTE**(2./3.)
2735     ZLM(J1D,JK) = Z2SQRT2*ZLWORK1/(ZLWORK2*SQRT(ZLWORK2))
2736   END DO
2737
2738 ZLM(:,JK)=( 0.5* (MAX(ZLMDN(:,JK),1.E-10)**(-2./3.)+MAX(ZLMUP(:,JK),1.E-10)**(-2./3.)) )**(-1.5)
2739 ZLM(:,JK)=MAX(ZLM(:,JK),XLINI)
2740
2741 !*       8.  end of the loop on the vertical levels
2742 !            --------------------------------------
2743 !
2744 END DO
2745 !
2746 !-------------------------------------------------------------------------------
2747 !
2748 !*       9.  boundaries
2749 !            ----------
2750 !
2751 ZLM(:,IKE)  =ZLM(:,IKE-1)
2752 ZLM(:,IKE+1)=ZLM(:,IKE-1)
2753 ZLMUP(:,IKE)  =ZLMUP(:,IKE-1)
2754 ZLMUP(:,IKE+1)=ZLMUP(:,IKE-1)
2755 ZLMDN(:,IKE)  =ZLMDN(:,IKE-1)
2756 ZLMDN(:,IKE+1)=ZLMDN(:,IKE-1)
2757 !
2758   DO JK=1,IKU
2759     PLM  (:,:,JK)   = RESHAPE(ZLM  (:,JK), (/ IIU,IJU /) )
2760   END DO
2761
2762 END SUBROUTINE BL89
2763
2764 END  MODULE MODULE_BL_MFSHCONVPBL
2765