老夫聊发少年狂，试试Python气象绘图。
现有数据为①记录能谱密度的二进制spect.dat；②记录均一化波数、均一化角波数和波长的aol.txt（奇怪的是aol.txt自动换成5个数值一行）。

import numpy as np
import struct
import matplotlib.pyplot as plt

def getbin(filename):
    with open(filename,mode='rb') as f:
        tmp1=f.read()
        fsize=f.tell()
    tmp11=list(struct.unpack('f'*26*479,tmp1)) 
#python没有二进制格式，所以需要用struck.unpack来把字符串转译成实数，f指float/float32
#26层，479点
    tmp12=np.reshape(tmp11,(26,479))
#从一维列表变形成二维列表，注意python是行优先
    sp=np.mean(tmp12,axis=0)
#按层平均
    return sp
#传出动能谱密度

def getticks(filename):
    with open(filename,mode='r') as f:
        data=f.readlines()
    d1=[]
    for line in data:
        arr=list(map(float,line.split()))
        d1.extend(arr)
#不知道为什么，fortran存储的文本文件变成了5个数值一行，所以用f.readlines()先把文本中所
#有的值读到data中
#然后对data中每一行按空格分割(split),转化为实数(float),针对整个列表实施操作，转成新列表
    d2=np.reshape(d1,(3,479))
    return  d2[1][:],d2[2][:]
#传出标准化角波数和波长

def main():
    sp_noskeb=getbin("spectra_noskeb.dat")
    sp_skeb=getbin("spectra_skeb.dat")
    wvnum,wvlen=getticks("aol.txt")
    wvlen=wvlen/1000.
    line1=wvnum**(-5/3)/(10.0**2.9)
    line2=wvnum**(-3)/(10.0**9.1)
#绘制log图上斜率为-5/3和-3的直线数值
    fig, ax =plt.subplots()
    ax.loglog(wvnum,sp_noskeb,'b-',label="noskeb")
    ax.loglog(wvnum,line1,'k:')
    ax.loglog(wvnum,line2,'k:')
#绘制对数曲线图
    ax2=ax.twiny()
#双横轴必须操作
    ax2.loglog(wvlen,sp_skeb,'r--',label="skeb")
    plt.gca().invert_xaxis()
#带数据翻转横轴（横轴默认从小到大排列，而非按数组次序显示）
    ax.legend(loc='upper right', shadow=True, fontsize='medium')
    ax2.legend(loc='upper center', shadow=True, fontsize='medium')
#绘制图例
    fig.savefig('kes.eps',dpi=300)
#保存eps格式图形
    plt.show()
#在屏幕上显示图形
if __name__ == '__main__':
    main()

最终的图形如下：