!> @file
!> @brief Subroutine that computes boundary layer fields.
!
!> Computes constant mass mean fields
!>
!> This routine computes constant mass (boundary layer)
!> fields.  The fields are a mean over layers parameter DPBND
!> (pascals) thick.  There are NBND constant mass layers, each
!> DPBND thick starting from the surface up.  Computed boundary 
!> layer fields are pressure, temperature, specific humidity,
!> relative humidity, U and V winds, vertical velocity,
!> and precipitable water.  Given these fundamental variables
!> other fields may be computed.
!>
!> @note If you change parameter NBND in this routine 
!>       don't forget to change it also in the calling
!>       subprogram, MISCLN.
!>
!> @param[out] PBND - Layer mean pressure in NBND boundary layers (NBL).
!> @param[out] TBND - Layer mean temperature in NBL.
!> @param[out] QBND - Layer mean specific humidity in NBL.
!> @param[out] RHBND - Layer mean relative hum. (QBND/QSBND) in NBL.
!> @param[out] UBND - Layer mean U wind component in NBL.
!> @param[out] VBND - Layer mean V wind component in NBL.
!> @param[out] WBND - Layer mean W wind component in NBL.
!> @param[out] OMGBND - Layer mean vertical velocity in NBL.
!> @param[out] PWTBND - Layer precipitable water in NBL.
!> @param[out] QCNVBND - Layer moisture convergence in NBL.
!> @param[out] LVLBND - ETA layer at midpoint of NBL.
!>
!> Program History
!> - 93-01-29  RUSS TREADON
!> - 93-05-07  RUSS TREADON - ADDED DOC BLOCK AND MORE COMMENTS.
!> - 93-06-19  RUSS TREADON - ADDED LVLBND TO PARAMETER LIST.
!> - 96-03-07  MIKE BALDWIN - CHANGE PWTR CALC TO INCLUDE CLD WTR
!>                            SPEED UP CODE
!> - 98-06-16  T BLACK      - CONVERSION FROM 1-D TO 2-D
!> - 98-08-18  MIKE BALDWIN - CHANGE QSBND TO RHBND IN CALL,
!>                            COMPUTE RH OVER ICE
!> - 98-12-22  MIKE BALDWIN - BACK OUT RH OVER ICE
!> - 00-01-04  JIM TUCCILLO - MPI VERSION 
!> - 02-01-15  MIKE BALDWIN - WRF VERSION
!> - 20-11-10  JESSE MENG   - USE UPP_PHYSICS MODULE
!> - 21-08-20  Wen Meng     - Retrict computation fro undefined points.
!> - 21-09-02  Bo Cui       - Decompose UPP in X direction.
!>     
!> @author Russ Treadon W/NP2 @date 1993-01-29
      SUBROUTINE BNDLYR(PBND,TBND,QBND,RHBND,UBND,VBND,       &
                        WBND,OMGBND,PWTBND,QCNVBND,LVLBND)

!     
!
      use vrbls3d,    only: pint, q, uh, vh, pmid, t, omga, wh, cwm
      use masks,      only: lmh
      use params_mod, only: d00, gi, pq0, a2, a3, a4
      use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, modelname,      &
                            jsta_m, jend_m, im, nbnd, spval, ista_2l, iend_2u, ista_m, iend_m, ista, iend
      use physcons_post,   only: con_rd, con_rv, con_eps, con_epsm1
      use gridspec_mod, only: gridtype
      use upp_physics, only: FPVSNEW
!
      implicit none
!
!     DECLARE VARIABLES.
!
      real,PARAMETER :: DPBND=30.E2
      integer, dimension(ista:iend,jsta:jend,NBND),intent(inout) :: LVLBND
      real,    dimension(ista:iend,jsta:jend,NBND),intent(inout) :: PBND,TBND,  &
               QBND,RHBND,UBND,VBND,WBND,OMGBND,PWTBND,QCNVBND

      REAL Q1D(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),V1D(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),              &
           U1D(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),QCNV1D(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U)
!
      REAL, ALLOCATABLE :: PBINT(:,:,:),QSBND(:,:,:)
      REAL, ALLOCATABLE :: PSUM(:,:,:), QCNVG(:,:,:)
      REAL, ALLOCATABLE :: PVSUM(:,:,:),NSUM(:,:,:)
!
      integer I,J,L,IE,IW,LL,LV,LBND
      real DP,QSAT,PV1,PV2,PMV,RPSUM,RPVSUM,PMIN,PM,DELP,PMINV,DELPV
      real es
!
!*****************************************************************************
!     START BNDLYR HERE
!
      ALLOCATE (PBINT(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,NBND+1))
      ALLOCATE (QSBND(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,NBND))
      ALLOCATE (PSUM(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,NBND))
      ALLOCATE (QCNVG(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LM))
      ALLOCATE (PVSUM(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,NBND))
      ALLOCATE (NSUM(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,NBND))
!
!     LOOP OVER HORIZONTAL GRID.  AT EACH MASS POINT COMPUTE
!     PRESSURE AT THE INTERFACE OF EACH BOUNDARY LAYER.
!     
!$omp  parallel do private(i,j)
      DO J=JSTA,JEND
        DO I=ISTA,IEND
          PBINT(I,J,1) = PINT(I,J,NINT(LMH(I,J))+1)
        ENDDO
      ENDDO
!
      DO LBND=2,NBND+1
!$omp  parallel do private(i,j)
        DO J=JSTA,JEND
          DO I=ISTA,IEND
            PBINT(I,J,LBND) = PBINT(I,J,LBND-1) - DPBND
          ENDDO
        ENDDO
      ENDDO

!          COMPUTE MOISTURE CONVERGENCE FOR EVERY LEVEL
      DO L=1,LM
!$omp  parallel do private(i,j)
        DO J=JSTA_2L,JEND_2U
          DO I=ISTA_2L,IEND_2U
            Q1D(I,J) = Q(I,J,L)
            U1D(I,J) = UH(I,J,L)
            V1D(I,J) = VH(I,J,L)
           ENDDO
        ENDDO
        CALL CALMCVG(Q1D,U1D,V1D,QCNV1D)
!$omp  parallel do private(i,j)
        DO J=JSTA,JEND
          DO I=ISTA,IEND
            QCNVG(I,J,L)=QCNV1D(I,J)
          ENDDO
        ENDDO
      ENDDO

!     
!     LOOP OVER HORIZONTAL.  AT EACH MASS POINT COMPUTE 
!     MASS WEIGHTED LAYER MEAN P, T, Q, U, V, OMEGA, 
!     WAND PRECIPITABLE WATER IN EACH BOUNDARY LAYER FROM THE SURFACE UP.
!     
!!$omp+ private(dp,pm,qsat)
!!$omp  parallel do private(i,j,lbnd,l,ie,iw,dp,pm,qsat,pv1,pv2,pmv)
      DO LBND=1,NBND
!$omp  parallel do private(i,j)
        DO J=JSTA,JEND
          DO I=ISTA,IEND
            PBND(I,J,LBND)   = D00
            TBND(I,J,LBND)   = D00
            QBND(I,J,LBND)   = D00
            QSBND(I,J,LBND)  = D00
            RHBND(I,J,LBND)  = D00
            UBND(I,J,LBND)   = D00
            VBND(I,J,LBND)   = D00
            WBND(I,J,LBND)   = D00
            OMGBND(I,J,LBND) = D00
            LVLBND(I,J,LBND) = 0
            NSUM(I,J,LBND)   = 0
            PSUM(I,J,LBND)   = D00
            PVSUM(I,J,LBND)  = D00
            PWTBND(I,J,LBND) = D00
            QCNVBND(I,J,LBND)= D00
          ENDDO
        ENDDO
!
!!$omp  parallel do private(i,j,l,dp,pm,es,qsat)
        DO L=1,LM
!$omp  parallel do private(i,j,dp,pm,es,qsat)
          DO J=JSTA,JEND
            DO I=ISTA,IEND
!
              PM = PMID(I,J,L)
              IF(PM<SPVAL)THEN
              IF((PBINT(I,J,LBND)   >= PM).AND.              & 
                 (PBINT(I,J,LBND+1) <= PM)) THEN
                DP = PINT(I,J,L+1) - PINT(I,J,L)
                PSUM(I,J,LBND)   = PSUM(I,J,LBND)   + DP
                NSUM(I,J,LBND)   = NSUM(I,J,LBND)   + 1
                LVLBND(I,J,LBND) = LVLBND(I,J,LBND) + L
                TBND(I,J,LBND)   = TBND(I,J,LBND)   + T(I,J,L)*DP
                QBND(I,J,LBND)   = QBND(I,J,LBND)   + Q(I,J,L)*DP
                OMGBND(I,J,LBND) = OMGBND(I,J,LBND) + OMGA(I,J,L)*DP
                IF(gridtype == 'A')THEN
                  UBND(I,J,LBND)  = UBND(I,J,LBND) + UH(I,J,L)*DP
                  VBND(I,J,LBND)  = VBND(I,J,LBND) + VH(I,J,L)*DP
                END IF
                WBND(I,J,LBND)    = WBND(I,J,LBND) + WH(I,J,L)*DP
                QCNVBND(I,J,LBND) = QCNVBND(I,J,LBND) + QCNVG(I,J,L)*DP
                PWTBND(I,J,LBND)  = PWTBND(I,J,LBND)           &
                                  + ( Q(I,J,L)+CWM(I,J,L))*DP*GI
                IF(MODELNAME == 'GFS')THEN
                  ES   = min(FPVSNEW(T(I,J,L)),PM)
                  QSAT = CON_EPS*ES/(PM+CON_EPSM1*ES)
                ELSE
                  QSAT = PQ0/PM*EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
                END IF 
                QSBND(I,J,LBND) = QSBND(I,J,LBND) + QSAT*DP
              ENDIF
              ELSE  !undeined grids
                PBND(I,J,LBND)=SPVAL
                TBND(I,J,LBND)=SPVAL
                UBND(I,J,LBND)=SPVAL
                VBND(I,J,LBND)=SPVAL
                WBND(I,J,LBND)=SPVAL
                OMGBND(I,J,LBND)=SPVAL
                QCNVBND(I,J,LBND)=SPVAL
                PWTBND(I,J,LBND)=SPVAL
                QBND(I,J,LBND)=SPVAL
                QSBND(I,J,LBND)=SPVAL
                RHBND(I,J,LBND)=SPVAL
              ENDIF
            ENDDO
          ENDDO
        ENDDO


        IF(gridtype=='E')THEN
          CALL EXCH(PINT(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,1))
          DO L=1,LM
            CALL EXCH(PINT(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,L+1))
!$omp  parallel do private(i,j,ie,iw,dp,pv1,pv2,pmv)
            DO J=JSTA_M,JEND_M
              DO I=ISTA_M,IEND_M
                IE = I+MOD(J,2)
                IW = I+MOD(J,2)-1
                PV1 = 0.25*(PINT(IW,J,L)    + PINT(IE,J,L)    &
                           +PINT(I,J+1,L)   + PINT(I,J-1,L))
                PV2 = 0.25*(PINT(IW,J,L+1)  + PINT(IE,J,L+1)  &
                           +PINT(I,J+1,L+1) + PINT(I,J-1,L+1))
                DP  = PV2-PV1
                PMV = 0.5*(PV1+PV2)
                IF((PBINT(IW,J,LBND)>=PMV).AND.        &
                   (PBINT(IW,J,LBND+1)<=PMV)) THEN
                  PVSUM(I,J,LBND) = PVSUM(I,J,LBND) + DP
                  UBND(I,J,LBND)  = UBND(I,J,LBND)  + DP* UH(I,J,L)
                  VBND(I,J,LBND)  = VBND(I,J,LBND)  + DP*VH(I,J,L)
                ENDIF
!
              ENDDO
            ENDDO
          ENDDO
        ELSE IF (gridtype=='B')THEN
          CALL EXCH(PINT(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,1))
          DO L=1,LM
            CALL EXCH(PINT(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,L+1))
!$omp  parallel do private(i,j,ie,iw,dp,pv1,pv2,pmv)
            DO J=JSTA_M,JEND_M
              DO I=ISTA_M,IEND_M
                IE = I+1
                IW = I
                PV1 = 0.25*(PINT(IW,J,L)     + PINT(IE,J,L)      &
                           +PINT(IW,J+1,L)   + PINT(IE,J+1,L))
                PV2 = 0.25*(PINT(IW,J,L+1)   + PINT(IE,J,L+1)    &
                           +PINT(IW,J+1,L+1) + PINT(IE,J+1,L+1))
                DP  = PV2-PV1
                PMV = 0.5*(PV1+PV2)
                IF((PBINT(IW,J,LBND)>=PMV).AND.        &
                   (PBINT(IW,J,LBND+1)<=PMV)) THEN
                  PVSUM(I,J,LBND) = PVSUM(I,J,LBND)+DP
                  UBND(I,J,LBND)  = UBND(I,J,LBND)+UH(I,J,L)*DP
                  VBND(I,J,LBND)  = VBND(I,J,LBND)+VH(I,J,L)*DP
                ENDIF
!
              ENDDO
            ENDDO
          ENDDO
        END IF
       
      ENDDO            ! end of lbnd loop
!
!!$omp+ private(rpsum)
!$omp  parallel do private(i,j,lbnd,rpsum,rpvsum)
      DO LBND=1,NBND
        DO J=JSTA,JEND
          DO I=ISTA,IEND
            IF(PSUM(I,J,LBND)/=0..AND.TBND(I,J,LBND)<SPVAL)THEN
              RPSUM           = 1./PSUM(I,J,LBND)
              LVLBND(I,J,LBND)= LVLBND(I,J,LBND)/NSUM(I,J,LBND)
              PBND(I,J,LBND)  = (PBINT(I,J,LBND)+PBINT(I,J,LBND+1))*0.5
              TBND(I,J,LBND)  = TBND(I,J,LBND)*RPSUM
              QBND(I,J,LBND)  = QBND(I,J,LBND)*RPSUM
              QSBND(I,J,LBND) = QSBND(I,J,LBND)*RPSUM
              OMGBND(I,J,LBND)= OMGBND(I,J,LBND)*RPSUM
              IF(gridtype=='A')THEN
               UBND(I,J,LBND)  = UBND(I,J,LBND)*RPSUM
               VBND(I,J,LBND)  = VBND(I,J,LBND)*RPSUM
              END IF 
              WBND(I,J,LBND)     = WBND(I,J,LBND)*RPSUM
              IF(QCNVBND(I,J,LBND)<SPVAL) &
              QCNVBND(I,J,LBND)  = QCNVBND(I,J,LBND)*RPSUM
            ENDIF
          ENDDO
        ENDDO

        IF(gridtype=='E' .or. gridtype=='B')THEN
          DO J=JSTA_M,JEND_M
            DO I=ISTA_M,IEND_M
              IF(PVSUM(I,J,LBND)/=0.)THEN
                RPVSUM         = 1./PVSUM(I,J,LBND)
                UBND(I,J,LBND) = UBND(I,J,LBND)*RPVSUM
                VBND(I,J,LBND) = VBND(I,J,LBND)*RPVSUM
              ENDIF
            ENDDO
          ENDDO
        END IF 
      ENDDO
!
!  IF NO ETA MID LAYER PRESSURES FELL WITHIN A BND LYR,
!   FIND THE CLOSEST LAYER TO THE BND LYR AND ASSIGN THE VALUES THERE
!
!!$omp+ private(delp,dp,l,pm,pmin,qsat)
!$omp  parallel do private(i,j,lbnd,l,ll,ie,iw,pminv,delp,dp,pm,pmin,es, &
!$omp&                     qsat,pmv,delpv,lv)
      DO LBND=1,NBND
        DO J=JSTA,JEND
          DO I=ISTA,IEND
            IF(PSUM(I,J,LBND)==0..AND.PBND(I,J,LBND)<SPVAL)THEN
              L    = LM
              PMIN = 9999999.
              PBND(I,J,LBND) = (PBINT(I,J,LBND)+PBINT(I,J,LBND+1))*0.5
!
              DO LL=1,LM
                PM   = PMID(I,J,LL)
                DELP = ABS(PM-PBND(I,J,LBND))
                IF(DELP<PMIN)THEN
                  PMIN = DELP
                  L    = LL
                ENDIF
              ENDDO
!
              DP = PINT(I,J,L+1)-PINT(I,J,L)
              PM = PMID(I,J,L)
              LVLBND(I,J,LBND) = L
              TBND(I,J,LBND)   = T(I,J,L)
              QBND(I,J,LBND)   = Q(I,J,L)
              IF(gridtype == 'A')THEN
                UBND(I,J,LBND) = UH(I,J,L)
                VBND(I,J,LBND) = VH(I,J,L)
              END IF 
              WBND(I,J,LBND)    = WH(I,J,L)
              QCNVBND(I,J,LBND) = QCNVG(I,J,L)
              IF(MODELNAME == 'GFS' .OR. MODELNAME == 'FV3R')THEN
                ES   = FPVSNEW(T(I,J,L))
                ES   = MIN(ES,PM)
                QSAT = CON_EPS*ES/(PM+CON_EPSM1*ES)
              ELSE
                QSAT=PQ0/PM*EXP(A2*(T(I,J,L)-A3)/(T(I,J,L)-A4))
              END IF  
              QSBND(I,J,LBND)  = QSAT
              OMGBND(I,J,LBND) = OMGA(I,J,L)
              PWTBND(I,J,LBND) = (Q(I,J,L)+CWM(I,J,L))*DP*GI
            ENDIF
!
!   RH, BOUNDS CHECK
!
            IF(QSBND(I,J,LBND)/=0..AND.QBND(I,J,LBND)<SPVAL)THEN
            RHBND(I,J,LBND) = QBND(I,J,LBND)/QSBND(I,J,LBND)
            IF (RHBND(I,J,LBND)>1.0) THEN
              RHBND(I,J,LBND) = 1.0
              QBND(I,J,LBND)  = RHBND(I,J,LBND)*QSBND(I,J,LBND)
            ENDIF
            IF (RHBND(I,J,LBND)<0.01) THEN
              RHBND(I,J,LBND) = 0.01
              QBND(I,J,LBND)  = RHBND(I,J,LBND)*QSBND(I,J,LBND)
            ENDIF
            ENDIF
          ENDDO
        ENDDO
!
        IF(gridtype == 'E')THEN
          DO J=JSTA_M,JEND_M
            DO I=ISTA_M,IEND_M
              IF(PVSUM(I,J,LBND)==0.)THEN
              LV = LM
              PMINV = 9999999.
              IE = I+MOD(J,2)
              IW = I+MOD(J,2)-1
! 
!           PINT HALOS UPDATED ALREADY
!
              DO LL=1,LM
                PMV = 0.125*(PINT(IW,J,LL)    + PINT(IE,J,LL)   +      &
                             PINT(I,J+1,LL)   + PINT(I,J-1,LL)  +      &
                             PINT(IW,J,LL+1)  + PINT(IE,J,LL+1) +      &
                             PINT(I,J+1,LL+1) + PINT(I,J-1,LL+1))
                DELPV = ABS(PMV-PBND(I,J,LBND))
                IF(DELPV<PMINV)THEN
                  PMINV = DELPV
                  LV    = LL
                ENDIF
              ENDDO
!
              UBND(I,J,LBND) = UH(I,J,LV)
              VBND(I,J,LBND) = VH(I,J,LV)
            ENDIF
           ENDDO
         ENDDO     
!        END IF

        ELSE IF(gridtype=='B')THEN
          DO J=JSTA_M,JEND_M
            DO I=ISTA_M,IEND_M
              IF(PVSUM(I,J,LBND)==0.)THEN
                LV=LM
                PMINV=9999999.
                IE=I+1
                IW=I
!
!           PINT HALOS UPDATED ALREADY
!
                DO LL=1,LM
                  PMV=0.125*(PINT(IW,J,LL)+PINT(IE,J,LL)+      &
                        PINT(IW,J+1,LL)+PINT(IE,J+1,LL)+       &
                        PINT(IW,J,LL+1)+PINT(IE,J,LL+1)+       &
                        PINT(IW,J+1,LL+1)+PINT(IE,J+1,LL+1))
                  DELPV=ABS(PMV-PBND(I,J,LBND))
                  IF(DELPV<PMINV)THEN
                    PMINV=DELPV
                    LV=LL
                  ENDIF
                ENDDO

                UBND(I,J,LBND) = UH(I,J,LV)
                VBND(I,J,LBND) = VH(I,J,LV)
              ENDIF
            ENDDO
          ENDDO     
        END IF 
      ENDDO                ! end of lbnd loop
!
      DEALLOCATE (PBINT, QSBND, PSUM, PVSUM, QCNVG, NSUM)
!
!     END OF ROUTINE
!     
      RETURN
      END
!