CALPBLREGIME.f

Go to the documentation of this file.

 
      SUBROUTINE calpblregime(PBLREGIME)
 
!
      use vrbls3d,      only: uh, vh, pmid, t, q, pint, zmid, zint
      use vrbls2d,      only: ths, qs, smstav, twbs, qwbs, pblh
      use masks,        only: dx
      use params_mod,   only: p1000, capa, d608, h1, g, rd, cp
      use ctlblk_mod,   only: jsta, jend, spval, lm, jsta_m, jend_m, im,    &
                              jsta_2l, jend_2u, ista, iend, ista_m, iend_m,ista_2l,iend_2u
      use gridspec_mod, only: gridtype
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!     
!     INCLUDE,DERIVE,SET PARAMETERS.
!     
      REAL    , PARAMETER ::  VCONVC=1.
!     
!     DECLARE VARIABLES.
!     
      REAL,dimension(ista_2l:iend_2u,jsta_2l:jend_2u),intent(inout) ::  PBLREGIME
!
      integer I,J,IE,IW,ii,jj
      real APE,THV,THVX,GOVRTH,UMASS,VMASS,WSPD,TSKV,DTHV,RHOX,fluxc,tsfc,  &
           vconv,vsgd,br,thx
     
!
!     
!*************************************************************************
!     
!     INITIALIZE ARRAYS.
!
!$omp  parallel do private(i,j)
        DO j=jsta,jend
          DO i=ista,iend
            pblregime(i,j) = spval
          ENDDO
        ENDDO
!
!     COMPUTE BULK RICHARDSON NUMBER AS CODED IN WRF module_sf_sfclay
!
!!$omp  parallel do
!!$omp& private(uhkl,ulkl,vhkl,vlkl,rib,ubot,utop,vbot,vtop,
!!$omp&         betta,ricr,ustarr,wmin,tvhtop,ztop,
!!$omp&         wndsl,wndslp,betta,ricr,ustarr,wmin 
!!$omp&       ,IFRSTLEV
!!$omp&       ,ICALPBL
!!$omp&       ,LVLP
!!$omp&       ,RIF
!!$omp&       ,RIBP
!!$omp&       ,UBOT1
!!$omp&       ,VBOT1
!!$omp&       ,ZBOT1
!!$omp&       ,THVBOT1)
!
      IF(gridtype /= 'A')THEN
         call exch(uh(1,jsta_2l,lm))
         call exch(vh(1,jsta_2l,lm))
      END IF
             
      DO j=jsta_m,jend_m
        DO i=ista_m,iend_m
!
          IF(pmid(i,j,lm)<spval .AND. qs(i,j)<spval .AND. &
             smstav(i,j)<spval) THEN 
          ape  = (p1000/pmid(i,j,lm))**capa
          thx  = t(i,j,lm)*ape
          thvx = (q(i,j,lm)*d608+h1)*thx
          govrth = g/thx
          IF(gridtype == 'E')THEN
            ie=i+mod(j+1,2) 
            iw=i+mod(j+1,2)-1
            umass = (uh(i,j-1,lm)+uh(iw,j,lm)+uh(ie,j,lm)              &   
                  +  uh(i,j+1,lm))/4.0
            vmass = (vh(i,j-1,lm)+vh(iw,j,lm)+vh(ie,j,lm)              &
                  +  vh(i,j+1,lm))/4.0
            wspd= sqrt(umass*umass+vmass*vmass)
          ELSE IF(gridtype == 'B')THEN
            ie = i
            iw = i-1
            umass = (uh(iw,j-1,lm)+uh(iw,j,lm)+uh(ie,j-1,lm)              &   
                  +  uh(i,j,lm))/4.0
            vmass = (vh(iw,j-1,lm)+vh(iw,j,lm)+vh(ie,j-1,lm)              &
                  +  vh(i,j,lm))/4.0
            wspd= sqrt(umass*umass+vmass*vmass)  
          ELSE
            wspd = sqrt(uh(i,j,lm)*uh(i,j,lm)+vh(i,j,lm)*vh(i,j,lm))
          END IF
                                                                                 
          tskv = ths(i,j)*(1.+d608*qs(i,j)*smstav(i,j))
          dthv = (thvx-tskv)
!  Convective velocity scale Vc and subgrid-scale velocity Vsg
!  following Beljaars (1995, QJRMS) and Mahrt and Sun (1995, MWR)
!                                ... HONG Aug. 2001
!
          rhox  = pint(i,j,lm+1)/rd/(t(i,j,lm)*(q(i,j,lm)*d608+h1)) !density
          fluxc = max(-twbs(i,j)/rhox/cp - d608*tskv*qwbs(i,j)/rhox,0.)
          tsfc  = ths(i,j)*(pint(i,j,lm+1)/p1000)**capa
          vconv = vconvc*(g/tsfc*pblh(i,j)*fluxc)**.33
! VCONV comes from Beljaars only
          vsgd  = 0.32 * (max(dx(i,j)/5000.-1.,0.))**.33
          wspd  = sqrt(wspd*wspd+vconv*vconv+vsgd*vsgd)
          wspd  = max(wspd,0.1)
          br    = govrth*(zmid(i,j,lm)-zint(i,j,lm+1))*dthv/(wspd*wspd)
     
          IF(br < 0.0) THEN
            pblregime(i,j) = 4.0
          ELSE IF(br == 0.0) THEN
            pblregime(i,j) = 3.0
          ELSE IF(br < 0.2) THEN
            pblregime(i,j) = 2.0
          ELSE
            pblregime(i,j) = 1.0 
          END IF
 
!         ii=im/2
!         jj=(jsta+jend)/2
!         if(i==ii.and.j==jj)print*,'Debug: CALPBLREGIME ',i,j,br,     &  
!         PBLREGIME(I,J)
          END IF !end IF PMID 
   
        ENDDO
      ENDDO
!      
!     END OF ROUTINE.
!     
      RETURN
      END