123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165 |
- #!/bin/bash
- #SBATCH -n 1
- #SBATCH -t 12:00:00
- #SBATCH -J post_perturbation
- #SBATCH -o slurm_post-%j.out
- #SBATCH -e slurm_post-%j.err
- set -o nounset
- set -o errexit
- set -x
- # -----------------------------------------------------------------------------
- # Author - F. Massonnet
- # Purpose- Process perturbations of DFS5.2 forcing and adds them
- # to the actual forcing.
- # What the script does:
- # 1) Interpolate from daily to 3hourly
- # 2) Add 29 February if necessary
- # 3) Add to the true forcing
- # -----------------------------------------------------------------------------
- # Which variable of the DFS5.2 forcing set has to be perturbed?
- # Possibilities: t2, q2, u10, v10, qlw, qsw, snow or precip
- var=qsw
- alpha=1.0 # VERY IMPORTANT: strength of the perturbation. 1 = same as interannual variability. 0.5 = half of it.
- # First and end years for which a perturbation has to be created.
- yearb=1973
- yeare=1973
- yearbp=1979 # First and end years defining the reference period on which
- yearep=2015 # the perturbations were created (must match those in preprocess.bash and in mkpert.R)
- nmb=1 # members to loop over. nmb = first member; nme = last one.
- # Usually the first member ("fc0") is the true forcing so it should *NOT* be perturbed.
- # Hence put nmb=1 and nme=25 if you want 25 perturbed forcings
- nme=25
- workdir=/esnas/scratch/$USER/TMP/TMP_${var}_26904/ # Where all perturbations are recorded
- # This folder was defined during the execution of preprocess.bash
- outtag=DFS5.2 # Name of the forcing after perturbations are created
- # Can be the same (DFS5.2) or modified (perturbed-DFS5.2 for instance)
- outdir=/esarchive/releases/fg/ocean/$outtag/
- cd $workdir
- echo "Workdir is $workdir"
- case ${var} in
- t2)
- min=100.0 # Min and max values allowed
- max=400.0
- freq=3hour # Frequency of availability
- ntim=2920 # Number of time steps in a year
- fvar=${var} # Name of the variable in the NetCDF
- ;;
- q2)
- min=0.0
- max=0.1
- freq=3hour
- ntim=2920
- fvar=${var}
- ;;
- u10)
- min=-100.0
- max=100.0
- freq=3hour
- ntim=2920
- fvar=${var}
- ;;
- v10)
- min=-100.0
- max=100.0
- freq=3hour
- ntim=2920
- fvar=${var}
- ;;
- qsw)
- min=0.0
- max=1000.0
- freq=1day
- ntim=365
- fvar=radsw
- ;;
- qlw)
- min=0.0
- max=1000.0
- freq=1day
- ntim=365
- fvar=radlw
- ;;
- snow)
- min=0.0
- max=0.001
- freq=1day
- ntim=365
- fvar=${var}
- ;;
- precip)
- min=0.0
- max=0.01
- freq=1day
- ntim=365
- fvar=${var}
- ;;
- *)
- echo "Variable $var unknown"
- exit
- esac
- for year in `seq ${yearb} ${yeare}`
- do
- for mm in `seq $nmb $nme`
- do
- if [ $mm == 0 ]
- then
- echo "WARNING!!!!"
- echo "Usually, at BSC, member 0 is booked to point towards the true forcing"
- echo "By creating a perturbed forcing named fc00, and then copying to the source"
- echo "directory, you are going to erase a symbolic link that points to the true forcing"
- echo "and thereby erase the true forcing!!"
- echo "Since this is highly dangerous, this script is aborting."
- echo "Contact francois.massonnet@bsc.es for further questions"
-
- exit
- fi
- m=$(printf "%02d" $mm)
- # There is a small trick here. If we have 3 days of data and ask for 3-hourly interpolation, we will have only 17 points and not 24.
- # We have in fact (ndays - 1) * 8 + 1. So the trick is to append the last day to the data twice and remove the last time frame.
- #
- # Extract the last time frame
- ncks -F -O -d time,365,365 pert_${var}_DFS5.2_${year}_fc${m}_ref${yearbp}-${yearep}.nc tmp.${year}.${m}.nc
- # Append it
- ncrcat -F -O pert_${var}_DFS5.2_${year}_fc${m}_ref${yearbp}-${yearep}.nc tmp.${year}.${m}.nc pert_${var}_DFS5.2_${year}_fc${m}.nc.1
- # Set time axis
- cdo settaxis,${year}-01-01,00:00,1day pert_${var}_DFS5.2_${year}_fc${m}.nc.1 pert_${var}_DFS5.2_${year}_fc${m}.nc.2
- # Interpolate in time
- cdo inttime,${year}-01-01,00:00,${freq} pert_${var}_DFS5.2_${year}_fc${m}.nc.2 pert_${var}_DFS5.2_${year}_fc${m}.nc.3
- # Remove the last time frame
- ncks -F -O -d time,1,${ntim} pert_${var}_DFS5.2_${year}_fc${m}.nc.3 pert_${var}_DFS5.2_${year}_fc${m}.nc.4
- # Add the desired fraction "alpha" of the perturbation to the true forcing
- cdo add -mulc,${alpha} pert_${var}_DFS5.2_${year}_fc${m}.nc.4 ${var}_DFS5.2_${year}.nc ${var}_fc${m}_DFS5.2_${year}.nc.0
-
- # Physical bounds
- cdo setrtoc,-10000000000,${min},${min} ${var}_fc${m}_DFS5.2_${year}.nc.0 ${var}_fc${m}_DFS5.2_${year}.nc.1
- cdo setrtoc,${max},10000000000,${max} ${var}_fc${m}_DFS5.2_${year}.nc.1 ${var}_fc${m}_DFS5.2_${year}.nc.2
- # Set time units to allow nice reading
- cdo settunits,years ${var}_fc${m}_DFS5.2_${year}.nc.2 ${outdir}/${var}_fc${m}_${outtag}_${year}.nc
- # Add description in the header
- ncatted -O -a description,${fvar},a,c,"Perturbed version of DFS5.2 variable ${fvar} for year ${year} (member fc${m}). Strength of perturbation is ${alpha} times the year-to-year differences estimated over the ${yearbp}-${yearep} reference period. For more details: francois.massonnet@bsc.es" ${outdir}/${var}_fc${m}_${outtag}_${year}.nc
- chmod 777 ${outdir}/${var}_fc${m}_${outtag}_${year}.nc
- rm -f tmp.${year}.${m}.nc pert_${var}_DFS5.2_${year}_fc${m}.nc.? ${var}_fc${m}_DFS5.2_${year}.nc.?
- done
- done
- echo "SCRIPT POSTPROCESS.BASH FINISHED"
|