In [13]:
'''
Code for plotting a few of the SBDART examples, with a function to read in select sbdart output

Written by Ben Trabing
For ATS622
contact: btrabing@colostate.edu
April 1 2019
'''

import numpy as np
from matplotlib import pyplot as plt
from pylab import *
In [14]:
def read_sbdart(filename,iout=1):
    '''
    Function to read in select output from SBDART
    Requires:
    import numpy as np
    from pylab import *
    
    Input:
    filename = file location as string
    iout = output type specified in SBDART input file
    iout = 1,10,11 are currently only options that are included
    
    Example:
    WL,FFV,TOPDN,TOPUP,TOPDIR,BOTDN,BOTUP,BOTDIR = read_sbdart('/Users/admin/sbdart/sbchk.1',iout=1)
    '''
    if iout == 1:
        '''                   
        WL,FFV,TOPDN,TOPUP,TOPDIR,BOTDN,BOTUP,BOTDIR
        
        WL    = wavelength                         (microns)
        FFV   = filter function value
        TOPDN = total downward flux at ZOUT(2) km  (w/m2/micron)
        TOPUP = total upward flux at ZOUT(2) km    (w/m2/micron)
        TOPDIR= direct downward flux at ZOUT(2) km (w/m2/micron)
        BOTDN = total downward flux at ZOUT(1) km  (w/m2/micron)
        BOTUP = total upward flux at  ZOUT(1) km   (w/m2/micron)
        BOTDIR= direct downward flux at ZOUT(1) km (w/m2/micron)
        '''
        f = loadtxt(filename,delimiter="/t",skiprows=3,dtype='str')

        WL=[]
        FFV=[]
        TOPDN=[]
        TOPUP=[]
        TOPDIR=[]
        BOTDN=[]
        BOTUP=[]
        BOTDIR=[]

        for row in f:
            c = str(row).replace('  ',' ')
            c = c.split(' ')
            
            #If multiple files are included in the same file, such as in example 3 of SBDART
            if c[0]=='"tbf' or c[1]=='':
                continue
                
            WL.append(float(c[1]))
            FFV.append(float(c[2]))
            TOPDN.append(float(c[3]))
            TOPUP.append(float(c[4]))
            TOPDIR.append(float(c[5]))
            BOTDN.append(float(c[6]))
            BOTUP.append(float(c[7]))
            BOTDIR.append(float(c[8]))
            
        return np.array(WL),np.array(FFV),np.array(TOPDN),np.array(TOPUP),np.array(TOPDIR),np.array(BOTDN),np.array(BOTUP),np.array(BOTDIR)
        
    elif iout == 10:
        '''
        WLINF,WLSUP,FFEW,TOPDN,TOPUP,TOPDIR,BOTDN,BOTUP,BOTDIR

        WLINF = lower wavelength limit             (microns)
        WLSUP = upper wavelength limit             (microns)
        FFEW  = filter function equivalent width   (microns)
        TOPDN = total downward flux at ZOUT(2) km  (w/m2)
        TOPUP = total upward flux at ZOUT(2) km    (w/m2)
        TOPDIR= direct downward flux at ZOUT(2) km (w/m2)
        BOTDN = total downward flux at ZOUT(1) km  (w/m2)
        BOTUP = total upward flux at  ZOUT(1) km   (w/m2)
        BOTDIR= direct downward flux at ZOUT(1) km (w/m2)
        '''
        f = loadtxt(filename,delimiter="/t",skiprows=0,dtype='str')

        WLINF=[]
        WLSUP=[]
        FFEW=[]
        TOPDN=[]
        TOPUP=[]
        TOPDIR=[]
        BOTDN=[]
        BOTUP=[]
        BOTDIR=[]

        for row in f:
            c = str(row).replace('   ',' ')
            c = c.replace('  ',' ')
            c = c.replace('  ',' ')
            c = c.split(' ')

            WLINF.append(float(c[1]))
            WLSUP.append(float(c[2]))
            FFEW.append(float(c[3]))            
            TOPDN.append(float(c[4]))
            TOPUP.append(float(c[5]))
            TOPDIR.append(float(c[6]))
            BOTDN.append(float(c[7]))
            BOTUP.append(float(c[8]))
            BOTDIR.append(float(c[9]))
        return np.array(WLINF),np.array(WLSUP),np.array(FFEW),np.array(TOPDN),np.array(TOPUP),\
                np.array(TOPDIR),np.array(BOTDN),np.array(BOTUP),np.array(BOTDIR)
        
        
    elif iout ==11:
        '''
        zz    = level altitudes                 (km)
        pp    = level pressure                  (mb)
        fxdn  = downward flux (direct+diffuse)  (W/m2)
        fxup  = upward flux                     (W/m2)
        fxdir = downward flux, direct beam only (W/m2)
        dfdz  = radiant energy flux divergence  (mW/m3)
        heat  = heating rate                    (K/day)
        
        '''
        f = loadtxt(filename,delimiter="/t",skiprows=1,dtype='str')

        Z=[]
        P=[]
        FXD=[]
        FXU=[]
        FXDD=[]
        DFDZ=[]
        HEAT=[]

        for row in f:
            c = str(row).replace('  ',' ')
            c = c.split(' ')
            Z.append(float(c[1]))
            P.append(float(c[2]))
            FXD.append(float(c[3]))
            FXU.append(float(c[4]))
            FXDD.append(float(c[5]))
            DFDZ.append(float(c[6]))
            HEAT.append(float(c[7]))
            
        return np.array(Z),np.array(P),np.array(FXD),np.array(FXU),np.array(FXDD),np.array(DFDZ),np.array(HEAT)

Example 1: Shortwave Fluxes¶

In [3]:
WL,FFV,TOPDN,TOPUP,TOPDIR,BOTDN,BOTUP,BOTDIR = read_sbdart('/Users/admin/Downloads/sbdart/sbchk.1',iout=1)
WL2,FFV2,TOPDN2,TOPUP2,TOPDIR2,BOTDN2,BOTUP2,BOTDIR2 = read_sbdart('/Users/admin/Downloads/sbdart/sbchk.11',iout=1)
In [4]:
fig = plt.figure(figsize=(12,8))
ax1 = plt.subplot(111)


#Slices based on resolution defined in the file for the three cases output
ax1.plot(WL,BOTDN,lw=3,color='k',ls='solid',label='Downwelling at Surface')
ax1.plot(WL,TOPDN,lw=3,color='r',ls='solid',label='Downwelling at TOA')
ax1.plot(WL2,BOTDN2,lw=3,color='gray',ls='solid',label='Downwelling at Surface')
ax1.plot(WL2,TOPDN2,lw=3,color='blue',ls='solid',label='Downwelling at TOA')

ax1.plot(WL,BOTUP,lw=3,color='k',ls='dashed',label='Upwelling at Surface')
ax1.plot(WL,TOPUP,lw=3,color='r',ls='dashed',label='Upwelling at TOA')
ax1.plot(WL2,BOTUP2,lw=3,color='gray',ls='dashed',label='Upwelling at Surface')
ax1.plot(WL2,TOPUP2,lw=3,color='blue',ls='dashed',label='Upwelling at TOA')




ax1.legend(loc=0,frameon=True,prop={'size':14})

#ax1.set_title('Example 3: Thermal IR Spectral Output',fontsize=14)
ax1.set_ylabel('Fluxes (W m$^{-2}$ $\mu$m$^{-1}$ )',fontsize=12)
ax1.set_xlabel('Wavelength ($\mu$m)',fontsize=12)

ax1.set_ylim(0)
ax1.set_xlim(.25,1)



ax1.tick_params(axis='x',labelsize=12)
ax1.tick_params(axis='y',labelsize=12)


#plt.savefig('SW.pdf',bbox_inches='tight')



plt.show()

Example 3: Thermal IR Specral Output¶

In [15]:
WL,FFV,TOPDN,TOPUP,TOPDIR,BOTDN,BOTUP,BOTDIR = read_sbdart('/Users/admin/sbdart/sbchk.3',iout=1)
In [16]:
fig = plt.figure(figsize=(12,8))
ax1 = plt.subplot(111)


#Slices based on resolution defined in the file for the three cases output
ax1.plot(WL[0:161],BOTDN[0:161],lw=3,color='C0',ls='solid',label=r'$\tau$=0')
ax1.plot(WL[161:322],BOTDN[161:322],lw=3,color='C1',ls='solid',label=r'$\tau$=1')
ax1.plot(WL[-161:],BOTDN[-161:],lw=3,color='C3',ls='solid',label=r'$\tau$=5')


ax1.legend(loc=0,frameon=True,prop={'size':14})

ax1.set_title('Example 3: Thermal IR Spectral Output',fontsize=14)
ax1.set_ylabel('Surface Downward Flux (W m$^{-2}$ $\mu$m$^{-1}$ )',fontsize=12)
ax1.set_xlabel('Wavelength ($\mu$m)',fontsize=12)

ax1.set_ylim(0)
ax1.set_xlim(4,20)



ax1.tick_params(axis='x',labelsize=12)
ax1.tick_params(axis='y',labelsize=12)


plt.savefig('Downwelling_IR_SURFACE.pdf',bbox_inches='tight')



plt.show()
In [7]:
fig = plt.figure(figsize=(12,8))
ax1 = plt.subplot(111)


#Slices based on resolution defined in the file for the three cases output
ax1.plot(WL[0:161],TOPUP[0:161],lw=3,color='C0',ls='solid',label=r'$\tau$=0')
ax1.plot(WL[161:322],TOPUP[161:322],lw=3,color='C1',ls='solid',label=r'$\tau$=1')
ax1.plot(WL[-161:],TOPUP[-161:],lw=3,color='C3',ls='solid',label=r'$\tau$=5')


ax1.legend(loc=0,frameon=True,prop={'size':14})

ax1.set_title('Example 3: Thermal IR Spectral Output',fontsize=14)
ax1.set_ylabel('TOA Upward Flux (W m$^{-2}$ $\mu$m$^{-1}$ )',fontsize=12)
ax1.set_xlabel('Wavelength ($\mu$m)',fontsize=12)

ax1.set_ylim(0)
ax1.set_xlim(4,20)



ax1.tick_params(axis='x',labelsize=12)
ax1.tick_params(axis='y',labelsize=12)


#plt.savefig('UPwelling_IR_TOA.pdf',bbox_inches='tight')



plt.show()

Example 4: Demonstration of Nakajima and King¶

In [17]:
WLINF,WLSUP,FFEW,TOPDN,TOPUP,TOPDIR,BOTDN,BOTUP,BOTDIR = read_sbdart('/Users/admin/sbdart/sbchk.4',iout=10)
In [18]:
def gtau(var,b):
    #Use this function to get the reflectances for the tau values which cannot be found through easy indexing
    tau = []
    for p in np.arange(0,7):
        tau.append(var[p::7][b])
    return np.array(tau)
In [19]:
fig = plt.figure(figsize=(12,8))
ax1 = plt.subplot(111)

#These are defined in example 4 but we can not get from the output files
reff = [2,4,8,16,32,64,128]
taus = [0,1,2,4,8,16,32,64,128]


ax1.scatter(TOPUP[::2]/TOPDN[::2],TOPUP[1::2]/TOPDN[1::2],color='k',s=100) # Just plot all the points

for p in np.arange(0,6):
    ax1.plot((TOPUP[::2]/TOPDN[::2])[p::7],(TOPUP[1::2]/TOPDN[1::2])[p::7],color='k',ls='solid',lw=4)
    ax1.text((TOPUP[::2]/TOPDN[::2])[p::7][-1]+.065, (TOPUP[1::2]/TOPDN[1::2])[p::7][-1], r'r$_{eff}$=%s $\mu$m'%reff[p],\
             horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')

for p in np.arange(0,9):
    ax1.plot(gtau((TOPUP[::2]/TOPDN[::2]),p),gtau((TOPUP[1::2]/TOPDN[1::2]),p),color='k',ls='dashed',lw=4)
    ax1.text(gtau((TOPUP[::2]/TOPDN[::2]),p)[-1], gtau((TOPUP[1::2]/TOPDN[1::2]),p)[-1]-.02, r'$\tau$=%s'%taus[p],\
             horizontalalignment='center',verticalalignment='center',fontsize=12,color='r')


ax1.set_title('Example 4: Demonstration of Nakajima and King',fontsize=14)
ax1.set_xlabel('Reflection Function ($\lambda$=.55 $\mu$m)',fontsize=14)
ax1.set_ylabel('Reflection Function ($\lambda$=2.16 $\mu$m)',fontsize=14)


ax1.set_xlim(-.05,1.1)
ax1.set_ylim(-.05,.6)


ax1.xaxis.set_ticks_position('both')
ax1.yaxis.set_ticks_position('both')


ax1.tick_params(axis='y',labelsize=12)
ax1.tick_params(axis='x',labelsize=12)


plt.savefig('Nakajima_King_Example.pdf',bbox_inches='tight')


plt.show()

Example of Heating Rate Profile¶

In [11]:
#Inclass Example
Z,P,FXD,FXU,FXDD,DFDZ,HEAT = read_sbdart('/Users/admin/sbdart/class',iout=11)
In [12]:
fig = plt.figure(figsize=(8,8))
ax1 = plt.subplot(111)

#Plot the heating Rates
ax1.plot(HEAT,Z,color='C0',lw=4)

ax1.set_xlabel('Heating Rate (K day$^{-1}$)',fontsize=12)
ax1.set_ylabel('Height (km)',fontsize=12)

ax1.set_ylim(0,20)
plt.show()
In [ ]: