module diag_pbl

contains
    subroutine get_hpbl_old(hpbla, kpbl, theta_v, z, z_surf, kl, cor, ustar)
        implicit none
        real, intent(in):: cor, ustar, z_surf
        integer, intent(in):: kl
        real, intent(in), dimension(KL):: theta_v, z
        real, intent(out):: hpbla
        integer, intent(out):: kpbl

        !local variabls
        real :: dz_low, hdyn, dz_hdyn, dz_conv
        integer:: k, kpbld, kpblc



        dz_low = z(kl) - z_surf
        hdyn = AMIN1(dz_low , 0.5E0 * ustar/cor)

        kpblc = kl
        kpbld = kl

        do k = kl-1,1,-1
            dz_hdyn =  dz_low - HDYN
            dz_conv = theta_v(k) - theta_v(kl)
            if(kpbld.EQ.kl .AND. dz_hdyn.GE.0.E0) kpbld = k
            if(kpblc.EQ.kl .AND. dz_conv.GT.0.E0) kpblc = k
        enddo
        kpbl = MIN0(kpblc, kpbld, KL-2)
        hpbla = z(kpbl) - z_surf
    end subroutine get_hpbl_old

    subroutine get_hpbl(hpbla, kpbl, theta_v, z, z_surf, kl, cor, ustar)
        implicit none
        real, intent(in):: cor, ustar, z_surf
        integer, intent(in):: kl
        real, intent(in), dimension(KL):: theta_v, z
        real, intent(out):: hpbla
        integer, intent(out):: kpbl

        !local variabls
        real :: dz_low, hdyn, dz_hdyn, dz_conv
        integer:: k, kpbld, kpblc



        dz_low = z(kl) - z_surf
        hdyn = AMIN1(dz_low , 0.5E0 * ustar/cor)

        kpblc = kl
        kpbld = kl

        do k = kl-1,1,-1
            dz_hdyn =  dz_low - HDYN
            dz_conv = theta_v(k) - theta_v(kl)
            if(kpbld.EQ.kl .AND. dz_hdyn.GE.0.E0) kpbld = k
            if(kpblc.EQ.kl .AND. dz_conv.GT.0.E0) kpblc = k
        enddo
        kpbl = MIN0(kpblc, kpbld, KL-2)
        hpbla = z(kpbl) - z_surf
    end subroutine get_hpbl

    subroutine diag_pblh_inmcm(z,thetav,lat,zs,ustar,kl,hpbl)
    ! This subroutine calculates height of the PBL above ground level according to INMCM algorithm
    ! input:
    ! z - array with heights above sea level
    ! zs - height of surface above sea level
    ! thetav - array with virtual potential temperatures
    ! lat - latitude
    ! ustar - ustar at the lowest model level
    ! kl - number of vertical levels
    !!! IMPORTANT:  In INMCM model vertical levels start form the top of the atmosphere, so arrays should be organized accordingly
    ! output:
    ! hpbl - height of PBL

    implicit none
    integer,intent(in)::kl
    real,intent(in),dimension(kl)::z,thetav
    real,intent(in)::lat,zs,ustar
    real,intent(out)::hpbl
    real,parameter::omega=7.2921/(10**5)
    real,parameter::cormin=5.0/(10**5)
    real cor,yzkl,hdyn,ydelz,ydeltv
    integer KPBL,KPBLC,KPBLD,k

    cor = 2 * omega * sin (lat * 3.14 / 180)
    cor = max(cormin,abs(cor))

    yzkl = z(kl)-zs
    hdyn = min(yzkl,0.5 * ustar/cor)

    ydelz = yzkl - hdyn

    if(ydelz.ge.0)then
    KPBLD = kl-1
    else
    KPBLD = kl
    endif

    KPBLC = kl
    do k=kl-1,1,-1
    ydeltv = thetav(k) - thetav(kl)
    if(KPBLC.eq.kl.and.ydeltv.gt.0) KPBLC = k
    enddo

    KPBL = min(KPBLC,KPBLD,KL-2)

    call get_hpbl(hpbl, kpbl, thetav, z, zs, kl, cor, ustar)
    !hpbl = z(KPBL) - zs
    return
    end subroutine diag_pblh_inmcm



    subroutine diag_pblh_rib(theta,z,u,kl,zs,hpbl,nkpbl)
    ! This subroutine calculates PBL depth according to (Troen and Mahrt 1986) as the lowest level where Rib>Ric
    ! Ric varies between 0.15 and 0.5
    ! Rib = (g/theta0)*(theta(z)-theta(s))*z/u(z)**2
    !input:
    ! theta - array with (virtual) potential temperature
    ! z - array with heights above sea level
    ! u - array with wind speed
    ! zs - height of surface above sea level
    ! kl - number of vertical levels
    ! output:
    ! hpbl - PBL height above ground level
    !!! IMPORTANT:  In INMCM model vertical levels start form the top of the atmosphere, so arrays should be organized accordingly
    implicit none
    integer,intent(in)::kl
    real,intent(in),dimension(kl)::z,theta,u
    real,intent(in)::zs
    real,intent(out)::hpbl
    integer,intent(out)::nkpbl
    real,parameter:: g = 9.81
    real,parameter:: Ric = 0.25
    real,parameter:: theta0 = 300.0
    real,parameter:: upper_bound = 6000 !upper boundary of PBLH
    real,dimension(kl)::Rib
    real du
    integer k,KPBL

    KPBL=kl

    do k=kl-1,1,-1
    if(z(k).lt.upper_bound)then

    du = u(k)-u(kl)
    if(du.eq.0) du = 0.1

    Rib(k) = (g / theta0) * (theta(k) - theta(kl)) * (z(k) - z(kl)) / (du**2)

    if(Rib(k).ge.Ric.and.KPBL.eq.kl)then
     KPBL = k
     if(KPBL.le.kl-2)then
     hpbl = z(k) - (z(k)-z(k+1))*(Rib(k)-Ric)/(Rib(k)-Rib(k+1))
     else
     hpbl = z(k)
     endif
    endif

    endif
    enddo

    if(KPBL.eq.kl) hpbl=z(kl)

    hpbl = hpbl - zs
    nkpbl = kl - KPBL + 1

    end subroutine diag_pblh_rib

    subroutine diag_pblh_rh(t,z,q,p,kl,zs,hpbl,nkpbl)
        ! This subroutine calculates PBL height as height where minimal gradient of relative humidity is found
        !input:
        ! t - air temperature (not potential)
        ! z - heights above sea level
        ! q - specific humidity
        ! p - pressure
        ! u - wind speed
        ! zs - height of surface above sea level
        ! kl - number of vertical levels; index of the lowest level
        ! output:
        ! hpbl - PBL height above ground level
        ! nkpbl - number of levels inside PBL
        !!! IMPORTANT:  In INMCM model vertical levels start form the top of the atmosphere, so arrays should be organized accordingly
    implicit none
        integer,intent(in)::kl
        real,intent(in),dimension(kl)::t,z,q,p
        real,intent(in)::zs
        real,intent(out)::hpbl
        integer,intent(out)::nkpbl
        real,dimension(kl):: Es,rh
        real,parameter::E0 = 611.0
        real,parameter:: upper_bound = 6000 !maximum PBLH
        integer i,j,k,KPBL
        real mindrh,drhdz

          Es(:) = E0*10**(7.45*(t(:)-273.15)/(235+t(:)-273.15))
          rh(:) = 100*q(:)*p(:)/(0.622*Es(:))

        KPBL=kl

        mindrh = 100.0
        do k=kl-1,1,-1
            if(z(k).lt.upper_bound)then

                drhdz = rh(k) - rh(k+1)

                if(drhdz.lt.mindrh)then
                    KPBL = k
                    mindrh = drhdz
                endif
            endif
        enddo

        hpbl = z(KPBL) - zs
        nkpbl = kl - KPBL + 1
    end subroutine diag_pblh_rh

end module diag_pbl