subroutine compute_state( z, rho, p, T, theta, bvsq, & grav, R_gas, gamma_gas, i_type, i_unstable, N ) c------------------------------------------------------------------------------ c Given either the temperature or potential temperature profile, compute_state c integrates the hydrostatic balance to produce profiles of all other c thermodynamic quantities. The flag i_type controls whether temperature or c potential temperature is specified. The density at the lowest grid point c (rho(1)) is a required input. This routine should be compiled in double c precision via the -r8 flag or equivalent. c c Inputs: c z(1:N) - Altitude array. c rho(1) - Density at the lowest grid point. c T(1:N) - Temperature array. Required if i_type=0. c theta(1:N) - Potential temperature array. Required if i_type=1. c grav - Gravitational acceleration. c R_gas - Gas constant for air. c gamma_gas - Specific heat ratio for air. c i_type - 0 for specified temperature, c 1 for specified potential temperature. c i_unstable - 0 unstable atmosphere is converted to neutrally stable. c 1 unstable atmosphere is left alone. c N - Number of points in the profiles. c c Outputs: c rho(2:N) - Density profile. Only the first value (rho(1)) was an input. c p(1:N) - Pressure profile. c T(1:N) - Temperature array. Output if i_type=1, otherwise input. c theta(1:N) - Potential temperature array. Output if i_type=0, c otherwise input. c bvfsq(1:N) - Buoyancy frequency squared array. c------------------------------------------------------------------------------ implicit none real z(N), rho(N), p(N), T(N), theta(N), bvsq(N) real grav, R_gas, gamma_gas real r1, r2, cp, p_0, dtheta(500), & alpha, beta, del_T, del_z, del_theta, ratio, dzm, dzp, dz, & wtp, wtm, dthdzm, dthdzp, dthdz, theta_mid integer i, k, m, L, N, i_type, i_unstable, i_count, & k_start, k_stop, k_left, k_right character*1 stab_type(10) r1 = ( gamma_gas - 1.0 )/gamma_gas r2 = 1.0/r1 cp = R_gas*r2 c *** When i_type=0 you must pass in a valid entry for rho(1) as well c *** as the T(1:N) array. c *** When i_type=1 you must pass in a valid entry for rho(1) as well c *** as the theta(1:N) array. c *** In either case, we set theta(1) = T(1) if( i_type .eq. 0 ) then theta(1) = T(1) else T(1) = theta(1) end if p_0 = rho(1)*R_gas*T(1) p(1) = p_0 c *** i_type=0 for specified temperature, =1 for potential temperature if( i_type .eq. 0 ) then c *** Integrate the hydrostatic relation for specified temperature do k=2, N del_T = T(k) - T(k-1) del_z = z(k) - z(k-1) if( del_T .eq. 0.0 ) then p(k) = p(k-1)*exp(-grav/(R_gas*T(k))*del_z) else beta = del_T/del_z p(k) = p(k-1)*(T(k)/T(k-1))**(-grav/(beta*R_gas)) end if rho(k) = p(k)/(R_gas*T(k)) theta(k) = T(k)*(p(k)/p_0)**(-r1) end do end if c *** Remove unstable regions if the i_unstable flag is set i_count = 0 if( i_unstable .eq. 0 ) then do k=2, N dtheta(k) = theta(k) - theta(k-1) if( dtheta(k) .lt. 0.0 ) i_count = i_count + 1 end do if( i_count .eq. 0 ) goto 10 m = 1 if( dtheta(2) .gt. 0 ) then stab_type(1) = 'S' else stab_type(1) = 'U' end if do k=3, N if( dtheta(k) .gt. 0 .and. stab_type(m) .eq. 'U' ) then m = m + 1 stab_type(m) = 'S' end if if( dtheta(k) .lt. 0 .and. stab_type(m) .eq. 'S' ) then m = m + 1 stab_type(m) = 'U' end if end do c print *, (stab_type(i),i=1,m) c *** Fix an unstable region that starts at the ground. if( dtheta(2) .lt. 0.0 ) then do k=3, N if( dtheta(k) .gt. 0.0 ) goto 20 end do print *, 'Entire profile is untable' stop 20 k_stop = k - 1 do k=1, k_stop theta(k) = theta(k_stop) dtheta(k) = 0.0 end do end if c *** Loop to identify and repair any unstable egions on the profile c *** interior. do L=1, 10 do k=2, N if( dtheta(k) .lt. 0.0 ) goto 30 end do goto 10 30 k_start = k-1 do k=k_start+1, N if( dtheta(k) .ge. 0.0 ) goto 40 end do 40 k_stop = k-1 theta_mid = 0.5*( theta(k_start) + theta(k_stop) ) do k=k_start-1, 1, -1 if( theta(k) .le. theta_mid ) goto 50 end do 50 k_left = k+1 do k=k_stop+1, N if( theta(k) .ge. theta_mid ) goto 60 end do 60 k_right = k-1 do k=k_left, k_right theta(k) = theta_mid if( k .gt. 1 ) dtheta(k) = theta(k) - theta(k-1) end do if( k_right+1 .le. N ) then dtheta(k_right+1) = theta(k_right+1) - theta(k_right) end if end do end if 10 continue if( i_type .eq. 1 .or. i_count .gt. 0 ) then c *** Integrate the hydrostatic relation for specified potential temperature do k=2, N del_z = z(k) - z(k-1) del_theta = theta(k) - theta(k-1) if( del_theta .eq. 0.0 ) then T(k) = T(k-1) - grav/cp*del_z else alpha = del_theta/del_z ratio = theta(k)/theta(k-1) T(k) = ratio*T(k-1) - & grav*theta(k)/(cp*alpha)*alog(ratio) end if p(k) = p_0*(T(k)/theta(k))**r2 rho(k) = p(k)/(R_gas*T(k)) end do end if c *** Compute the buoyancy frequency squared do k=2, N-1 dzm = z(k ) - z(k-1) dzp = z(k+1) - z(k ) dz = dzm + dzp wtp = dzm/dz wtm = 1.0 - wtp dthdzm = ( theta(k ) - theta(k-1) )/dzm dthdzp = ( theta(k+1) - theta(k ) )/dzp dthdz = wtm*dthdzm + wtp*dthdzp bvsq(k) = grav/theta(k)*dthdz if( k .eq. 2 ) then dthdz = 2.0*dthdzm - dthdz bvsq(k-1) = grav/theta(k-1)*dthdz end if if( k .eq. N-1 ) then dthdz = 2.0*dthdzp - dthdz bvsq(k+1) = grav/theta(k+1)*dthdz end if end do c if( i_count .ge. 2 ) then c print *, 'Info: ', (stab_type(i),i=1,m), i_count, c & 'unstable atmosphere points converted to neutral' c do k=1, N c write(99,99) z(k), theta(k), bvsq(k) c99 format(1p,3e12.4) c end do c stop c end if return end