Blame | Last modification | View Log | Download | RSS feed
#!/bin/csh
# -------------------------------------------------------------------
# Set some variables and paths
# -------------------------------------------------------------------
# Jump to specified step
set step="help"
set param1=
set param2=
set param3=
if ( $#argv == 1 ) then
set step=$1
endif
if ( $#argv == 2 ) then
set step=$1
set param1=$2
endif
if ( $#argv == 3 ) then
set step=$1
set param1=$2
set param2=$3
endif
if ( $#argv == 4 ) then
set step=$1
set param1=$2
set param2=$3
set param3=$4
endif
set file1=${param1}
set file2=${param2}
# Set base directory for programmes
set bdir=${DYN_TOOLS}/inversion/
# Set name of the parameter file, coastline file and sample directory
set parafile="${PWD}/inversion.param"
set coastfile="${bdir}/prep/coastline.dat"
set sampledir="/net/bio/atmosdyn/erainterim/cdf/2006/01/"
# ---------------------------------------------------------------------------
# Installation, help and sample
# ---------------------------------------------------------------------------
# Get parameter file, if not yet available
if ( ! -f ${parafile} ) then
cp ${bdir}/inversion.param .
endif
# Extract parameters from parameter file
set progex=${bdir}/inversion.perl
set date=`${progex} ${parafile} date | awk '{ print $2}'`
set nofiter=`${progex} ${parafile} n_of_iteration | awk '{ print $2}'`
set save=`${progex} ${parafile} save_iteration | awk '{ print $2}'`
set idir=`${progex} ${parafile} inp_dir | awk '{ print $2}'`
set rdir=`${progex} ${parafile} run_dir | awk '{ print $2}'`
set odir=`${progex} ${parafile} out_dir | awk '{ print $2}'`
# Create the needed directories
if ( ${step} == "inst" ) then
if ( ! -d ${idir} ) mkdir ${idir}
if ( ! -d ${rdir} ) mkdir ${rdir}
if ( ! -d ${odir} ) mkdir ${odir}
endif
# Create the needed directories
if ( ${step} == "help" ) then
echo
echo "Installation"
echo " inst: Creates the input, run and output directory; copy template parameter file"
echo
echo "Sample case study"
echo " sample: Copy all files for a sample case study"
echo
echo "Preparing input files [prep]"
echo " prep0: Calculate S file with PV and TH [ P -> S ]"
echo " prep1: Interpolate onto height levels [ P,Z -> H ]"
echo " prep2: Rotate into local cartesian coordinate system [ H -> R ]"
echo " prep3: Add TH,PV,NSQ and RHO to the data file [ -> R ]"
echo " prep4: Define modified and anomaly PV field and boundary values [ -> R ]"
echo " prep5: Reduce the domain size and split the input files [ R -> ORG,MOD,ANO,REF]"
echo " prep6: Add the reference profile [MOD -> REF]"
echo " prep7: Add coastlines to REF file [-> REF]"
echo " prep8: Move the files to the run directory [ORG,MOD,ANO,REF -> ]"
echo
echo "Perform the PV inversion [pvin]"
echo " pvin1: Add NSQ, TH, RHO, and PV [ -> MOD]"
echo " pvin2: Change Ertel's PV anomaly into one of quasi-geostrophic PV [MOD -> ANO]"
echo " pvin3: Inversion of quasi-geostrophic PV anomaly [ANO -> ANO]"
echo " pvin4: Subtract anomaly from MOD file [MOD-ANO -> MOD]"
echo " pvin5: Keep iterative steps if save flag is set"
echo
echo "Postprocessing [post]"
echo " post1: Copy needed files from input and run directory [P,MOD -> ]"
echo " post2: Rotate from quasi-cartesian coordinate frame to lat/lon system [ MOD -> GEO ]"
echo " post3: Bring modified fields back to P file [ P+GEO -> P ]"
echo " post4: Calculate S file with PV and TH [ P -> S ]"
echo " post5: Make clean"
echo " post6: Calculate difference to original P and S file"
echo
echo "Diagnostic Tools"
echo " diag1: Check the consitency of the boundary conditions"
echo " diag2: Calculate the quasi-geostrophic PV (diag2 [ORG|MOD|ANO])"
echo " diag3: Get the difference between two files (diag3 [ORG|MOD|ANO] [ORG|MOD|ANO])"
echo " diag4: Calculate geopotential with hydrostatic equation (diag4 [P|ORG|MOD|ANO])"
echo " diag5: Q vector analysis; geostrophic wind balance [ORG|MOD|ANO] [GEO|DIV_UV|QVEC|DIV_Q]"
echo
echo "Modifying Pv anomalies"
echo " mod1: Stretching and amplitude changes for anomalies [MOD]"
endif
# Copy sample files (if specified)
if ( ${step} == "sample" ) then
\cp ${sampledir}/P20060116_18 ${idir}
\cp ${sampledir}/eraint_cst ${idir}
\cp ${sampledir}/Z20060116_18 ${idir}
\cp ${sampledir}/pressure_cst ${idir}
endif
# ---------------------------------------------------------------------------
# Preparatory steps
# ---------------------------------------------------------------------------
# Change to data directory
cd ${idir}
# Step 0: Calculate S file with PV and TH (not in prep mode)
if ( ${step} == "prep0" ) then
\rm -f S${date}
p2s P${date} TH PV
endif
# Step 1: Interpolate onto height levels
if ( ${step} == "prep1" | ${step} == "prep" ) then
echo "P${date}" >! fort.10
echo "Z${date}" >> fort.10
echo "H${date}" >> fort.10
${bdir}/inversion.perl ${parafile} p2z >> fort.10
echo "U U P${date}" >> fort.10
echo "V V P${date}" >> fort.10
echo "T T P${date}" >> fort.10
echo "Q Q P${date}" >> fort.10
${bdir}/prep/p2z
\rm -f fort.10
endif
# Step 2: Rotate into local cartesian coordinate system
if ( ${step} == "prep2" | ${step} == "prep" ) then
echo "H${date}" >! fort.10
echo "R${date}" >> fort.10
${bdir}/inversion.perl ${parafile} rotate_grid >> fort.10
echo "5" >> fort.10
echo "ORO" >> fort.10
echo "U.V" >> fort.10
echo "P" >> fort.10
echo "T" >> fort.10
echo "Q" >> fort.10
${bdir}/prep/rotate_grid
\rm -f fort.10
endif
# Step 3: Add TH,PV,NSQ and RHO to the data file
if ( ${step} == "prep3" | ${step} == "prep" ) then
echo "TH" >! fort.10
echo "R${date}" >> fort.10
echo "R${date}" >> fort.10
echo "5 " >> fort.10
${bdir}/prep/z2s
echo "PV" >! fort.10
echo "R${date}" >> fort.10
echo "R${date}" >> fort.10
echo "5 " >> fort.10
${bdir}/prep/z2s
echo "NSQ" >! fort.10
echo "R${date}" >> fort.10
echo "R${date}" >> fort.10
echo "5 " >> fort.10
${bdir}/prep/z2s
echo "RHO" >! fort.10
echo "R${date}" >> fort.10
echo "R${date}" >> fort.10
echo "5 " >> fort.10
${bdir}/prep/z2s
\rm -f fort.10
endif
# Step 4: Set the modified PV field and boundary values
if ( ${step} == "prep4" | ${step} == "prep" ) then
echo "R${date}" >! fort.10
${bdir}/inversion.perl ${parafile} def_anomaly >> fort.10
${bdir}/prep/def_anomaly
\rm -f fort.10
endif
# Step 5: Reduce the domain size and split the input files
if ( ${step} == "prep5" | ${step} == "prep" ) then
echo "PV PV R${date} ORG_${date}" >! fort.10
echo "U U R${date} ORG_${date}" >> fort.10
echo "V V R${date} ORG_${date}" >> fort.10
echo "TH TH R${date} ORG_${date}" >> fort.10
echo "Q Q R${date} ORG_${date}" >> fort.10
echo "P P R${date} ORG_${date}" >> fort.10
echo "T T R${date} ORG_${date}" >> fort.10
echo "PV_FILT PV_AIM R${date} MOD_${date}" >> fort.10
echo "U U R${date} MOD_${date}" >> fort.10
echo "V V R${date} MOD_${date}" >> fort.10
echo "Q Q R${date} MOD_${date}" >> fort.10
echo "P P R${date} MOD_${date}" >> fort.10
echo "T T R${date} MOD_${date}" >> fort.10
echo "NSQ NSQ R${date} MOD_${date}" >> fort.10
echo "RHO RHO R${date} MOD_${date}" >> fort.10
echo "TH TH R${date} MOD_${date}" >> fort.10
echo "PV_ANOM PV R${date} ANO_${date}" >> fort.10
echo "TH_ANOM TH R${date} ANO_${date}" >> fort.10
echo "UU_ANOM U R${date} ANO_${date}" >> fort.10
echo "VV_ANOM V R${date} ANO_${date}" >> fort.10
echo "ORO ORO R${date} REF_${date}" >> fort.10
echo "X X R${date} REF_${date}" >> fort.10
echo "Y Y R${date} REF_${date}" >> fort.10
echo "LAT LAT R${date} REF_${date}" >> fort.10
echo "LON LON R${date} REF_${date}" >> fort.10
echo "CORIOL CORIOL R${date} REF_${date}" >> fort.10
${bdir}/prep/cutnetcdf
\rm -f fort.10
endif
# Step 6: Add the reference profile
if ( ${step} == "prep6" | ${step} == "prep" ) then
\rm -f fort.10
echo "MOD_${date}" >! fort.10
echo "REF_${date}" >> fort.10
${bdir}/inversion.perl ${parafile} ref >> fort.10
${bdir}/prep/ref_profile
\rm -f fort.10
endif
# Step 7: Add coastlines to REF file
if ( ${step} == "prep7" | ${step} == "prep" ) then
\rm -f fort.10
echo \"REF_${date}\" >! fort.10
echo \"${coastfile}\" >> fort.10
${bdir}/inversion.perl ${parafile} coastline >> fort.10
${bdir}/prep/coastline
\rm -f fort.10
endif
# Step 8: Move the files to the run directory
if ( ${step} == "prep8" | ${step} == "prep" ) then
\mv MOD_${date} MOD_${date}_cst ${rdir}
\mv ORG_${date} ORG_${date}_cst ${rdir}
\mv ANO_${date} ANO_${date}_cst ${rdir}
\mv REF_${date} REF_${date}_cst ${rdir}
\rm -f fort.10
endif
# ---------------------------------------------------------------------------
# Inversion
# ---------------------------------------------------------------------------
# Change to data directory
cd ${rdir}
# Start loop
set count=0
loop:
# Step 1: Add NSQ, TH, RHO, and PV to MOD file, take grid from REF file
if ( ${step} == "pvin1" | ${step} == "pvin" ) then
echo "NSQ" >! fort.10
echo "MOD_${date}" >> fort.10
echo "REF_${date}" >> fort.10
echo "5 " >> fort.10
${bdir}/pvin/z2s
echo "RHO" >! fort.10
echo "MOD_${date}" >> fort.10
echo "REF_${date}" >> fort.10
echo "5 " >> fort.10
${bdir}/pvin/z2s
echo "TH" >! fort.10
echo "MOD_${date}" >> fort.10
echo "REF_${date}" >> fort.10
echo "5 " >> fort.10
${bdir}/pvin/z2s
echo "PV" >! fort.10
echo "MOD_${date}" >> fort.10
echo "REF_${date}" >> fort.10
echo "5 " >> fort.10
${bdir}/pvin/z2s
\rm -f fort.10
endif
# Step 2: Change Ertel's PV anomaly into an anomaly of quasi-geostrophic PV
if ( ${step} == "pvin2" | ${step} == "pvin" ) then
echo "MOD_${date}" >! fort.10
echo "REF_${date}" >> fort.10
echo "ANO_${date}" >> fort.10
${bdir}/pvin/pv_to_qgpv
\rm -f fort.10
endif
# Step 3: Inversion of quasi-geostrophic PV anomaly with Neumann boundary
if ( ${step} == "pvin3" | ${step} == "pvin" ) then
echo "ANO_${date}" >! fort.10
echo "REF_${date}" >> fort.10
${bdir}/pvin/inv_cart
\rm -f fort.10
endif
# Step 4: Prepare the output of the inversion for next iteration step
if ( ${step} == "pvin4" | ${step} == "pvin" ) then
\rm -f fort.10
echo "MOD_${date}" >! fort.10
echo "ANO_${date}" >> fort.10
${bdir}/inversion.perl ${parafile} prep_iteration >> fort.10
${bdir}/pvin/prep_iteration
\rm -f fort.10
endif
# Step 5: Keep iterative steps if save flag is set
if ( ${step} == "pvin5" | ${step} == "pvin" ) then
if ( "${save}" == "yes" ) then
set pre=''
if ( ${count} < 10 ) then
set pre='0'
endif
\cp MOD_${date} MOD_${date}_${pre}${count}
\cp ANO_${date} ANO_${date}_${pre}${count}
endif
endif
# End loop for iterations
if ( ${step} == "pvin" ) then
@ count = ${count} + 1
if ( ${count} < ${nofiter} ) goto loop
endif
# ---------------------------------------------------------------------------
# Postprocessing
# ---------------------------------------------------------------------------
# Change to output directory
cd ${odir}
# Step 1: Copy needed files from input and run directory
if ( ${step} == "post1" | ${step} == "post" ) then
ln -sf ${rdir}/MOD_${date} .
ln -sf ${rdir}/MOD_${date}_cst .
ln -sf ${rdir}/ORG_${date} .
ln -sf ${rdir}/ORG_${date}_cst .
\cp ${idir}/P${date} .
set cfn = `getcfn P${date}`
\cp ${idir}/${cfn} .
endif
# Step 2: Rotate from quasi-cartesian coordinate frame to lat/lon system
if ( ${step} == "post2" | ${step} == "post" ) then
echo "MOD_${date}" >! fort.10
echo "GEO.MOD_${date}" >> fort.10
${bdir}/inversion.perl ${parafile} rotate_lalo >> fort.10
echo "3" >> fort.10
echo "T" >> fort.10
echo "U.V" >> fort.10
echo "P" >> fort.10
${bdir}/post/rotate_lalo
\rm -f fort.10
echo "ORG_${date}" >! fort.10
echo "GEO.ORG_${date}" >> fort.10
${bdir}/inversion.perl ${parafile} rotate_lalo >> fort.10
echo "3" >> fort.10
echo "T" >> fort.10
echo "U.V" >> fort.10
echo "P" >> fort.10
${bdir}/post/rotate_lalo
\rm -f fort.10
endif
# Step 3: Bring modified fields back to P file
if ( ${step} == "post3" | ${step} == "post" ) then
\rm -f fort.10
echo "P${date}" >! fort.10
echo "GEO.MOD_${date}" >> fort.10
echo "GEO.ORG_${date}" >> fort.10
${bdir}/inversion.perl ${parafile} add2p >> fort.10
${bdir}/post/add2p
\rm -f fort.10
endif
# Step 4: Calculate S file with PV and TH
if ( ${step} == "post4" | ${step} == "post" ) then
\rm -f S${date}
p2s P${date} TH PV
endif
# Step 5: Make clean
if ( ${step} == "post5" | ${step} == "post" ) then
\rm -f MOD_${date} MOD_${date}_cst
\rm -f GEO_${date} GEO_${date}_cst
endif
# Step 6: Get difference of original and modified P,S files
if ( ${step} == "post6" | ${step} == "post" ) then
\rm -f levs
echo "950." >! levs
echo "925." >> levs
echo "900." >> levs
echo "875." >> levs
echo "850." >> levs
echo "825." >> levs
echo "800." >> levs
echo "775." >> levs
echo "750." >> levs
echo "725." >> levs
echo "700." >> levs
echo "675." >> levs
echo "650." >> levs
echo "625." >> levs
echo "600." >> levs
echo "575." >> levs
echo "550." >> levs
echo "525." >> levs
echo "500." >> levs
echo "475." >> levs
echo "450." >> levs
echo "425." >> levs
echo "400." >> levs
echo "375." >> levs
echo "350." >> levs
echo "325." >> levs
echo "300." >> levs
echo "275." >> levs
echo "250." >> levs
echo "225." >> levs
echo "200." >> levs
echo "175," >> levs
echo "150." >> levs
echo "125." >> levs
echo "100." >> levs
echo " 75." >> levs
echo " 50." >> levs
\rm -f vars
echo "U P" >! vars
echo "V P" >> vars
echo "T P" >> vars
echo "Z P" >> vars
echo "TH S" >> vars
echo "PV S" >> vars
nput2p ${date} pr_cst levs vars
ncks -A -v PS P${date} L${date}
\mv -f L${date} L${date}.1
\cp levs vars ${idir}/
cd ${idir}
nput2p ${date} pr_cst levs vars
ncks -A -v PS P${date} L${date}
\rm levs
\rm vars
\mv L${date} ${odir}/L${date}.0
cd ${odir}
ncdiff L${date}.1 L${date}.0 L${date}
\rm -f L${date}.1
\rm -f L${date}.0
endif
# ---------------------------------------------------------------------------
# Diagnotic Tools
# ---------------------------------------------------------------------------
# Step 1: Check the consistency of the boundary conditions (diag1)
if ( ${step} == "diag1" ) then
cd ${rdir}
\rm -f fort.10
echo "ANO_${date}" >! fort.10
echo "REF_${date}" >> fort.10
${bdir}/diag/check_boundcon
\rm -f fort.10
endif
# Step 2: Calculate the quasi-geostrophic PV (diag2 [ORG|MOD|ANO]
if ( ${step} == "diag2" ) then
cd ${rdir}
\rm -f fort.10
echo "${file1}_${date}" >! fort.10
echo "REF_${date}" >> fort.10
${bdir}/diag/calc_qgpv
\rm -f fort.10
endif
# Step 3: Get difference between two files (diag3 [ORG|MOD|ANO] - [ORG|MOD|ANO])
if ( ${step} == "diag3" ) then
cd ${rdir}
\rm -f fort.10
echo "${file1}_${date}" >! fort.10
echo "${file2}_${date}" >> fort.10
echo "DIA_${date}" >> fort.10
${bdir}/diag/difference
\rm -f fort.10
endif
# Step 4: Calculate geopotential with hydrostatic equation (diag4)
if ( ${step} == "diag4" ) then
cd ${odir}
\rm -f fort.10
echo "P${date}" >! fort.10
echo "P${date}" >> fort.10
${bdir}/diag/hydrostatic
\rm -f fort.10
endif
# Step 5: Q vector analysis / geostrophic balance check
if ( ${step} == "diag5" ) then
cd ${rdir}
\rm -f fort.10
echo "${param3}" >! fort.10
echo "${file1}" >> fort.10
echo "${file2}" >> fort.10
echo 1 >> fort.10
${bdir}/diag/qvec_analysis
\rm -f fort.10
endif
# ---------------------------------------------------------------------------
# Modifying tools
# ---------------------------------------------------------------------------
if ( ${step} == "mod1" ) then
cd ${idir}
\rm -f fort.10
set commandstr=$2
echo "R${date}" >! fort.10
echo \"${commandstr}\" >> fort.10
${bdir}/spec/modify_anomaly
\rm -f fort.10
endif