!NCLFORTSTART
SUBROUTINE DCOMPUTEPWHI(p, tv, qv, ht, pwhi, nx, ny, nz, nzh)
    USE wrf_constants, ONLY : RD

    IMPLICIT NONE

    !f2py threadsafe
    !f2py intent(in,out) :: pwhi

    INTEGER, INTENT(IN) :: nx, ny, nz, nzh
    REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: p, tv, qv
    REAL(KIND=8), DIMENSION(nx,ny,nzh), INTENT(IN) :: ht
    REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: pwhi

!NCLEND

    INTEGER :: i, j, k
    pwhi = 0

    !$OMP PARALLEL
    DO k=1,nz
        !$OMP DO COLLAPSE(2) SCHEDULE(runtime)
        DO j=1,ny
            DO i=1,nx
                IF (tv(i,j,k) <= 273.15)    THEN
                    pwhi(i,j) = pwhi(i,j) + ((p(i,j,k)/(RD*tv(i,j,k)))*qv(i,j,k)*(ht(i,j,k+1) - ht(i,j,k)))
                ENDIF
            END DO
        END DO
        !$OMP END DO
    END DO
    !$OMP END PARALLEL

    RETURN

END SUBROUTINE DCOMPUTEPWHI


!NCLFORTSTART
SUBROUTINE DCOMPUTEPWLOW(p, tv, qv, ht, pwlow, nx, ny, nz, nzh)
    USE wrf_constants, ONLY : RD

    IMPLICIT NONE

    !f2py threadsafe
    !f2py intent(in,out) :: pwlow

    INTEGER, INTENT(IN) :: nx, ny, nz, nzh
    REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: p, tv, qv
    REAL(KIND=8), DIMENSION(nx,ny,nzh), INTENT(IN) :: ht
    REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: pwlow

!NCLEND

    INTEGER :: i, j, k
    pwlow = 0

    !$OMP PARALLEL
    DO k=1,nz
        !$OMP DO COLLAPSE(2) SCHEDULE(runtime)
        DO j=1,ny
            DO i=1,nx
                IF (tv(i,j,k) > 273.15)    THEN
                    pwlow(i,j) = pwlow(i,j) + ((p(i,j,k)/(RD*tv(i,j,k)))*qv(i,j,k)*(ht(i,j,k+1) - ht(i,j,k)))
                ENDIF
            END DO
        END DO
        !$OMP END DO
    END DO

    !$OMP END PARALLEL

    RETURN

END SUBROUTINE DCOMPUTEPWLOW
    
!NCLFORTSTART
SUBROUTINE DCOMPUTEPW(p, tv, qv, ht, pw, nx, ny, nz, nzh)
    USE wrf_constants, ONLY : RD

    IMPLICIT NONE

    !f2py threadsafe
    !f2py intent(in,out) :: pw

    INTEGER, INTENT(IN) :: nx, ny, nz, nzh
    REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: p, tv, qv
    REAL(KIND=8), DIMENSION(nx,ny,nzh), INTENT(IN) :: ht
    REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: pw

!NCLEND

    INTEGER :: i, j, k
    !REAL(KIND=8),PARAMETER :: R=287.06

    pw = 0

    !$OMP PARALLEL

    DO k=1,nz
        !$OMP DO COLLAPSE(2) SCHEDULE(runtime)
        DO j=1,ny
            DO i=1,nx
                pw(i,j) = pw(i,j) + ((p(i,j,k)/(RD*tv(i,j,k)))*qv(i,j,k)*(ht(i,j,k+1) - ht(i,j,k)))
            END DO
        END DO
        !$OMP END DO
    END DO

    !$OMP END PARALLEL

    RETURN
END SUBROUTINE DCOMPUTEPW