wrftladj/module_cu_du.F

   1 MODULE module_cu_du
   2
   3    USE module_wrf_error
   4
   5    REAL    , PARAMETER :: cincap = -10.
   6    REAL    , PARAMETER :: capemin = 10.
   7    REAL    , PARAMETER :: dpthmin = 1000.
   8    REAL    , PARAMETER :: alpha = 0.00002
   9    REAL    , PARAMETER :: eps = 0.5
  10    REAL    , PARAMETER :: Vfall = 5.
  11
  12 !--------------------------------------------------------------------
  13
  14 CONTAINS
  15
  16    SUBROUTINE DUCU(                                          &
  17               ids,ide, jds,jde, kds,kde                      &
  18              ,ims,ime, jms,jme, kms,kme                      &
  19              ,its,ite, jts,jte, kts,kte                      &
  20              ,DT,KTAU,DX                                     &
  21              ,rho,RAINCV,NCA, PRATEC                         &  ! add PRATEC by zhuxiao
  22              ,U,V,TH,T,W,dz8w,Z,Pcps,pi                      &
  23              ,W0AVG                                          &
  24              ,CP,RD,RV,G,XLV                                 &  ! constant variable
  25              ,EP2,SVP1,SVP2,SVP3,SVPT0                       &  ! constant variable
  26              ,STEPCU,CU_ACT_FLAG,warm_rain,CUTOP,CUBOT       &
  27              ,QV                                             &
  28             ! optionals
  29              ,RTHCUTEN,RQVCUTEN                              &
  30                                                              )
  31
  32 !-------------------------------------------------------------
  33    IMPLICIT NONE
  34 !-------------------------------------------------------------
  35    INTEGER,      INTENT(IN   ) ::                            &
  36                                   ids,ide, jds,jde, kds,kde, &
  37                                   ims,ime, jms,jme, kms,kme, &
  38                                   its,ite, jts,jte, kts,kte
  39
  40    INTEGER,      INTENT(IN   ) :: STEPCU
  41    LOGICAL,      INTENT(IN   ) :: warm_rain
  42
  43    REAL,         INTENT(IN   ) :: XLV
  44    REAL,         INTENT(IN   ) :: CP,RD,RV,G,EP2
  45    REAL,         INTENT(IN   ) :: SVP1,SVP2,SVP3,SVPT0
  46
  47    INTEGER,      INTENT(IN   ) :: KTAU
  48
  49    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         , &
  50           INTENT(IN   ) ::                                   &
  51                                                           U, &
  52                                                           V, &
  53                                                           W, &
  54                                                          TH, &
  55                                                           T, &
  56                                                          QV, &
  57                                                        dz8w, &
  58                                                           z, &
  59                                                        Pcps, &
  60                                                         rho, &
  61                                                          pi
  62 !
  63    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         , &
  64           INTENT(INOUT) ::                                   &
  65                                                       W0AVG
  66
  67    REAL,  INTENT(IN   ) :: DT, DX
  68 !
  69    REAL, DIMENSION( ims:ime , jms:jme ),                     &
  70           INTENT(INOUT) ::                           RAINCV  &
  71                                                     ,PRATEC
  72
  73    REAL,    DIMENSION( ims:ime , jms:jme ),                  &
  74             INTENT(INOUT) ::                            NCA
  75
  76    REAL, DIMENSION( ims:ime , jms:jme ),                     &
  77           INTENT(OUT) ::                              CUBOT, &
  78                                                       CUTOP
  79
  80    LOGICAL, DIMENSION( ims:ime , jms:jme ),                  &
  81           INTENT(INOUT) :: CU_ACT_FLAG
  82
  83 !
  84 ! Optional arguments
  85 !
  86
  87    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),           &
  88          OPTIONAL,                                           &
  89          INTENT(INOUT) ::                                    &
  90                                                    RTHCUTEN, &
  91                                                    RQVCUTEN
  92
  93 !
  94 ! LOCAL VARS
  95
  96    LOGICAL :: flag_qr, flag_qi, flag_qs
  97
  98    REAL, DIMENSION( kts:kte ) ::                             &
  99                                                         U1D, &
 100                                                         V1D, &
 101                                                         T1D, &
 102                                                        TH1D, &
 103                                                        DZ1D, &
 104                                                         Z1D, &
 105                                                        QV1D, &
 106                                                         P1D, &
 107                                                       RHO1D, &
 108                                                     W0AVG1D
 109
 110    REAL, DIMENSION( kts:kte )::                              &
 111                                                       DQVDT, &
 112                                                       DTHDT
 113
 114    REAL    :: PPRATE,TST,tv,PRS,RHOE,W0,SCR1,DXSQ,RTHCUMAX
 115
 116    INTEGER :: i,j,k,i_start,i_end,j_start,j_end,sz,NTST,ICLDCK
 117 !
 118    DXSQ=DX*DX
 119
 120    NTST=STEPCU
 121    ICLDCK=MOD(KTAU,NTST)
 122    IF(ICLDCK.EQ.0 .or. KTAU .eq. 1) then
 123 !
 124 !  Keep away from specified and relaxation zone (should be for just specified and nested bc)
 125    sz = 1
 126    i_start=max(ids+sz,its)
 127    i_end=min(ide-1-sz,ite)
 128    j_start=max(jds+sz,jts)
 129    j_end=min(jde-1-sz,jte)
 130 !
 131      DO J = j_start, j_end
 132        DO I= i_start, i_end
 133
 134             DO k=kts,kte
 135                DQVDT(k)=0.
 136                DTHDT(k)=0.
 137             ENDDO
 138             RAINCV(I,J)=0.
 139             PRATEC(I,J)=0.
 140             CUTOP(I,J)=KTS
 141             CUBOT(I,J)=KTE+1
 142 !
 143 ! assign vars from 3D to 1D
 144
 145             DO K=kts,kte
 146                U1D(K) =U(I,K,J)
 147                V1D(K) =V(I,K,J)
 148                T1D(K) =T(I,K,J)
 149                TH1D(K) =TH(I,K,J)
 150                RHO1D(K) =rho(I,K,J)
 151                QV1D(K)=QV(I,K,J)
 152                IF ( QV1D(K) .LT. 1.E-08 ) QV1D(K) = 1.E-08
 153                P1D(K) =Pcps(I,K,J)
 154                W0AVG1D(K) =W0AVG(I,K,J)
 155                DZ1D(k)=dz8w(I,K,J)
 156                Z1D(k)=z(I,K,J)
 157             ENDDO
 158             CALL DUCU1D(I, J,                       &
 159                  U1D,V1D,T1D,QV1D,P1D,DZ1D,Z1D,     &
 160                  W0AVG1D,DT,DX,DXSQ,RHO1D,TH1D,     &
 161                  XLV,CP,RD,RV,G,                    &
 162                  EP2,SVP1,SVP2,SVP3,SVPT0,          &
 163                  DQVDT,DTHDT,                       &
 164                  PPRATE,NCA,NTST,                   &
 165                  CUTOP,CUBOT,                       &
 166                  ids,ide, jds,jde, kds,kde,         &
 167                  ims,ime, jms,jme, kms,kme,         &
 168                  its,ite, jts,jte, kts,kte)
 169             IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
 170               DO K=kts,kte
 171                  RTHCUTEN(I,K,J)=DTHDT(K)
 172                  RQVCUTEN(I,K,J)=DQVDT(K)
 173               ENDDO
 174               PRATEC(I,J)=PPRATE
 175               RAINCV(I,J)=PPRATE*DT
 176             ENDIF
 177        ENDDO
 178      ENDDO
 179    ENDIF
 180 !
 181    END SUBROUTINE DUCU
 182 ! ****************************************************************************
 183 !-----------------------------------------------------------
 184    SUBROUTINE DUCU1D (I, J,                           &
 185                       U0,V0,T0,QV0,P0,DZQ,Z,W0AVG1D,       &
 186                       DELT,DX,DXSQ,rhoe,TH0,               &
 187                       XLV,CP,RD,RV,G,                      &
 188                       EP2,SVP1,SVP2,SVP3,SVPT0,            &
 189                       DQVDT,DTHDT,                         &
 190                       PPRATE,NCA,NTST,                     &
 191                       CUTOP,CUBOT,                         &
 192                       ids,ide, jds,jde, kds,kde,           &
 193                       ims,ime, jms,jme, kms,kme,           &
 194                       its,ite, jts,jte, kts,kte)
 195 !-----------------------------------------------------------
 196 !
 197       IMPLICIT NONE
 198 !-----------------------------------------------------------
 199       INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde, &
 200                                 ims,ime, jms,jme, kms,kme, &
 201                                 its,ite, jts,jte, kts,kte, &
 202                                 I,J,NTST
 203
 204 !
 205       REAL, DIMENSION( kts:kte ),                          &
 206             INTENT(IN   ) ::                           U0, &
 207                                                        V0, &
 208                                                        T0, &
 209                                                       TH0, &
 210                                                       QV0, &
 211                                                        P0, &
 212                                                      rhoe, &
 213                                                       DZQ, &
 214                                                         Z, &
 215                                                   W0AVG1D
 216 !
 217       REAL,  INTENT(IN   ) :: DELT,DX,DXSQ
 218 !
 219
 220       REAL,  INTENT(IN   ) :: XLV,CP,RD,RV,G
 221       REAL,  INTENT(IN   ) :: EP2,SVP1,SVP2,SVP3,SVPT0
 222
 223 !
 224       REAL, DIMENSION( kts:kte ), INTENT(INOUT) ::         &
 225                                                     DQVDT, &
 226                                                     DTHDT
 227
 228       REAL,    DIMENSION( ims:ime , jms:jme ),             &
 229             INTENT(INOUT) ::                          NCA
 230
 231       REAL, DIMENSION( ims:ime , jms:jme ),                &
 232             INTENT(OUT) ::                          CUBOT, &
 233                                                     CUTOP
 234       REAL,  INTENT(OUT  ) :: PPRATE
 235 !
 236 !...DEFINE LOCAL VARIABLES...
 237 !
 238       REAL, DIMENSION( kts:kte ) :: cond,h,hs,qs,x
 239       REAL    :: buoy,cape,cin,dh,dq,dt,dtm,ep,es, &
 240                  evap,hp,mp,qp,qsp,rrk,rrkp, &
 241                  tadp,tdp,zb,zg,zi,zt
 242       INTEGER :: ipos,isat,k,kb,ki,kt
 243 !
 244 !...DEFINE PROFILES
 245       DO k=kts,kte
 246         h(k)=cp*t0(k)+g*z(k)+xlv*qv0(k)
 247         es=1000.*svp1*EXP(svp2*(t0(k)-svpt0)/(t0(k)-svp3))
 248         qs(k)=ep2*es/(p0(k)-es)
 249         hs(k)=cp*t0(k)+g*z(k)+xlv*qs(k)
 250         x(k)=xlv*xlv*qs(k)/(cp*rv*t0(k)*t0(k))
 251         dthdt(k)=0.
 252         dqvdt(k)=0.
 253       ENDDO
 254       pprate=0.
 255       zg=z(1)-0.5*dzq(1)
 256 !
 257 !...LOOP OVER PARCELS
 258       loop_origin: DO ki=kts,kte
 259         hp=h(ki)
 260         qp=qv0(ki)
 261         mp=alpha*rhoe(ki)*dzq(ki)
 262         zi=z(ki)
 263         buoy=0.
 264         cape=0.
 265         cin=0.
 266         dtm=0.
 267         isat=0
 268         ipos=0
 269         kt=0
 270         kb=0
 271         zt=0.
 272         zb=0.
 273         cond=0.
 274 !
 275 !...LIFT PARCEL
 276         loop_lift: DO k=ki+1,kte
 277           tadp=t0(ki)+(g/cp)*(z(ki)-z(k))
 278           ep=p0(k)*qv0(ki)/(ep2+qv0(ki))
 279           tdp=(svpt0-(svp3/svp2)*ALOG(0.001*ep/svp1))/(1.-(1./svp2)*ALOG(0.001*ep/svp1))
 280           IF(tadp.GE.tdp)THEN
 281 !         unsaturated
 282             IF(isat.EQ.1)THEN
 283               print *,i,j,'sounding warning: unsat above sat'
 284             ENDIF
 285             dt=tadp-t0(k)
 286             cond(k)=0.
 287           ELSE
 288 !         saturated
 289             IF(isat.EQ.0)THEN
 290               kb=k
 291               zb=z(k)-0.5*dzq(k)
 292             ENDIF
 293             isat=1
 294             dh=hp-hs(k)
 295             dt=(dh/cp)/(1.+x(k))
 296             qsp=qs(k)+(dh/xlv)*x(k)/(1.+x(k))
 297 !...CONDENSATE PRODUCED
 298             cond(k)=mp*(qp-qsp)
 299             qp=qsp
 300           ENDIF
 301           buoy=buoy+g*dt*dzq(k)/t0(k)
 302           cape=max(cape,buoy)
 303           IF(buoy.GE.cincap)cin=min(cin,buoy)
 304           IF(dt .GE. 0.)THEN
 305             kt=k
 306             zt=z(k)+0.5*dzq(k)
 307           ELSE IF(dt .LT. 0. .AND. dtm .GE. 0.)THEN
 308 ! cloud top is level closest to parcel temperature
 309             IF(abs(dt) .LT. abs(dtm))THEN
 310               kt=k
 311               zt=z(k)+0.5*dzq(k)
 312             ENDIF
 313           ENDIF
 314           dtm=dt
 315 ! continue lifting until buoyancy is gone
 316           IF(buoy.LT.cincap)THEN
 317             EXIT loop_lift
 318           ENDIF
 319           IF(buoy.GT.0.)THEN
 320 !         positive area detected
 321             ipos=1
 322           ENDIF
 323           IF(k.EQ.1)THEN
 324             kt=k
 325             zt=z(k)+0.5*dzq(k)
 326             zi=z(ki)
 327             print *,'sounding warning: cloud top at model top'
 328           ENDIF
 329         ENDDO loop_lift
 330 !
 331 !...CHECK FOR CLOUD
 332         IF(isat.EQ.0)THEN
 333 !       no cloud from lifting - no convection
 334           CYCLE loop_origin
 335         ENDIF
 336         IF(zt-zb.LE.dpthmin)THEN
 337 !       not more than one cloud level - no convection
 338           CYCLE loop_origin
 339         ENDIF
 340         IF(ipos.EQ.0)THEN
 341 !       no buoyancy in cloud - no convection
 342           CYCLE loop_origin
 343         ENDIF
 344         IF(cape.LE.capemin)THEN
 345 !       not enough cape
 346           CYCLE loop_origin
 347         ENDIF
 348 !
 349 !...IF CHECK FOR CLOUD SUCCESSFUL
 350 !
 351 !...DETRAINMENT
 352         dh=hp-hs(kt)
 353         dt=(dh/cp)/(1.+x(kt))
 354         dq=qs(kt)+(dh/xlv)*x(kt)/(1.+x(kt))-qv0(kt)
 355         dthdt(kt)=dthdt(kt)+mp*(th0(kt)/t0(kt))*dt/(rhoe(kt)*dzq(kt))
 356         dqvdt(kt)=dqvdt(kt)+mp*dq/(rhoe(kt)*dzq(kt))
 357 !
 358 !...SUBSIDENCE
 359         loop_subsidence: DO k=kt-1,ki,-1
 360           dthdt(k)=dthdt(k)+mp*(th0(k+1)-th0(k))/(rhoe(k)*dzq(k))
 361           dqvdt(k)=dqvdt(k)+mp*(qv0(k+1)-qv0(k))/(rhoe(k)*dzq(k))
 362         ENDDO loop_subsidence
 363 !
 364 !...RAINFALL AND EVAPORATION
 365         rrkp=0.
 366         loop_rainfall: DO k=kt,1,-1
 367           rrk=rrkp+cond(k)
 368           evap=dzq(k)*rrkp/Vfall*eps*(qs(k)-qv0(k))
 369
 370 ! restrict evap to below cloud base
 371           IF(k.GE.kb) evap=0.
 372           evap=min(rrk,evap)
 373           rrk= rrk-evap
 374
 375           dqvdt(k)=dqvdt(k)+evap/(rhoe(k)*dzq(k))
 376           dthdt(k)=dthdt(k)-(xlv/cp)*(th0(kt)/t0(kt))*evap/(rhoe(k)*dzq(k))
 377           rrkp=rrk
 378         ENDDO loop_rainfall
 379         pprate=pprate+rrkp
 380
 381       ENDDO loop_origin
 382
 383 !-----------------------------------------------------------------------
 384    END SUBROUTINE DUCU1D
 385 ! ***********************************************************************
 386 !====================================================================
 387    SUBROUTINE ducuinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN,         &
 388                      RQICUTEN,RQSCUTEN,NCA,W0AVG,P_QC,P_QR,         &
 389                      SVP1,SVP2,SVP3,SVPT0,                          &
 390                      P_FIRST_SCALAR,restart,allowed_to_read,        &
 391                      ids, ide, jds, jde, kds, kde,                  &
 392                      ims, ime, jms, jme, kms, kme,                  &
 393                      its, ite, jts, jte, kts, kte                   )
 394 !--------------------------------------------------------------------
 395    IMPLICIT NONE
 396 !--------------------------------------------------------------------
 397    LOGICAL , INTENT(IN)           ::  restart,allowed_to_read
 398    INTEGER , INTENT(IN)           ::  ids, ide, jds, jde, kds, kde, &
 399                                       ims, ime, jms, jme, kms, kme, &
 400                                       its, ite, jts, jte, kts, kte
 401    INTEGER , INTENT(IN)           ::  P_QC,P_QR,P_FIRST_SCALAR
 402
 403    REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::       &
 404                                                           RTHCUTEN, &
 405                                                           RQVCUTEN, &
 406                                                           RQCCUTEN, &
 407                                                           RQRCUTEN, &
 408                                                           RQICUTEN, &
 409                                                           RQSCUTEN
 410
 411    REAL ,   DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: W0AVG
 412
 413    REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: NCA
 414
 415    INTEGER :: i, j, k, itf, jtf, ktf
 416    REAL, INTENT(IN)    :: SVP1,SVP2,SVP3,SVPT0
 417
 418    jtf=min0(jte,jde-1)
 419    ktf=min0(kte,kde-1)
 420    itf=min0(ite,ide-1)
 421
 422    IF(.not.restart)THEN
 423
 424       DO j=jts,jtf
 425       DO k=kts,ktf
 426       DO i=its,itf
 427          RTHCUTEN(i,k,j)=0.
 428          RQVCUTEN(i,k,j)=0.
 429       ENDDO
 430       ENDDO
 431       ENDDO
 432
 433       IF (P_QC .ge. P_FIRST_SCALAR) THEN
 434          DO j=jts,jtf
 435          DO k=kts,ktf
 436          DO i=its,itf
 437             RQCCUTEN(i,k,j)=0.
 438          ENDDO
 439          ENDDO
 440          ENDDO
 441       ENDIF
 442
 443       DO j=jts,jtf
 444       DO i=its,itf
 445          NCA(i,j)=-100.
 446       ENDDO
 447       ENDDO
 448
 449       DO j=jts,jtf
 450       DO k=kts,ktf
 451       DO i=its,itf
 452          W0AVG(i,k,j)=0.
 453       ENDDO
 454       ENDDO
 455       ENDDO
 456
 457    endif
 458
 459    END SUBROUTINE ducuinit
 460
 461 END MODULE module_cu_du
 462