function MixedLayerDepth(Nz_in,zi_in,z_in,rho_in,thetav_in) result (mld_output)
    ! Touzé-Peiffer et al. (2022) mixed layer height algorithm (https://doi.org/10.1175/JAMC-D-21-0048.1).
    !
    ! Note that following Touzé-Peiffer et al., virtual potential temperature is based on water vapor mixing
    ! ratio r only, ignoring a contribution of liquid water: θυ = θ(1 + 0.61r). 
    !
    ! Modified: Jan Kazil (set z_min = 20,
    !                      changed depth -> height in comments, expanded comments)
    ! Modified: Peter Blossey (changed to vector-valued function,
    !                          second output is mixed-layer-average thetav), 25 Apr 2024
    ! Modified: Peter Blossey (convert to function, handle non-zero 
    !      surface altitude, integrate from z_min including partial
    !      levels, use interface altitudes as input), 14 Dec 2023
    ! Author: jan.kazil@noaa.gov 2023-07-26 15:09:13 -06:00
    implicit none

    integer, intent(in) :: Nz_in
    real, intent(in) :: zi_in(Nz_in+1) ! interface heights, zi_in(1) is surface height
    real, intent(in) :: z_in(Nz_in+1) ! model level heights
    real, intent(in) :: rho_in(Nz_in), thetav_in(Nz_in) ! rho_in(1) is at height 0.5*SUM(zi_in(1:2))
    real, dimension(2) :: mld_output ! holds mixed layer height in (1) and mixed layer mean thetav in (2).
    
    integer :: k, kb
    real :: z1, z2, f1, f2, weight_this_level, z_surface
    real :: integral(Nz_in), weight(Nz_in)

    real, parameter :: z_min = 20. !m
    real, parameter :: epsilon = 0.2 !K 

    logical, parameter :: debug = .false.

    z_surface = zi_in(1)

    mld_output = 0. ! height above surface

    weight = 0.
    integral = 0.

    do k = 1,Nz_in
       if(zi_in(k+1) < z_surface + z_min) then
          cycle
       elseif(zi_in(k+1) > z_surface + z_min) then
          kb = MAX(1,k-1)

          ! accumulate integrals of rho*dz and rho*thetav*dz
          weight_this_level = rho_in(k)*(zi_in(k+1) - MAX(z_min,zi_in(k)))
          weight(k) = weight(kb) + weight_this_level
          integral(k) = integral(kb) + weight_this_level*thetav_in(k)

          if(debug) then
             write(*,922) k, weight(k), integral(k), integral(k)/weight(k) + epsilon, thetav_in(k)
922          format('k/wgt/integral/thresh/thetav = ',I5,2E10.2,2F8.2)
          end if

          ! check for top of mixed-layer
          if(thetav_in(k).ge. integral(k)/weight(k) + epsilon) then

             if(weight(kb).eq.0.) then
                ! should not happen...
                mld_output(1) = z_in(k) - z_surface
                
                ! Return vertically-averaged theta_v within mixed layer in element 2
                mld_output(2) = integral(kb)/weight(kb)
             else
                ! interpolate between level k and k-1 to find height.
                z1 = z_in(k) ! level k
                z2 = z_in(kb) ! level k-1

                f1 = thetav_in(k) - ( integral(k)/weight(k) + epsilon )
                f2 = thetav_in(kb) - ( integral(kb)/weight(kb) + epsilon )

                ! Return height in element 1 of mld_output
                mld_output(1) = z1 - f1*(z2-z1)/(f2-f1) - z_surface

                ! Add portion of grid cell below the Mixed Layer height to integral of theta_v
                !   -- assuming uniform value of theta_v within each grid cell.
                weight_this_level = rho_in(k)*(mld_output(1) - MAX(z_min,zi_in(k)))
                weight(k) = weight(kb) + weight_this_level
                integral(k) = integral(kb) + weight_this_level*thetav_in(k)

                ! Return vertically-averaged theta_v within mixed layer in element 2
                mld_output(2) = integral(k)/weight(k)
             end if

             exit ! break out of loop, return function value
          end if
       end if
    end do

    if(debug) write(*,*) 'MLH = ', mld_output(1), ', ML-avg theta_v = ', mld_output(2)

  end function MixedLayerDepth