#kinewdc-d.ss -- kinematic for  Divergence Correction
#                differential part.
#
# this is to set the horizontal divergence and vorticity. Even though 
# divergence is corrected in radfills, it has to be reset 
# after using cdffill. 
# In the next step (kinewcdc-i.ss or kinewsdc-i.ss) the column 
# divergence correction will be perform.
#
#-------------------------------------------------------------

#set -x

case $# in
  0) ;;
  *) echo "usage kine-d.ss < stdin > stdout"; exit 1; ;;
esac

#
#

#-------------------------------------------------------------------
# Divergence Correction and Vertical Mass Current:
#
# from smothed u & v creates(see bellow):
# unc_div_hvel -- uncorrected divergence of horizontal velocity 
# vor          -- total vorticity
# unc_div_hmf -- uncorrected divergence of the horizontal mass flux
# div_hmf     -- divergence of the horizontal mass flux
# hmc         -- horizontal mass current = area integral of div_hmf
# vmf         -- vertical integration of hmc / harea_max
# dmf         -- detrained mass flux = horizontal area average of div_hmf
#                                    = area integral of div_hmf / harea_max
#----------------------------------------------------------------------
# A) Colum-BY-Column Divergence Correction
#   1) From density and the uncorrected divergence of the horizontal 
#      velocity(unc_div_hu), calculate the uncorrected divergence
#      of the horizontal mass flux (unc_div_hmf);  hmf = rho*hu
#   2) perform a column-by-column divergence correction on (unc_div_hmf)
#      to obtain the colum-corrected divergence of the horizontal
#      mass flux (div_hmf). 
#   3) The area integral of div_hmf can be though as a line integral
#      of the horizontal mass flux. Hence, we can think of such integral
#      as being proportional to the net horizontal mass flux at height Z.
#      In fact, multiplying the value of this integral by the area  of
#      integration*, results in the net mass per unit time exported from
#      the system at height Z. Therefore, we call this integral the 
#      horizontal mass current(hmc).
#
# B) System Divergence Correction of the hmc
#    This is performed to compensate for numerical errors during the area 
#    integration. 
#   
#    Since hmc is the area integral of div_hmf, dividing hmc by the horizontal
#    area of the system results in the average of div_hmf. We call this
#    average of the divergence of the horizontal mass flux, the 
#    detrained mass flux(dmf)
#
# C) The twice corrected dmf is finally used to compute the 
#    vertical mass flux(vmf) via vertical integration from bottom. 
#
#------------------------------------------------------------
#
#Require:cdffmean, cdfinteg_exp, cdfdefint_exp

SCALE=1e-3;

#
# Calculate the uncurrected divergence of horizontal velocity
# and vorticity fields.
#
cdfderiv dudx u x |\
cdfderiv dvdy v y |\
cdfderiv dudy u y |\
cdfderiv dvdx v x |\
cdfmath "dudx dvdy + unc_div_hvel = " |\
cdfmath "dvdx dudy - vor = "  |\

# Start Colum divergence correction ..
#
# Calculate the uncurrected horizontal mass current density
#
cdfmath "unc_div_hvel rho * unc_div_hmf =" |\
#
#calculate constant for the column-by-column correction.
#
cdffmean z .m unc_div_hmf |\
#
#make the correction
#
cdfmath "unc_div_hmf unc_div_hmf.m - div_hmf =" |\
cdfmath "div_hmf rho / div_hvel ="  |\
#
#Colum divergence correction .. done!
#

#
#extract fields of interest
#
cdfextr u v vor div_hvel div_hmf dbz |\

#
#create name for integrated fields by scaling down the extracted fields
#
cdfmath " u ${SCALE} * u_ave = " |\
cdfmath " v ${SCALE} * v_ave = " |\
cdfmath " div_hmf ${SCALE} * hmc = " |\
cdfmath " vor ${SCALE} * cir = " |\

#
#creating field area of reflectivity 
#
cdfmath " dbz ${SCALE} 0 ? refa =" |\

#
#extract fields to integrete 
#
cdfextr u_ave v_ave hmc cir refa |\

#
#creating scaled auxiliar area fields 
#
cdfmath " u_ave ${SCALE} 0 ? a_u = v_ave ${SCALE} 0 ? a_v = " |\
cdfmath " hmc ${SCALE} 0 ? a_hmc = " |\

#
# Perform areal integration --
#
cdfdefint_exp x |\
cdfdefint_exp y |\

#
# Get average values for the wind velocity components
#
cdfmath "u_ave a_u / u_ave = v_ave a_v / v_ave = " |\

#
# get maximun area of the system
#
cdfmax z .max a_hmc |\
cdfmath "a_hmc.max hsysarea.max =" |\
cdfextr -r a_hmc.max |\

#Make System divergence correction to compensate for numerical errors
#during areal integration.
#
cdffmean z .m hmc |\
cdfmath "hmc hmc.m - hmc =" |\
cdfextr -r hmc.m |\
#
#System Divergence Correction .. done!
#

#
#get detrained mass flux in (kg/m^3)/ks
#
cdfmath "hmc hsysarea.max / dmf = " |\

#Get vertical mass flux in (kg/m^3)(m/ks)
# the minus comes from the mass continuity eq.
cdfinteg_exp -b vmf dmf z  |\
cdfmath "vmf -1000. * hsysarea.max / vmf = " |\

#Scale dmf and vmf
cdfmath "dmf 1e2 * dmf = vmf 1e2 * vmf = " 



#end