grapesinput2grads.sh

grapesinput2grads.sh

#!/bin/bash

#########################################################################
# Author: 凌科
# Create Time: 2025/10/28
# Description: 从grapesinput提取变量，按照grads的要求再保存到grapesinput.dat中，并生成grapesinput.ctl文件
# Version: 1.0
#########################################################################
set -e


s_we=1
e_we=1163
s_sn=1
e_sn=543
s_vert=1
e_vert=65
num_soil_layers=5


max_i=$(( e_we - s_we + 1 ))
max_j=$(( e_sn - s_sn + 1 ))
var1_max_k=$(( (e_vert+1)-(s_vert-1)+1 ))
var2_max_k=$num_soil_layers

# 平面字节数
ixjx4=$(( max_i * max_j * 4 ))
# 一层大小，加前4字节长度和后4字节长度
layer_size=$(( ixjx4 + 8 ))

output_dir="tmp_output"
layer_seq=0

# 检查参数个数
if [ $# -lt 1 ]; then
    echo "usage: $0 grapesinput"
    exit 1
fi

input_file=$1




##########################################################
# 读取时间信息，年、月、日、时、分，生成grapesinput.ctl文件是会用到
read_integer() {

        input_file=$1
        offset=$2

        # 读取4个字节的十六进制字符串
        hex_string=$(xxd -s "$offset" -l 4 -p "$input_file")
        integer_value=$(printf "%d" "0x$hex_string")

        echo $integer_value
}

year=$(read_integer $input_file 4)
month=$(read_integer $input_file 8)
day=$(read_integer $input_file 12)
hour=$(read_integer $input_file 16)
minute=$(read_integer $input_file 20)

echo "year: $year"
echo "month: $month"
echo "day: $day"
echo "hour: $hour"
echo "minute: $minute"

##########################################################
# 创建输出目录
mkdir -p $output_dir

##########################################################
# 生成长度信息文件

len_bin_file="$output_dir/len.bin"

# 将整数格式化为 8 位十六进制字符串（4 字节）
hex_string=$(printf "%08X" "$layer_size")

# 将十六进制字符串转换为二进制数据并写入文件
echo -n "$hex_string" | xxd -r -p > $len_bin_file


###########################################################
copy_layer() {

        start_pos=$1
        k=$2
        max_k=$3
        layer_seq=$4

        # 创建临时目录
        tmp_dir="$output_dir/$layer_seq"
        mkdir -p $tmp_dir


        layer_pids=()
        # 由于grapesinput文件中是i k j顺序，提取一层，即i j平面，先提取j个i，再合并
        for ((j = 0; j < max_j; j++)); do

                # 由于grapesinput里数据最小单元是4字节real类型，所以bs就为4
                # 拷贝X轴1维数据，所以叫bar
                bar_start_pos=$(( start_pos + k * max_i * 4 + j * max_k * max_i * 4 ))
                skip_size=$(( bar_start_pos / 4 ))
                count_size=$(( max_i * 4 /4 ))

                bar_filename=$(printf "%06d" "$j")
                output_file="$tmp_dir/$bar_filename"

                dd if="$input_file" of="$output_file" bs=4 skip=$skip_size count=$count_size &
                layer_pids+=($!)

        done

        # 等待执行完成
        for pid in "${layer_pids[@]}"; do
            wait $pid
        done

        # 合并
        sync
        cat $len_bin_file $tmp_dir/* $len_bin_file > $output_dir/LAYER_$(printf "%06d" "$layer_seq")
        sync

        # 删除
        rm -rf $tmp_dir
}



###########################################################
# 从grapesinput文件中提取变量信息
# 19 个三维变量

var1_size=$(( var1_max_k * ixjx4 + 8 ))

for ((n = 0; n < 19; n++)); do

        for ((k = 0; k < $var1_max_k; k++)); do
                # 三维变量起始位置
                start_pos=$(( 28 + n * var1_size + 4 ))

                # 并行写入
                copy_layer $start_pos $k $var1_max_k $layer_seq &

                # 自增1
                layer_seq=$(( layer_seq + 1 ))
        done

done



###########################################################
# 2 个三维变量

var2_size=$(( var2_max_k * ixjx4 + 8 ))

for ((n = 0; n < 2; n++)); do

        for ((k = 0; k < $var2_max_k; k++)); do
                # 三维变量起始位置
                start_pos=$(( 28 + 19 * var1_size + n * var2_size + 4 ))

                # 并行写入
                copy_layer $start_pos $k $var2_max_k $layer_seq &

                # 自增1
                layer_seq=$(( layer_seq + 1 ))
        done


done

###########################################################
# 40 个二维变量

output_file=$output_dir/LAYER_$(printf "%06d" "$layer_seq")
start_pos=$(( 28 + 19 * var1_size + 2 * var2_size ))
length=$(( 40 * layer_size ))

skip_size=$(( start_pos / 4 ))
count_size=$(( length / 4 ))
dd if="$input_file" of="$output_file" bs=4 skip=$skip_size count=$count_size


###########################################################
# 等待执行完成
wait


# 合并
cat $output_dir/LAYER_* > $output_dir/grapesinput.dat
sync

# 删除
rm $output_dir/LAYER_*
rm $len_bin_file

cat > ${output_dir}/${input_file}.ctl << EOF
dset ^grapesinput.dat
options  template  sequential big_endian
title grapesinput
undef -9.99E+33
xdef   $e_we  linear   110.33    0.005
ydef   $e_sn  linear    21.65    0.005
zdef   $var1_max_k linear 1 1
tdef   1 linear $(printf "%02d" $hour)Z${day}$(date -d "$year-$month-01" +"%b" | tr '[:upper:]' '[:lower:]')${year}  1hr
vars 61
dthdz $var1_max_k 0 dthdz
dpdz $var1_max_k 0 dpdz
dpwdz $var1_max_k 0 dpwdz
thref $var1_max_k 0 potential tempertature
piref $var1_max_k 0 piref
zsx $var1_max_k 0 zsx
zsy $var1_max_k 0 zsy
pzsx $var1_max_k 0 pzsx
pzsy $var1_max_k 0 pzsy
zz $var1_max_k 0 zz
pip $var1_max_k 0 perturbed PI
pi $var1_max_k 0 PI
u $var1_max_k 0 U wind component
v $var1_max_k 0 V wind component
wzet $var1_max_k 0 wzet
thp $var1_max_k 0 perturbed potential tempertature
th $var1_max_k 0 potential tempertature
moist_2 $var1_max_k 0 moist_2
dtdz $var1_max_k 0 dtdz
tslb 5 0 tslb
smois 5 0 smois
zs 0 0 zs
zst 0 0 zst
fisx 0 0 fisx
fisy 0 0 fisy
pfisx 0 0 pfisx
pfisy 0 0 pfisy
tsk 0 0 tsk
ps 0 0 ps
psea 0 0 psea
u10 0 0 u10
v10 0 0 v10
t2 0 0 t2
q2 0 0 q2
ht 0 0 ht
tmn 0 0 tmn
xland 0 0 xland
lu_index 0 0 lu_index
snowc 0 0 snowc
xice 0 0 xice
xlat 0 0 xlat
xlong 0 0 xlong
soil_type 0 0 soil_type
vegfra 0 0 vegfra
uslope 0 0 uslope
vslope 0 0 vslope
wslope 0 0 wslope
albedo 0 0 albedo
oc12d 0 0 oc12d
var2d 0 0 var2d
oa1 0 0 oa1
oa2 0 0 oa2
oa3 0 0 oa3
oa4 0 0 oa4
ol1 0 0 ol1
ol2 0 0 ol2
ol3 0 0 ol3
ol4 0 0 ol4
radd 0 0 radd
radl 0 0 radl
bkdd 0 0 bkdd
endvars
EOF