;******************************************************************************************* ; Liquid water and ice retrievals ; This retrieval requires the arm mwr. ;******************************************************************************************* ;condensed_water_parameterization_cirrus_snow3_sb, cld_base, ice_micro_source, $ ;liq_micro_source, locale, cld_top,lwp_in,j pro condensed_water_parameterization_cirrus_snow3_sb, cld_base, ice_micro_source, $ liq_micro_source, locale, cld_top,lwp,time_index,site ;******************************************************************************************* ; Commons ;******************************************************************************************* common cond_input, lwp_modis,ze_input,temp,vdop,vap,z,num_layers,optical_depth common cond_output, iwc, lwc, frac common cond_output2, re_liq, re_ice common cond_output3, cirrus_extinction common cond_output4, stratus_extinction common microphys_params, N0_output,N0_output_error,D0_output,D0_output_error,$ am_output,aa_output,bm_output,ba_output,alpha_output,beta_output common pwr_law_coeff, az, av, am, bz, bv, bm common pwr_law_coeff_bi2_st, am_alg, aa_alg, bm_alg, ba_alg, alpha_alg, beta_alg common Sb, frac_err_az,frac_err_bz,frac_err_av,frac_err_bv,frac_err_am,frac_err_bm, $ r_am_bm,r_am_av,r_am_bv,r_am_az,r_am_bz,r_bm_av,r_bm_bv,r_bm_az,r_bm_bz,r_av_bv,$ r_bv_az,r_bv_bz,r_az_bz, r_av_az, r_av_bz common Sa, r_iwc_dmass,frac_var_iwc,frac_var_dmass common Se, frac_err_z,frac_err_vdq,r_Z_vdq common seedit, seed ;******************************************************************************************* ; Set some variables ;******************************************************************************************* ; missing data flag missing=-8888 ; clear sky flag clearsky=-9999 ; Do you want print statements? print_flag=0 ; If the mwr lwp is less than 20g/m2 then the value is considered unreliable. ; Another value will be used. mwr_threshold=20. ; Find the number of elements in the height array numheights=n_elements(z) ; Calculate height interval (used to be hardcoded as 90) ddz=z[20]-z[19] ;****************************************************************************************** ; Create necessary arrays ;****************************************************************************************** ccm3_ice_frac=make_array(/float,numheights,value=0) ccm3_liquid_profile=make_array(/float,numheights,value=0) liquid_wt=fltarr(numheights) lwc=fltarr(numheights) iwc=fltarr(numheights) re_liq=fltarr(numheights) re_ice=fltarr(numheights) cirrus_extinction=fltarr(numheights) stratus_extinction=fltarr(numheights) N0_output=fltarr(numheights) N0_output_error=fltarr(numheights) D0_output=fltarr(numheights) D0_output_error=fltarr(numheights) am_output=fltarr(numheights) aa_output=fltarr(numheights) bm_output=fltarr(numheights) ba_output=fltarr(numheights) alpha_output=fltarr(numheights) beta_output=fltarr(numheights) ;******************************************************************************************* ; If rain is being fasely determined as cloud base, move cloud base above the rain. ;******************************************************************************************* ; fast forward height up to cloud base i=0 while z[i] le cld_base do i=i+1 ; Check Ze of cloud base. If it is gt 20 then fast forward and put cloud base above ; the 20 Ze heights. This is considered rain. if 10.*alog10(ze_input[i]) gt 20. then begin while 10.*alog10(ze_input[i]) gt 20. and i lt n_elements(ze_input)-1 do i=i+1 cld_base=z[i] endif ;******************************************************************************************* ; Store the temp at cloud base ;******************************************************************************************* cld_base_temp=temp[i] ;******************************************************************************************* ; Put all the reflectivity under cloud base to clear sky value. ; Only remove these reflectivity if cld_base is above freezing. ; (only remove rain, leave in the snow) ;******************************************************************************************* ;if cld_base gt 0. then begin if cld_base gt 0. and cld_base_temp gt 273. then begin k=0 while z[k] lt cld_base do begin ze_input[k]=clearsky liq_micro_source[k]=7 ice_micro_source[k]=7 k=k+1 endwhile endif ze=ze_input ;******************************************************************************************* ; Liquid Retrieval ;******************************************************************************************* ; Calculate the fraction of the profile that should be liquid and ice according ; to the ccm3 model equations. The fraction of ice depends on temperature. ; The amount of liquid depends on the mwr vap and height. ; Calculate the lwp according to the ccm3. ccm3_lwp=0. for j=0,numheights-1 do begin if ze[j] gt missing then begin ; gives an ice and water percentage at this height based on temp only ; -10C=0% ice (0), -30C=100% ice (1) ccm3_ice_frac[j]=ccm3_cloud_iceliquid_fraction_param_function(temp[j]) ; gives the amount of liquid water at this height, function of height,vap. ccm3_liquid_profile[j]=(1.-ccm3_ice_frac[j])*ccm3_liquid_water_param_function(z[j], vap) ; sums the ccm3 lwp in the column ccm3_lwp=ccm3_lwp+(ccm3_liquid_profile[j]*ddz) endif endfor ; What fraction of the lwp occurs at temps colder than freezing based on the CCM3 profile? ; Calculate the amount of cold water and warm water in the ccm3 profile ; This depends on temperature only and must be above cloud base. warm_water=0. & total_sqrtZ_warm=0. & ccm3_cold_water=0. & total_3rtZ_warm=0. for j=0,numheights-1 do begin ; Sum warm liquid water if temp[j] gt 273. and ze[j] gt missing and z[j] ge cld_base then begin warm_water=warm_water+ccm3_liquid_profile[j] total_sqrtZ_warm=total_sqrtZ_warm+(sqrt(ze[j])) total_3rtZ_warm=total_3rtZ_warm+((ze[j])^(1./3.)) endif ; Sum cold liquid water if temp[j] le 273. and ze[j] gt missing and temp[j] gt 233. and z[j] ge cld_base then begin ccm3_cold_water=ccm3_cold_water+ccm3_liquid_profile[j] endif endfor ; If there is any liquid water in this profile calculate ; lwc, liq_micro_source, re_liq warm_water_path=0. cold_water_path=0. if ccm3_cold_water gt 0. or warm_water gt 0. then begin ; If the lwp is less than 20g/m2 then use the ccm3 lwp value if lwp lt mwr_threshold then lwp=ccm3_lwp warm_water_path=lwp*(warm_water/(total(ccm3_liquid_profile))) cold_water_path=lwp-warm_water_path ; so now distribute the warm water path vertically based on square root of Z for j=0,numheights-1 do begin ; Calculate liquid water content (LWC) for warm water if temp[j] gt 273. and ze[j] gt -8888. then begin lwc[j]=((sqrt(ze[j]))/total_sqrtZ_warm)*(warm_water_path/ddz) liq_micro_source[j]=5 if lwc[j] lt 1.e-4 then lwc[j]=0.0 ; Calculate liquid water content (LWC) for cold water endif else if ze[j] gt -8888. and temp[j] gt 243. then begin lwc[j]=(cold_water_path*(ccm3_liquid_profile[j]/ccm3_cold_water))/ddz liq_micro_source[j]=6 if lwc[j] lt 1.e-4 then lwc[j]=0.0 endif ; Calculate liquid effective radius if 10.*alog10(ze[j]) lt 0.0 and ze[j] gt -8888. and lwc[j] ge 1.e-4 then begin re_liq[j]=19.5*exp(0.0384*(10.*alog10(ze[j]))) re_liq[j]=re_liq[j]+(0.20*re_liq[j]) endif else begin if ze[j] gt -8888. and lwc[j] ge 1.e-4 then re_liq[j]=15. endelse ;if re_liq[j] gt 15 then begin ; print,time_index,ze[j],lwc[j],re_liq[j],'re_liq' ;endif endfor endif ; if there is any liquid then calculate lwp and effective radius ;****************************************************************************************** ; Plot the liquid retrieval for a certain profile ;****************************************************************************************** ;if time_index eq 2 then begin ;if time_index ge 195 and time_index le 234 then begin ; plot_inputs,time_index,ccm3_ice_frac,ccm3_liquid_profile,ze,temp,z,$ ; lwc,re_liq,vap,lwp,ccm3_lwp,cold_water_path,warm_water_path ;endif ;******************************************************************************************* ; Ice Retrieval ; If there is reflectivity and temp is lt 0C then ice retrieval is done. ;******************************************************************************************* ; ; now determine the ice water content from the regression data. Assume all water warmer ; than -15c and the regression below applies to all ice at colder temps ; Calculate iwc, re_ice, and ice_micro_source coef=[-8.434e-009,6.721e-008,2.409e-010] bz=-2.4 & ba=1.0 & bm=1.8 & bi=0.6 ; Loop through the heights for j=0,numheights-1 do begin ; If temp is less than 0C and there is reflectivity ;if 1 eq 0 then begin ;turns off this ice retrieval if ze[j] gt -8888. and temp[j] lt 273. then begin ; If temp is 0 to -5C, might need to change ze ;if temp[j] gt 268. and 10.*alog10(ze[j]) gt 10. then begin ; zzz=10.^(10./10.) ;endif else begin zzz=ze[j] ;endelse ; Adjust vdop if negative if vdop[j] lt 0. then vdop[j]=0.5 ;big ice or snow v_temp=vdop[j]*(((zzz)))^((ba/((6.+bz)))) t_temp=611.*exp(((2.5e6)/461.5)*((1./273.)-(1./temp[j]))) ;vapor pressure fofvz=alog10((0.251)*((zzz*1.e-12)/0.00164)*((5135./(vdop[j]*100.))^1.18)) iwc[j]=10.^(((-2.22)+(0.638*fofvz)+((-7.18e-5)*t_temp))) ice_micro_source[j]=1 iwc[j]=iwc[j]*1.e6 ; converts to g/m3 ; Ice particle size depends on height re_ice[j]=ccm3_cloud_ice_re_param_function(exp(-z[j]/5500.)) re_ice[j]=re_ice[j]+(re_ice[j]*0.10) ; a test to decrease the optical depth 6/24/05 ; If there is no reflectivity or temp is gt then 0C then no ice here endif else begin iwc[j]=0.0 re_ice[j]=0. endelse endfor ;****************************************************************************************** ; Plot the ice retrieval for a certain profile ;****************************************************************************************** ;if time_index eq 2 then begin ;if time_index ge 2 and time_index le 4 then begin ; plot_ice_retrieval,time_index,ccm3_ice_frac,ccm3_liquid_profile,ze,temp,z,lwc,re_liq,vap,$ ; iwc,re_ice,pos,dmin_re,dmax_re,dmin_wc,dmax_wc,plot_top_height ;endif ;*************************************************************************************** ; Find the cloud layers in the profile and determine properties ;*************************************************************************************** layer_determination,missing,clearsky,ze,z,temp,j,cirrus,numlayers,layer_base_index,layer_top_index,$ max_dbz_temp,base_temp,layer_top_temp,base_height,top_height ; at this point we have retrievals of ice and liquid. It may be the case that the ; layer is really all ice and the liquid can be neglected or assumed not to exist. ; So, do the ice retrieval with the available information if the layer base is colder than ; freezing, the input lwp is less than the threshold ;************************************************************************************** ; Implement ciret2. ; Loop through the cloud layers and perform arm cirrus retrieval if the layer ; meets cirrus conditions. ;************************************************************************************** for i=0,numlayers-1 do begin ; If these conditions are met, do arm cirrus retrieval if max_dbz_temp[i] le 273. and layer_top_temp[i] le 253. and $ layer_base_index[i] gt 0 then begin ; we found a cirrus layer if n_elements(rho_air) eq 0 then rho_air=-9999. for kk=layer_base_index[i], layer_top_index[i] do begin if layer_top_index[i] gt layer_base_index[i] then begin top_dist=(top_height[i]-z[kk])/(top_height[i]-base_height[i]) endif else begin top_dist=ddz endelse aa=0.32783559+(0.010212126*(10.*alog10(ze[kk])))+(0.0011355369*temp[kk])+$ ((-0.037706691)*top_dist)+((-0.0023235938)*(vdop[kk]*100.)) ba=1.3906124+((-0.021529097)*(10.*alog10(ze[kk])))+((-0.0020046848)*temp[kk])+$ ((0.071951995)*top_dist)+((0.0043448790)*(vdop[kk]*100.)) if aa lt 0. or ba gt 2.0 then begin aa=0.3 & ba=1.9 endif if ba lt 1.1 then begin aa=0.08 & ba=1.1 endif am=0.0031 & bm=2.26 if rho_air le 0. then begin if float(temp[kk]) ge 253. then pressure=50000. if float(temp[kk]) lt 253. then pressure=40000. if float(temp[kk]) lt 243. then pressure=30000. if float(temp[kk]) lt 233. then pressure=25000. if float(temp[kk]) lt 223. then pressure=20000. rho_air=1.e-3*(pressure/(287.04*float(temp[kk]))) ;print,'calculated rho_air',i,'layer of cirrus' endif az=(0.19)*((6./(3.14159*0.92))^2)*((am)^2) bz=(2.*bm)-6. ab=(0.04394+0.06049)/2. bb=(0.970+0.831)/2. dyn_visc=(1.7e-5)*(10.) ; the 1000/100 converts from kg/m*s to g/cm*s kin_visc=dyn_visc/rho_air g=981. bv=(bb*(bm+2.-ba))-1. av=ab*kin_visc*((2.*am*g/(rho_air*(kin_visc^2)*aa))^bb) frac_err_z=0.58 & frac_err_vdq=0.5 ; these are estimates from known error in the observations frac_var_iwc=2.0 & frac_var_dmass=1.0 frac_var_bv=0.5 & frac_var_bm=0.5 & frac_var_av=5. & frac_var_am=5. & frac_var_az=5. & frac_var_bz=0.5 frac_err_az=0.3 & frac_err_bz=0.3 & frac_err_av=0.3 & frac_err_bv=0.3 & frac_err_am=0.3 & frac_err_bm=0.3 r_am_bm=0.01 & r_am_av=0.01 & r_am_bv=0.01 & r_am_az=0.01 & r_am_bz=0.01 & r_bm_av=0.01 & r_bm_bv=0.01 & r_bm_az=0.01 r_bm_bz=0.01 & r_av_bv=0.01 & r_bv_az=0.01 & r_bv_bz=0.01 & r_az_bz=0.01 r_av_az=0.01 & r_av_bz=0.01 & r_iwc_dmass=0.01 & frac_var_iwc=0.01 & frac_var_dmass=0.01 & r_Z_vdq=0.01 ciret2_optimal_iwc_Dmass_ap_norm_loglininp, rho_air*1000., 10.*alog10(ze[kk]), vdop[kk], N0_ret, lambda_ret, iwc_ret, conc_ret, $ count_median_length_ret,mass_mean_length_ret, v_mass_mean_ret, N0_ret_err, lambda_ret_err, $ iwc_ret_err, conc_ret_err, $ count_median_length_ret_err,mass_mean_length_ret_err, v_mass_mean_ret_err, $ iwc_exact_ret,dmass_exact_ret, $ dge_ret, dge_ret_err, effective_radius_ret, effective_radius_ret_err, ext_coeff_ret, ext_coeff_ret_err,param_fname iwc_LI=protat_cirrus_regression((10.*alog10(ze[kk])), temp[kk], locale,'k') if abs(iwc_LI-(iwc_ret))/iwc_LI gt 2. then begin N0_ret_err=-9999. dge_ret=-9999. iwc_ret=-9999 iwc[kk]=iwc_LI re_ice[kk]=40. ice_micro_source[kk]=3 endif ; get extinction and write to output arrays if N0_ret_err gt -9999. then begin ; check to see if the cirrus is not tropical if locale ne 'trop' then begin cirrus_extinction[kk]=ext_coeff_expo_zv_bvbzbm(ze[kk]*1.d-12, vdop[kk]*100., double(bm), double(bv),$ double(bz), double(am), double(av), double(az), double(rho_air)) endif else begin cirrus_extinction[kk]=0.0 endelse if dge_ret gt 0.0 then begin iwc[kk]=iwc_ret re_ice[kk]=effective_radius_ret ice_micro_source[kk]=4 N0_output[kk]=N0_ret N0_output_error[kk]=N0_ret_err D0_output[kk]=1./lambda_ret D0_output_error[kk]=1./lambda_ret_err am_output[kk]=am aa_output[kk]=aa bm_output[kk]=bm ba_output[kk]=ba alpha_output[kk]=0. beta_output[kk]=0. endif endif endfor ; layer base index to layer top index loop endif ;If this is a cirrus layer endfor ;end of loop through cloud layers ;****************************************************************** ; In cirrus layers, set lwc, liq_micro_source, and re_liq to 0 ; because there should be no liquid water in cirrus ; lwc (lwc_outp) liq_micro_source re_liq (rel_outp) ;****************************************************************** result=where(cirrus eq 1,count) if count gt 0 and total(lwc) gt 0 then begin for l=0,count-1 do begin top=layer_top_index[result[l]] base=layer_base_index[result[l]] if total(lwc[base:top]) gt 0 then begin lwc[base:top]=0 liq_micro_source[base:top]=8 re_liq[base:top]=0 endif endfor endif ;********************************************************** ; Look for a temp inversion below 4km where the temp ; range is 260-273 K *** don't do this for nsa *** ;********************************************************** if site ne 'nsa' then begin ; First find height, temp and mwr condition result=where(z le 4000 and temp ge 260 and temp le 273 and lwp gt mwr_threshold,count) ; If those conditions are satisfied, see if there are cloud layers contained in 0-4km height range. if count gt 0 then begin result2=where(base_height ge 0 and top_height le 4000,count2) ; If there are clouds, look for inversion if count2 gt 0 then begin ; Now look for a temperature inversion num_inv=0 for i=0,count-2 do begin ;print,count,i,z[result[i]],temp[result[i]] if temp[result[i]] lt temp[result[i+1]] then begin print,'inversion' num_inv=num_inv+1 endif endfor ; Found the temp inversion, get rid of the ice in these clouds if num_inv gt 0 then begin print,time_index,num_inv for j=0,count2-1 do begin result3=where(z ge base_height[j] and z le top_height[j],count3) if count3 gt 0 then begin iwc[result3]=0 ice_micro_source[result3]=8 re_ice[result3]=0 endif endfor endif endif endif endif ;don't do this for nsa ;*********************************************************** ; NSA ; If cloud base is below freezing, then only ice should be ; below the cloud base. Set the liquid values to 0 ;*********************************************************** if cld_base gt 0. and cld_base_temp lt 273. then begin result=where(z lt cld_base and ze gt missing,count) if count gt 0 then begin ;print,'j',time_index ;print,'lwc',lwc[result] ;print,'lms',liq_micro_source[result] ;print,'rel',re_liq[result] lwc[result]=0 liq_micro_source[result]=2 re_liq[result]=0 endif endif ;if time_index eq 219 then stop ;********************************************************** ; plot and see if these values were changed ;********************************************************** do_this='no' if do_this eq 'yes' then begin alldata=[dmin_wc,dmax_wc] dmax=max(alldata) dmin=min(alldata) plot,lwc,z,color=0,position=pos[1,*],xstyle=5,ystyle=5,$ /noerase,xrange=[dmin,dmax],/nodata,yrange=[z[0],plot_top_height] oplot,lwc,z,color=230,linestyle=3,thick=2 oplot,iwc,z,color=30,linestyle=3,thick=2 alldata=[dmin_re,dmax_re] dmax=max(alldata) dmin=min(alldata) plot,re_liq,z,color=0,position=pos[2,*],xstyle=5,ystyle=5,$ /noerase,xrange=[dmin,dmax],/nodata,yrange=[z[0],plot_top_height] oplot,re_liq,z,color=230,linestyle=3,thick=2 oplot,re_ice,z,color=30,linestyle=3,thick=2 write_gif,'/Users/u0079358/Documents/fifth_cloudret/flux_closureL/ice_retrieval/ice_retrieval.'+string(time_index,format='(I04)')+'.gif',tvrd() endif return end