program TestTouzePeiffer implicit none integer :: Nz, k, ierr, jj real, dimension(:), allocatable :: rho, thetav, zi, z, dz real :: delta_z, thetav0, thetav1, zstep, zsurface, zml, zml2 logical, parameter :: debug = .false. thetav0 = 290. thetav1 = 290.5 zstep = 400. ! step transition in thetav across this level zsurface = 0. ! ideally, this should correspond to a value of zi(k) do jj = 1,10 ! loop over different vertical resolutions. ! Each puts the top of the profile around 800m select case (jj) case (1) Nz = 20 delta_z = 40. case (2) Nz = 24 delta_z = 34. case (3) Nz = 40 delta_z = 20. case (4) Nz = 48 delta_z = 17. case (5) Nz = 80 delta_z = 10. case (6) Nz = 160 delta_z = 5. case (7) Nz = 320 delta_z = 2.5 case (8) Nz = 640 delta_z = 1.25 case (9) Nz = 1280 delta_z = 0.625 case (10) Nz = 2560 delta_z = 0.3125 end select ALLOCATE(rho(Nz),thetav(Nz),z(Nz),dz(Nz),zi(Nz+1),STAT=ierr) if(ierr.ne.0) STOP 'here' rho(:) = 1.2 thetav(:) = thetav0 zi(1) = zsurface do k = 1,Nz zi(k+1) = zi(k) + delta_z dz(k) = delta_z z(k) = 0.5*(zi(k)+zi(k+1)) ! midway between zi(k:k+1) ! modify thetav above zstep if(z(k).ge.zstep) thetav(k) = thetav1 if(debug) write(*,920) k, z(k), rho(k), thetav(k) 920 format('k/z/rho/thetav = ',I5,3F8.2) end do if(debug) write(*,*) zml = MixedLayerDepth(Nz,zi,z,rho,thetav) zml2 = MixedLayerDepthKazil(Nz,z,dz,rho,thetav) write(*,921) Nz, delta_z, zstep-zsurface, zml, zml2 921 format('Nz/deltaz/StepDepth/MLDepth(mod, Kazil) = ',I5,4F8.2) DEALLOCATE(rho,thetav,zi,z,dz,STAT=ierr) if(ierr.ne.0) STOP 'here2' end do contains real function MixedLayerDepth(Nz_in,zi_in,z_in,rho_in,thetav_in) ! Touze-Peiffer et al. (2022) mixed layer height algorithm (https://doi.org/10.1175/JAMC-D-21-0048.1). ! ! Note that following Touze-Peiffer et al., virtual potential temperature is based on water vapor mixing ! ratio r only, ignoring a contribution of liquid water: theta_v = theta*(1 + 0.61*r). ! ! Modified: jan.kazil@noaa.gov 2023-12-20 18:51:28 -07:00 (set z_min = 20) ! 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 depths, zi_in(1) is surface height real, intent(in) :: z_in(Nz_in+1) ! model level depths real, intent(in) :: rho_in(Nz_in), thetav_in(Nz_in) ! rho_in(1) is at height 0.5*SUM(zi_in(1: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) MixedLayerDepth = 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(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... MixedLayerDepth = z_in(k) - z_surface 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 ) MixedLayerDepth = z1 - f1*(z2-z1)/(f2-f1) - z_surface end if exit ! break out of loop, return function value end if end if end do if(debug) write(*,*) 'MLH = ', MixedLayerDepth end function MixedLayerDepth real function MixedLayerDepthKazil(z_n,z,dz,rho_z,theta_v_z) ! Touze-Peiffer et al. (2022) mixed layer height algorithm (https://doi.org/10.1175/JAMC-D-21-0048.1). ! ! Note that following Touze-Peiffer et al., virtual potential temperature is based on water vapor mixing ! ratio r only, ignoring a contribution of liquid water: theta_v = theta*(1 + 0.61*r). ! ! Modified: jan.kazil@noaa.gov 2023-12-20 18:51:28 -07:00 (set z_min = 20, cosmetic changes) ! Modified: Peter Blossey (convert to function), Dec 2023 ! Author: jan.kazil@noaa.gov 2023-07-26 15:09:13 -06:00 implicit none integer, intent(in) :: z_n integer :: t_i,x_i,y_i,z_i real :: a,b,z1,z2,f1,f2 real :: z_min,epsilon real :: integral,weight real :: thetav_mean_im1,thetav_mean_i,thetav_prime,thetav_mean_prime real :: integral_z(z_n),weight_z(z_n),ratio_z(z_n) ! ! Data (SI units) - not populated in this example ! real, intent(in) :: z(z_n) ! Mass level altitude real, intent(in) :: dz(z_n) ! Mass level thickness real, intent(in) :: rho_z(z_n) ! Air density real, intent(in) :: theta_v_z(z_n) ! Virtual potential temperature (DRY - see https://doi.org/10.1175/JAMC-D-21-0048.1) real :: mlh ! Mixed layer height (discrete) real :: mlh_con ! Mixed layer height (continuous) z_min = 20.0 ! m epsilon = 0.2 ! K ! ! Calculate the mixed layer depth (as in https://doi.org/10.1175/JAMC-D-21-0048.1) ! integral = 0.0 weight = 0.0 do z_i = 1, z_n if (z(z_i) < z_min) then cycle endif integral = integral + theta_v_z(z_i)*rho_z(z_i)*dz(z_i) weight = weight + rho_z(z_i)*dz(z_i) if (theta_v_z(z_i) >= integral/weight + epsilon) then mlh = z(z_i) exit endif enddo ! Calculate the mixed layer depth (following https://doi.org/10.1175/JAMC-D-21-0048.1), ! but with linear interpolation to make the mixed layer depth a continuous variable, ! as opposed to a discrete one (only defined on the vertical levels) ! Calculate the vertical integrals and their ratio in Eq. 1 (https://doi.org/10.1175/JAMC-D-21-0048.1) integral_z(:) = 0.0 weight_z(:) = 0.0 ratio_z(:) = 0.0 do z_i = 1, z_n if (z(z_i) < z_min) then cycle endif if (z_i == 1) then integral_z(z_i) = theta_v_z(z_i)*rho_z(z_i)*dz(z_i) weight_z(z_i) = rho_z(z_i)*dz(z_i) ratio_z(z_i) = integral_z(z_i)/weight_z(z_i) else integral_z(z_i) = integral_z(z_i-1) + theta_v_z(z_i)*rho_z(z_i)*dz(z_i) weight_z(z_i) = weight_z(z_i-1) + rho_z(z_i)*dz(z_i) ratio_z(z_i) = integral_z(z_i)/weight_z(z_i) endif enddo ! Using interpolation, determine the altitude at which Eq. 1 (https://doi.org/10.1175/JAMC-D-21-0048.1) holds do z_i = 1, z_n - 1 if (z(z_i) < z_min) then cycle endif if (theta_v_z(z_i+1) >= ratio_z(z_i+1) + epsilon) then z1 = z(z_i) z2 = z(z_i+1) f1 = theta_v_z(z_i ) - (ratio_z(z_i ) + epsilon) f2 = theta_v_z(z_i+1) - (ratio_z(z_i+1) + epsilon) a = (f2-f1)/(z2-z1) b = f1 - a*z1 mlh_con = -b/a exit endif enddo MixedLayerDepthKazil = mlh_con end function MixedLayerDepthKazil end program TestTouzePeiffer