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 =  z(k) - z_surf - 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(*):: 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


        write(*,*) 'here_hpbl'
        dz_low = z(kl) - z_surf
        hdyn = min(dz_low , 0.5E0 * ustar/cor)
        write(*,*) 'hdyn', hdyn
        kpblc = kl
        kpbld = kl

        do k = kl-1,1,-1
            dz_hdyn =  z(k) - z_surf - 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 + 1
        enddo
        kpbl = MIN0(kpblc, kpbld, kl-2)
        hpbla = z(kpbl) - z_surf
        write(*,*) 'hpbla1, kpbl' , hpbla, kpbl
    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(hpbl,nkpbl,theta,z,u,v,kl,zs, du2min)
        ! 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,v
        real,intent(in)::zs, du2min
        real,intent(out)::hpbl
        integer,intent(out)::nkpbl
        real,parameter:: g = 9.81
        real,parameter:: Ric = 0.5
        real,parameter:: upper_bound = 16000.0 !upper boundary of PBLH
        real::Rib
        real du, dv, dtvirt
        integer k,KPBL, kpbld, kpblc

        KPBL=kl
        kpbld = kl

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

                du = (u(k)) * (u(k)) + (v(k)) * (v(k))
                dtvirt = theta(k) - theta(kl)
                if(du <= du2min) du = 0.01

                Rib = (g * 2.0 / (theta(k) + theta(k+1)))  &
                        * dtvirt * (z(k) - z(kl)) / du

                kpbl=k
                if(Rib.ge.Ric) then
                    !write(*,*) 'Rib', k, rib, du, dtvirt
                        kpbl=k
                        hpbl = z(k)
                    exit
                end if
            endif
        enddo
        do k=kl-1,2,-1
            dtvirt = theta(k) - theta(kl)
            if (dtvirt > 0 .and. kpbld == kl ) kpbld = k
        end do

        kpbl  = min(kpbl, kpbld+1, kl-2)
        hpbl = z(kpbl)
        hpbl = hpbl - zs
        nkpbl = kpbl

    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