docs: def Akmaev_adjustment(theta, q, beta, n_k, theta_k, s_k, t_k):
      '''Single column only.'''
      L = q.size  # number of vertical levels
      # Akmaev step 1
      k = 1
      n_k[k-1] = 1
      theta_k[k-1] = theta[k-1]
      l = 2
      while True:
          # Akmaev step 2
          n = 1
          thistheta = theta[l-1]
          while True:
              # Akmaev step 3
              if theta_k[k-1] <= thistheta:
                  # Akmaev step 6
                  k += 1
                  break  # to step 7
              else:
                  if n <= 1:
                      s = q[l-1]
                      t = s*thistheta
                  # Akmaev step 4
                  if n_k[k-1] <= 1:
                  # lower adjacent level is not an earlier-formed neutral layer
                      s_k[k-1] = q[l-n-1]
                      t_k[k-1] = s_k[k-1] * theta_k[k-1]
                  # Akmaev step 5
                      #  join current and underlying layers
                  n += n_k[k-1]
                  s += s_k[k-1]
                  t += t_k[k-1]
                  s_k[k-1] = s
                  t_k[k-1] = t
                  thistheta = t/s
                  if k==1:
                      #  joint neutral layer is the first one
                      break  # to step 7
                  k -= 1
                  # back to step 3
          # Akmaev step 7
          if l == L:  # the scan is over
              break  # to step 8
          l += 1
          n_k[k-1] = n
          theta_k[k-1] = thistheta
          # back to step 2
  
      # update the potential temperatures
      while True:
          while True:
              # Akmaev step 8
              if n==1: # current model level was not included in any neutral layer
                  break  # to step 11
              while True:
                  #  Akmaev step 9
                  theta[l-1] = thistheta
                  if n==1:
                      break
                  #  Akmaev step 10
                  l -= 1
                  n -= 1
                  #  back to step 9
          # Akmaev step 11
          if k==1:
              break
          k -= 1
          l -= 1
          n = n_k[k-1]
          thistheta = theta_k[k-1]
          # back to step 8
      return theta