C NCLFORTSTART subroutine sndganalncl(nlevs,t,p,r,u,v,ctemp) implicit none C ; input integer nlevs real t(nlevs), p(nlevs),r(nlevs),u(nlevs),v(nlevs) C ; returned real ctemp C NCLEND call sndg_anal(nlevs,t,p,r,u,v,ctemp) return end c---------------------------- subroutine sndg_anal(nlevs,t,p,r,u,v,conv_temp) c c..Subroutine to analyze stability, positive area, level info from a sounding c c..Local variables KLCL, KCCL, KLFC, KEL refer to the vertical index above c..the levels of LCL, CCL, LFC, and EL respectively. In other words, the c..exact location of the LCL is between levels KLCL and KLCL-1 in the data. c parameter (max_levs=100, zero=273.15) parameter (Rd=287.05, Cp=1004.0, G=9.8, XMS2KTS=1.94) dimension t(*),p(*),r(*),u(*),v(*),the(max_levs) character*128 label logical is_lfc c c..Initialize some stuff c RCP = Rd/Cp CPR = Cp/Rd area_pos = 0. area_pos2 = 0. area_neg = 0. area_neg2 = 0. is_lfc = .false. prec_w = 0. plfc = 0. tlfc = 0. pel = 0. tel = 0. cape = 0. sreh = 0. vgpi = 0. alifted = 0. akindex = 0. totals = 0. sweat = 0. c c..First lets find all the vertical indicies of important levels c..to make things easier later. (get array index corresponding to ??? hPa) c k850 = 1 k700 = 2 k500 = 3 do k=1,nlevs if(p(k).ge.85000.0) k850=k if(p(k).ge.70000.0) k700=k if(p(k).ge.50000.0) k500=k if(p(k).ge.30000.0) k300=k enddo k850 = max(1,k850) k700 = max(2,k700) k300 = min(k300,nlevs) c C+---+------------------------------------c c..Compute Dewpoint, LCL pres, Theta, Theta-e, Wet Bulb Theta for each level c c print11,' Temp | Pres | W | Th | Th-e | Td | RH |' c + ,' T_LCL|P_LCL | Th-Wet ' 11 format(a,a) hwbz = -999. do k=1,nlevs theta = t(k) * (100000.0/p(k))**RCP dew_t = t_dew(p(k),r(k)) tlcl = t_lcl(t(k),dew_t) plcl = 100000.0 * (tlcl/theta)**CPR c write(6,*) 'k = ',k,' p = ',p(k),' t = ',t(k),' r = ',r(k), c & ' tlcl = ',tlcl the(k)=theta_e(p(k),t(k),r(k),tlcl) rh = e_sub_s(dew_t) / e_sub_s(t(k)) thwet = theta_wetb(the(k)) twet = satlft (thwet, p(k)) if (hwbz .lt. 0. .and. twet .le. 0.) then hwbz = std_atmos(p(k)) if (k .eq. 1) hwbz = 0. endif cc print21, t(k)-zero,p(k)*0.01,r(k)*1000.,theta-zero,the(k)-zero cc + ,dew_t-zero,rh*100.,tlcl-zero,plcl*0.01,thwet-zero 21 format(f6.1,1x,f6.1,1x,f5.2,1x,f6.2,1x,f6.2,1x + ,f6.2,1x,f5.1,1x,f6.2,1x,f6.1,1x,f6.1) enddo c print* c C+---+------------------------------------c c..Compute precipitable water (meters) in column c prec_w=prec_water(p(1),r(1),nlevs) c C+---+------------------------------------c c..Compute mixed layer info using lowest 110 hPa data c tot_pres = 0. tot_mixr = 0. tot_temp = 0. sum_delta = 0. avg_temp = t(1) avg_mixr = r(1) avg_pres = p(1) do k=2,nlevs if(p(k).ge.p(1)-11000.0) then delta_p = p(k-1) - p(k) sum_delta = sum_delta + delta_p tot_pres = tot_pres + delta_p*(p(k)+p(k-1))*0.5 tot_mixr = tot_mixr + delta_p*(r(k)+r(k-1))*0.5 tot_temp = tot_temp + delta_p*(t(K)+t(k-1))*0.5 endif enddo if (sum_delta .gt. 0.0) then avg_temp = tot_temp / sum_delta avg_mixr = tot_mixr / sum_delta avg_pres = tot_pres / sum_delta endif avg_th = avg_temp * (100000.0/avg_pres)**RCP avg_t_dew = t_dew(avg_pres,avg_mixr) c C+---+------------------------------------c c..Compute Lifted Index (C) using SELS method based on 100hPa thick c..mixed layer temp + 2 degrees C c Abandoned the +2 for consistency... JFB c c avg_temp2 = avg_temp + 2. avg_temp2 = avg_temp avg_th2 = avg_temp2 * (100000.0/avg_pres)**RCP tlcl2 = t_lcl(avg_temp2,avg_t_dew) plcl2 = 100000.0 * (tlcl2/avg_th2)**CPR thelcl2 = theta_e(plcl2,tlcl2,avg_mixr,tlcl2) parcel500 = compT_fr_The(thelcl2,50000.0,iup) t500 = t(k500) + ((t(k500+1)-t(k500))/(p(k500+1)-p(k500))) & * (50000.-p(k500)) alifted = t500 - parcel500 c c compute showalter index c if (k850 .eq. 1) then t850 = t(k850) r850 = r(k850) else t850 = t(k850) + ((t(k850+1)-t(k850))/(p(k850+1)-p(k850))) & * (85000.-p(k850)) r850 = r(k850) + ((r(k850+1)-r(k850))/(p(k850+1)-p(k850))) & * (85000.-p(k850)) endif th850 = t850 * (100000.0/85000.0)**RCP d850 = t_dew(85000.,r850) tlcl8 = t_lcl(t850,d850) plcl8 = 100000.0 * (tlcl8/th850)**CPR thelcl8 = theta_e(plcl8,tlcl8,r850,tlcl8) pcl500 = compT_fr_The(thelcl8,50000.0,iup) swi = t500 - pcl500 c C+---+------------------------------------c c..Find LCL, CCL, EL, and LFC from the mixed layer mixing data c..These findings are exact. One assumption hurting the c..computations: theta-w is assumed to vary linearly between c..two points on the environment temperature curve. This is c..used in y=mx + b fit to interpolate any temperature between c..the two endpoints. It is a reasonable assumption with c..small delta-pressure spacing and in the lower levels but c..gets more unreasonable as the spacing increases and c..go higher. Still, it beats the alternative - just assigning c..the nearest level info. c c c..First LCL - Lifting Condensation Level and save the c..theta-e and theta-w from here cause we can check positive c..and negative areas and LFC, EL based on these. c tlcl = t_lcl(avg_temp,avg_t_dew) plcl = 100000.0 * (tlcl/avg_th)**CPR do k =1,nlevs if(p(k).le.plcl) goto 100 enddo 100 continue klcl = min(k,nlevs) thelcl = theta_e(plcl,tlcl,avg_mixr,tlcl) thwlcl = theta_wetb(thelcl) c c..Second, find CCL - Convective Condensation Level and c..thus, convective temperature. Start at 500hPa and go c..downward along the temperature curve until it crosses c..the mixing ratio line that goes thru the LCL. Do so by c..finding the first time the mixing ratio of the environment c..air (assuming saturation) exceeds the mixing ratio of c..the mixed layer parcel (the same as the one going thru c..the LCL). c do k=k500,2,-1 tmixr = r_sub_s(p(k),t(k)) if(tmixr .ge. avg_mixr) goto 110 enddo 110 continue kccl = min(k+1,nlevs) tmixr = r_sub_s(p(kccl),t(kccl)) if(abs(tmixr - r_sub_s(p(kccl-1),t(kccl-1))) .gt. 1.e-5 ) then slope = (tmixr-avg_mixr) / (tmixr-r_sub_s(p(kccl-1),t(kccl-1))) if (slope .lt. 0.) slope = 0. else slope = 0.5 endif tccl = t(kccl) - slope*(t(kccl)-t(kccl-1)) pccl = p(kccl) - slope*(p(kccl)-p(kccl-1)) thccl = tccl * (100000.0/pccl)**RCP conv_temp = min(323.15,thccl * (p(1)/100000.0)**RCP) return end function theta_wetb(thetae_K) real*8 c(0:6), d(0:6) data c/-1.00922292e-10, -1.47945344e-8, -1.7303757e-6 + ,-0.00012709, 1.15849867e-6, -3.518296861e-9 + ,3.5741522e-12/ data d/0.00000000, -3.5223513e-10, -5.7250807e-8 + ,-5.83975422e-6, 4.72445163e-8, -1.13402845e-10 + ,8.729580402e-14/ x=min(475.0,thetae_K) if( x .le. 335.5 ) then answer = c(0)+x*(c(1)+x*(c(2)+x*(c(3)+x*(c(4)+x*(c(5)+ + x*c(6) ))))) else answer = d(0)+x*(d(1)+x*(d(2)+x*(d(3)+x*(d(4)+x*(d(5)+ + x*d(6) ))))) endif theta_wetb = answer + 273.15 return end function t_dew(pres_Pa,w_non) p = pres_Pa RR=w_non+1e-8 ES=P*RR/(.622+RR) ESLN=LOG(ES) T_Dew=(35.86*ESLN-4947.2325)/(ESLN-23.6837) return end function t_lcl(temp_K,tdew_K) tt = temp_K tttd= tdew_K denom= ( 1.0/(tttd-56.0) ) + (log(tt/tttd)/800.) t_lcl = ( 1.0 / denom ) + 56.0 return end function r_sub_s(pres_Pa,temp_K) c.. c.. this calls function e_sub_s which computes saturation c.. vapor pressure (Pa) and converts to sat. mixing ratio (kg/kg) c.. pres_Pa - pressure (pa) c.. temp_K - temperature (k) c.. es=e_sub_s(temp_K) r_sub_s=.622*es/(pres_Pa-es) c return end c c+---+-----------------------------------------------------------------+ function e_sub_s(temp_K) c.. c.. compute saturation vapor pressure (Pa) over liquid with c.. polynomial fit of goff-gratch (1946) formulation. (walko, 1991) c.. dimension c(0:8) data c/610.5851,44.40316,1.430341,.2641412e-1,.2995057e-3 + ,.2031998e-5,.6936113e-8,.2564861e-11,-.3704404e-13/ x=max(-80.,temp_K-273.16) e_sub_s = c(0)+x*(c(1)+x*(c(2)+x*(c(3)+x*(c(4)+x*(c(5) + +x*(c(6)+x*(c(7)+x*c(8)))))))) return end function std_atmos(pres_Pa) cc.. cc.. pres_Pa = Pressure in Pascals cc.. standard atmos height in meters is returned for given p cc.. c pr = pres_Pa*0.01 c height = 44307.692 * (1.0 - (pr/1013.25)**0.190) c c std_atmos=height c c return end function compT_fr_The(thelcl_K,pres_Pa,iup) c.. c.. pres_Pa = Pressure in Pascals c.. thelcl = Theta-e at LCL (units in Kelvin) c.. c.. Temperature (K) is returned given Theta-e at LCL c.. and a pressure. This describes a moist-adiabat. c.. This temperature is the parcel temp at level Pres c.. along moist adiabat described by theta-e. c.. guess= (thelcl_K - 0.5 * ( max(thelcl_K-270., 0.))**1.05) + * (pres_Pa/100000.0)**.2 epsilon=0.01 do iter=1,100 w1 = r_sub_s(pres_Pa,guess) w2 = r_sub_s(pres_Pa,guess+1.) tenu = theta_e(pres_Pa,guess,w1,guess) tenup = theta_e(pres_Pa,guess+1,w2,guess+1.) cor = (thelcl_K - tenu) / (tenup - tenu) guess = guess + cor if( (cor.lt.epsilon) .and. (-cor.lt.epsilon) ) then compT_fr_The=guess return endif enddo write(iup,*) & 'Convergence not reached in compT_fr_The, continuing.' thwlcl_K=theta_wetb(thelcl_K) compT_fr_The = thwlcl_K*((pres_Pa/100000.0)**0.286) return end function prec_water(pres_Pa, w_non, nlevels) c.. c.. pres_Pa = Pressure in Pascals c.. w_non = mixing ratio (non-dimensional = kg/kg) c.. returned precipitable water value in meters only below 300hPa c.. dimension pres_Pa(*), w_non(*) parameter (RHO_W=1000.0, G=9.8) sum=0. do k=nlevels,2,-1 if(pres_Pa(k).gt.30000.0) goto 10 enddo 10 continue do kk=k,2,-1 sum = sum + ( (w_non(kk)+w_non(kk-1))*0.5) + * abs(pres_Pa(kk)-pres_Pa(kk-1)) enddo prec_water = sum / (G*RHO_W) return end function theta_e(pres_Pa,temp_K,w_non,tlcl_K) c.. c.. The following code was based on Bolton (1980) eqn #43 c.. and claims to have 0.3 K maximum error within -35 < T < 35 C c.. pres_Pa = Pressure in Pascals c.. temp_K = Temperature in Kelvin c.. w_non = mixing ratio (non-dimensional = kg/kg) c.. tlcl_K = Temperature at Lifting Condensation Level (K) c.. pp = pres_Pa tt = temp_K rr = w_non + 1.e-8 c write(6,*) 'rr = ',rr,' w_non = ',w_non tlc = tlcl_K c power=(0.2854*(1.0 - (0.28*rr) )) xx = tt * (100000.0/pp)**power p1 = (3.376/tlc) - 0.00254 p2 = (rr*1000.0) * (1.0 + 0.81*rr) c write(6,*) 'p1 = ',p1,' p2 = ',p2 theta_e = xx * exp(p1*p2) return end function satlft (thw1, p1) parameter (cta = 273.16, akap = 0.28541) thw = thw1 - 273.15 p = p1 * .01 if (p .ne. 1000.) go to 5 satlft = thw return 5 continue pwrp = (p / 1000.) ** akap tone = (thw + cta) * pwrp - cta eone = wobf(tone) - wobf(thw) rate = 1. go to 15 10 continue rate = (ttwo - tone) / (etwo - eone) tone = ttwo eone = etwo 15 continue ttwo = tone - eone * rate pt = (ttwo + cta) / pwrp - cta etwo = pt + wobf(ttwo) - wobf(pt) - thw dlt = etwo * rate if (abs(dlt) .gt. 0.1) go to 10 satlft = ttwo - dlt return end c------------------------------------------------------------- function wobf (t) x = t - 20. if (x .gt. 0.) go to 10 pol = 1. + x * (-8.8416605e-03 & + x * ( 1.4714143e-04 + x * (-9.6719890e-07 & + x * (-3.2607217e-08 + x * (-3.8598073e-10))))) wobf = 15.130 / pol ** 4 return 10 continue pol = 1. + x * ( 3.6182989e-03 & + x * (-1.3603273e-05 + x * ( 4.9618922e-07 & + x * (-6.1059365e-09 + x * ( 3.9401551e-11 & + x * (-1.2588129e-13 + x * ( 1.6688280e-16))))))) wobf = 29.930 / pol ** 4 + 0.96 * x - 14.8 return end