c subroutine sndg_anal(t,p,r,u,v, & flminsou,frmaxsou,fbminsou,ftmaxsou,nlevs) 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