pro condensed_water_parameterization_cirrus_snow, cld_base

common cond_input, lwp,ze_input,temp,vdop,vap,z

common cond_output, iwc, lwc, frac
common cond_output2, re_liq, re_ice
common cond_output3, cirrus_extinction

common pwr_law_coeff, az, av, am, bz, bv, bm
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

ccm3_liquid_profile=fltarr(n_elements(z)) & liquid_wt=fltarr(n_elements(z))
lwc=fltarr(n_elements(z)) & iwc=fltarr(n_elements(z))
re_liq=fltarr(n_elements(z)) & re_ice=fltarr(n_elements(z))
cirrus_extinction=fltarr(n_elements(z))

common seedit, seed

;print, ze

if n_elements(seed) eq 0 then seed=1.e8*systime(/julian)

rms_flag=0

; the following moves us up to the laser cloud base to begin avoiding drizzle and snow and bugs

;print, 'cld base is ' ,cld_base, z[2],z[0]

if cld_base gt 0. then begin
k=0
while z[k] lt cld_base and k lt n_elements(z)-2 do begin

    ze_input[k]=-9999.
    ;print, 'in loop at ht ',z[k], ze_input[k]
    k=k+1
endwhile

;print, 'out of loop at ht ',z[k], ze_input[k], ze_input[k+1]

endif

ze=ze_input




for j=0,n_elements(z)-1 do begin
    if ze[j] gt -9999. then begin
ccm3_ice_frac=ccm3_cloud_iceliquid_fraction_param_function(temp[j])
;print, z[j], ' is height going into liquid param'
ccm3_liquid_profile[j]=(1.-ccm3_ice_frac)*ccm3_liquid_water_param_function(z[j], vap)
    endif
endfor

; What fraction of the lwp occurs at temps colder than freezing based on the CCM3 profile?

warm_water=0. & j=0. & total_sqrtZ_warm=0. & ccm3_cold_water=0.
for j=0,n_elements(z)-1 do begin
    if temp[j] gt 273. and ze[j] gt -9999. then begin
       warm_water=warm_water+ccm3_liquid_profile[j]
       total_sqrtZ_warm=total_sqrtZ_warm+(sqrt(ze[j]))
    endif

    if temp[j] le 273. and ze[j] gt -9999. then begin
       ccm3_cold_water=ccm3_cold_water+ccm3_liquid_profile[j]
    endif

endfor

    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,n_elements(z)-1 do begin
    if temp[j] gt 273. and ze[j] gt -9999. then begin
       lwc[j]=((sqrt(ze[j]))/total_sqrtZ_warm)*(warm_water_path/90.)
       ; if rms_flag eq 1 then begin
    ;  lwc[j]=lwc[j]+((randomn(seed,/normal))*lwc[j])
    ;  if lwc[j] le 0.0 then lwc[j]=((sqrt(ze[j]))/total_sqrtZ_warm)*(warm_water_path/90.)

       ;print,  ((randomn(seed,/normal))*lwc[j]) ,lwc[j], format='(e17.9,1x,e17.9)'
    ;   endif
       if lwc[j] lt 1.e-4 then lwc[j]=0.0

    endif else begin
         if ze[j] gt -9999. then lwc[j]=(cold_water_path*(ccm3_liquid_profile[j]/ccm3_cold_water))/90.
         if rms_flag eq 1 then begin
       lwc[j]=lwc[j]+((2.*randomn(seed,/normal))*lwc[j])
       if lwc[j] le 0.0 then lwc[j]=(cold_water_path*(ccm3_liquid_profile[j]/ccm3_cold_water))/90.
         endif

         if lwc[j] lt 1.e-4 then lwc[j]=0.0
    endelse

if 10.*alog10(ze[j]) lt 0.0 and ze[j] gt -9999. 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 -9999. and lwc[j] ge 1.e-4 then re_liq[j]=15.   ; was 8.4 6/29/05
endelse

rold=re_liq[j]

;   if rms_flag eq 1 then begin
;     re_liq[j]=re_liq[j]+((randomn(seed,/normal))*re_liq[j])
;     if re_liq[j] le 0.0 then re_liq[j]=rold
;   endif


endfor




; now convert the liquid profile to a weighting function
;    if lwp gt 0.0 and lwp lt 1.e4 then begin
;for j=0,n_elements(z)-1 do begin
;    if ze[j] gt -9999. then begin
;if total(ccm3_liquid_profile) gt 0.0 then begin
;     liquid_wt[j]=ccm3_liquid_profile[j]/total(ccm3_liquid_profile)
;endif else begin
;    liquid_wt[j]=0.0
;endelse
;lwc[j]=(liquid_wt[j]*(lwp))/90. ; assume 90m range gates;
;
;if 10.*alog10(ze[j]) lt 0.0 then begin
;    re_liq[j]=19.5*exp(0.0384*(10.*alog10(ze[j])))
;endif else begin
;    re_liq[j]=8.4
;endelse

;re_liq[j]=5.0

;print, liquid_wt[j],lwc[j],re_liq[j],((3./2.)*lwc[j]/re_liq[j])*90.,' is the liquid wt, lwc, and re'


;    endif
;endfor
;    endif

; 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

bz=-2.0 & ba=1.7
;IDL> print, r1   this is based on most data from 2000 iop with temp not in cc
;     0.509637     0.640465    0.0502367
;IDL> print, rmul1
 ;    0.660524
;IDL> print, mlr_coeff1
; 1.95209e-008
; 1.64473e-007
;-4.37153e-010
;IDL> print, mlr_const1
; 1.07064e-007

; 9/17/2004
;IDL> print, mlr_coeff1
; 8.76098e-008
; 2.67518e-010
;IDL> print, mlr_const1
;-1.52032e-008

;coef=-1.52032e-008 & coef=[coef,8.76098e-008] & coef=[coef,2.67518e-010]; & coef=[coef,-4.37153e-010]
coef=-8.434e-009 & coef=[coef,6.721e-008] & coef=[coef,2.409e-010]; & coef=[coef,-4.37153e-010] ; changed using King correction 6/10/2005
bz=-2.4 & ba=1.0 & bm=1.8 & bi=0.6

;print, ze

for j=0,n_elements(z)-1 do begin
    if ze[j] gt -9999. and temp[j] lt 273. and z[j] gt 2000. then begin


    if temp[j] gt 268. and 10.*alog10(ze[j]) gt 10. then begin
       zzz=10.^(10./10.)
    endif else begin
       zzz=ze[j]     ;0.631 ; remove 2 db for removal of liquid - this is a test 6/24/2005
       ;print, 'adjusting z in cond2'
    endelse


    if vdop[j] lt 0. then vdop[j]=0.5
    ;z_temp=((ze[j]))^0.6
    ;v_temp=vdop[j]*((zzz))^(ba/(6.+bz))
    v_temp=vdop[j]*(((zzz)))^((ba/((6.+bz))))
    ;v_temp2=vdop[j]*(((zzz/0.25)))^((ba/((6.+bz))))
    ;t_temp=611.*exp(((2.5e6)/461.5)*((1./273.)-(1./temp[j])))
    ;t_temp=temp[j]
    t_temp=611.*exp(((2.5e6)/461.5)*((1./273.)-(1./temp[j])))


    ;iwc[j]=coef[0]+(coef[1]*v_temp)+(coef[2]*t_temp)  ; commented out for new physical param 8/11/05
    ;iwc_test=coef[0]+(coef[1]*v_temp2)+(coef[2]*t_temp)


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)))

;print, zzz, vdop[j],iwc[j],iwc2, ' is old iwc and physical iwc'

    ;stop
    iwc[j]=iwc[j]*1.e6 ; converts to g/m3



    if iwc[j] lt 0. or lwp le 10. then begin
iwc[j]=liu_illingworth_cirrus_regression(10.*alog10(ze[j]), temp[j], 'mid')
;print, 'iwc lt 0 or lwp le 0 - using L&I in cond2'
;ice_fraction[l,j]=1.
;iwc_source[l,j]=landi_index
    endif







         if rms_flag eq 1 then begin
       iwc[j]=iwc[j]+((randomn(seed,/normal))*iwc[j])
       if iwc[j] le 0.0 then iwc[j]=coef[0]+(coef[1]*v_temp)+(coef[2]*t_temp)
         endif


    ;print, j,iwc[j],10.*alog10(zzz),10.*alog10(ze[j]),temp[j],vdop[j]

    ;re_ice[j]=15.
    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
    ;print, 're ice is ',re_ice[j],exp(-z[j]/5500.)

;   if rms_flag eq 1 then begin
;     re_ice[j]=re_ice[j]+((randomn(seed,/normal))*re_ice[j])
;     if re_ice le 0.0 then re_ice[j]=ccm3_cloud_ice_re_param_function(exp(-z[j]/5500.))
;   endif


    endif else begin
        iwc[j]=0.0
    re_ice[j]=0.

    endelse
endfor


; implement ciret2.  First see if profiles has cirrus.
k=0
while k lt n_elements(temp)-5 do begin


cirrus_flag=0 & max_dbz_temp=9999. & layer_top_temp=9999. & base_index=0 & top_index=0 & max_ze=-9999.
layer_base_index=-1
              while temp[k] gt 273. and k lt n_elements(temp)-2 do k=k+1
if k lt n_elements(temp)-2 and ze[k] gt -9999. then begin
;;;;;;;;;;;;;;;;;
                     while ze[k] gt -9999 and k lt n_elements(ze)-2 do k=k+1

if k lt n_elements(ze)-2 then begin
              while ze[k] eq -9999 and k lt n_elements(z)-2 do k=k+1    ; this gets us to cloud
while k lt n_elements(z)-2 and ze[k] gt -9999. do begin
kk=k
if ze[kk] gt -9999. then begin
    layer_top_temp=temp[kk] & layer_top_index=kk
    if layer_base_index eq -1 then layer_base_index=kk
       if ze[kk] gt max_ze then begin
    max_dbz_temp=temp[kk] & max_ze=ze[kk]
       endif
endif
k=k+1
endwhile

endif

;;;;;;;;;;;;;;;;;;;;
endif else begin    ; we hit the freezing level out of cloud

if k lt n_elements(ze)-2 then begin
              while ze[k] eq -9999 and k lt n_elements(z)-2 do k=k+1    ; this gets us to cloud
while k lt n_elements(z)-2 and ze[k] gt -9999. do begin
kk=k
if ze[kk] gt -9999. then begin
    layer_top_temp=temp[kk] & layer_top_index=kk
    if layer_base_index eq -1 then layer_base_index=kk
       if ze[kk] gt max_ze then begin
    max_dbz_temp=temp[kk] & max_ze=ze[kk]
       endif
endif
k=k+1
endwhile

endif

endelse

if max_dbz_temp le 273. and layer_top_temp le 240. and layer_base_index 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, layer_top_index do begin

if ze[kk] ne -9999. then begin

    if layer_top_index gt layer_base_index then begin
top_dist=(z[layer_top_index]-z[kk])/(z[layer_top_index]-z[layer_base_index])
    endif else begin
top_dist=90.
    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 3.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


rho_flag=0
if rho_air le 0. then begin

rho_flag=1

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])))

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=liu_illingworth_cirrus_regression(10.*alog10(ze[kk]), temp[kk], 'mid')

if abs(iwc_LI-(iwc_ret))/iwc_LI gt 5. then begin
    N0_ret_err=-9999.
    dge_ret=-9999.
    iwc_ret=-9999
    iwc[kk]=iwc_LI
    re_ice[kk]=40.
endif

; get extinction and write to output arrays

;Z=(10.^((double(dbz_in))/10.d))*1.d-12
;Vdq=(double(vdq__in))

       ;if N0_ret_err gt -9999. 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))
	
	if N0_ret_err eq -7777. then N0_ret_err=-9999.

if dge_ret gt 0.0 then begin
    iwc[kk]=iwc_ret
    re_ice[kk]=effective_radius_ret
endif
 ;      endif

endif   ; ze -9999
endfor   ; layer base index to layer top index

endif   ;if layer_base_index gt 0 then begin


k=k+1
endwhile    ;while k lt n_elements(temp)-5 do begin


;for j=0,n_elements(re_liq)-1 do begin
;if re_liq[j] gt 0.0 then re_liq[j]=re_liq[j]+(re_liq[j]*0.20)
;endfor

return
end