In [1]:
exp_id='jra25.gndg_UR4'
exp_test='src'
model_name='miroc5-iso'
myname='kanon'
gcmdir='/data29/kanon/'+model_name+'/'
home='/home/kanon/monitor/'+model_name+'/'

# DomeF
# 南緯77度19分01秒 東経39度42分12秒座標: 南緯77度19分01秒 東経39度42分12秒
df_lat=-77.3
df_lon=39.66
c_df='black'
s_df=25

jra25.gndg on isotope3

In [2]:
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.colors as mc
from matplotlib.ticker import AutoMinorLocator
import matplotlib.path as mpath
import cartopy.crs as ccrs
from cartopy.util import add_cyclic_point
import cartopy.feature as feature
from cmocean import cm as cmo
from netCDF4 import Dataset as NetCDFFile
import statsmodels.api as sm
import statsmodels.tools.eval_measures as smte
import os
import datetime
import paramiko
import pandas as pd
import datetime
import numpy as np
import matplotlib.dates as mdates
from cdo import *
cdo =Cdo()
%matplotlib inline
In [3]:
latmin=-90
latmax=-40
nlats=18
lonmin=-180
lonmax=180
nlons = 128 # 経度方向のデータ数
year=2003 
numofdays=365
In [4]:
# Compute a circle in axes coordinates, which we can use as a boundary
# for the map. We can pan/zoom as much as we like - the boundary will be
# permanently circular.
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
In [5]:
# For the size of the figures
plt.rcParams['figure.dpi'] = 300

Current date and time

In [6]:
now = datetime.datetime.now()
print (now.strftime("%Y-%m-%d %H:%M"))
2020-08-09 00:11

The following jobs are currently running on isotope3

In [7]:
client = paramiko.SSHClient()
client.load_system_host_keys()
client.connect('isotope3.iis.u-tokyo.ac.jp', username='kanon')

stdin, stdout, stderr = client.exec_command('qstat')
queue_status = stdout.readlines()
queue_status = [l.split() for l in queue_status]
def check_everything():
    if len(queue_status) != 0:
        queue_df = pd.DataFrame(queue_status[2:])
        queue_df.columns = ['JOB ID', 'JOB NAME', 'USERNAME', 'CPU TIME USE', 'STATUS', 'QUEUE']
        if any(queue_df['USERNAME'].str.startswith('kanon')):
            if any(queue_df['USERNAME'].str.startswith('kanon') & queue_df['JOB NAME'].str.startswith(exp_id[0:4])):
                return queue_df[(queue_df['USERNAME'] == 'kanon') & (queue_df['JOB NAME'] == exp_id)]
            else:
                print ('No '+ exp_id + ' job running on isotope3')
        else:
            print("no jobs running on isotope3")
        return None
check_everything()
no jobs running on isotope3
In [8]:
# Function to read the log file
def get_log_output(exp):
    client = paramiko.SSHClient()
    client.load_system_host_keys()
    client.connect('isotope3.iis.u-tokyo.ac.jp', username='kanon')
    stdin, stdout, stderr = client.exec_command("cd "+gcmdir+"out/"+exp+ "; cat y*/SYSOUT")
    return stdout

# Function to extract the date/time of the simulations (per model year)
def generate_dataframe_from_cat_SYSOUT(log):
    # Read the output of the SYSOUT files 
    log_dataframe = pd.read_table(log, 
                                  sep=r" at ", 
                                  infer_datetime_format=True, 
                                  names=["State", "Date"], 
                                  engine='python')
    columnsTitles = ["Date", "State"]
    log_dataframe = log_dataframe.reindex(columns=columnsTitles)
    # Add the run number to log_dataframe
    run_number = np.array(log_dataframe["State"])
    x=0
    for i in range(len(log_dataframe["State"])):
        if log_dataframe["State"][i].endswith('started'):
            x=x+1
            run_number[i]=x
        else:
            run_number[i]=x
    log_dataframe['Run Number'] = run_number
    # Set the date of the run as an index
    log_dataframe.set_index("Date", inplace = True)
    # Convert the date to be read by python
#   log_dataframe.set_index(pd.to_datetime(log_dataframe.index, utc=True), inplace=True)
    log_dataframe.index("2015-11-05")
#    log_dataframe.index = log_dataframe.index.tz_convert('Asia/Tokyo')

    return log_dataframe

# Function to calculate the mean time that it takes to perform one model year 
def compute_effective_throughput(log_dataframe):
    starts = log_dataframe[log_dataframe.State == "job started"]
    ends = log_dataframe[log_dataframe.State == "job end"]
    # Drop the duplicated starts:
    starts.drop_duplicates(subset="Run Number", keep="last")
    merged = pd.concat([starts, ends])
    groupby = merged.groupby("Run Number")
    run_diffs = {"Run Number": [], "Time Diff": []}
    for name, group in groupby:
        run_diffs["Run Number"].append(int(name))
        run_diffs["Time Diff"].append(group.index[-1] - group.index[0])
    diffs = pd.DataFrame(run_diffs).sort_values("Run Number").set_index("Run Number")
    DAY = datetime.timedelta(1)
    throughput = DAY / diffs.mean()
    return diffs.mean(), throughput, diffs

#walltime, throughput, diffs = compute_effective_throughput(
#    generate_dataframe_from_cat_SYSOUT(get_log_output(exp_id))
#df = pd.DataFrame.from_dict({"Walltime": walltime, "Throughput": throughput}, orient="index")
#df.columns = [exp_id]
#df
In [9]:
# Simulation Timeline
#log_dataframe = generate_dataframe_from_cat_SYSOUT(get_log_output(exp_id))
# Drop the last entry if it's start:
#if "job started" in log_dataframe.iloc[-1]["State"]:
 #   end_of_log = log_dataframe.iloc[:-1].tail(30)
#else:
#    end_of_log = log_dataframe.tail(30)
#end_groups = end_of_log.groupby("Run Number")
#f, ax = plt.subplots(1, 1, dpi=150, figsize=(15, 1.5))
#for name, group in end_groups:
#    bdate = group.index[0]
#    edate = group.index[1]
#    edate, bdate = [mdates.date2num(item) for item in (edate, bdate)]
    # The color is the same as the progressbar below, use the colormeter to figure it out.
#    ax.barh(0, edate - bdate, left=bdate, height=0.2, color=(217./255., 83./255., 79./255.),
 #           edgecolor="black")
#ax.set_ylim(-0.5, 0.5)
#for direction in ["top", "left", "right"]:
 #   ax.spines[direction].set_visible(False)
#ax.yaxis.set_visible(False)
#ax.xaxis_date()
#ax.xaxis.set_major_locator(mdates.HourLocator(interval=4))
#ax.xaxis.set_minor_locator(mdates.HourLocator(interval=1))
#ax.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M %d.%m.%y", tz=log_dataframe.index.tz))
In [10]:
# Get the start, end and current years
stdin, stdout, stderr = client.exec_command("cd "+gcmdir+"sh/"+exp_id+ "; cat agcm_"+exp_id+"_*.sh | grep 'setenv YSTR'")
#stdin, stdout, stderr = client.exec_command("cd "+gcmdir+"sh/"+exp_id+ "; grep 'setenv YSTR' "+exp_id+".e*")
#start_year = int(stdout.readlines()[0].split()[2])
start_year=1981
stdin, stdout, stderr = client.exec_command("cd "+gcmdir+"sh/"+exp_id+ "; cat agcm_"+exp_id+"_*.sh | grep 'setenv YINT'")
#stdin, stdout, stderr = client.exec_command("cd "+gcmdir+"sh/"+exp_id+ "; grep 'setenv YINT' "+exp_id+".e*")
#int_year= int(stdout.readlines()[0].split()[2])
int_year=30
#final_year = start_year + int_tear
stdin, stdout, stderr = client.exec_command("cd "+gcmdir+"sh/"+exp_id+ "; grep '@ YEAR =' "+model_name+".e* | tail -n 1")# | sed s/y//")
#current_year = int(stdout.readlines()[0].split()[3])
#current_year = int(current_year) - start_year
current_year=int(2010)-start_year
final_year = int_year - start_year

# Run finished?
# If the simulation is finished, stdout.readlines() = ['exit\n'], then its length is equal to 1. If not, its
# length is equal to 0.
stdin, stdout, stderr = client.exec_command("cd "+gcmdir+"sh/"+exp_id+ "; grep exit "+exp_id+".e*")
run_finished = len(stdout.readlines())

# Progress Bar
#years_per_day = np.round(df[exp_id].Throughput)
#years_left = final_year - (current_year + run_finished)
#days_left = years_left / years_per_day
#finishing_date = datetime.datetime.now() + datetime.timedelta(days=days_left)
#desc = "Done on "+finishing_date.strftime("%d %b, %Y")+ ": {percentage:3.0f}%"
#r_bar = " "+str(current_year+run_finished)+"/"+str(final_year)+ ", Throughput ~"+str(years_per_day)+"yrs/day"
#from tqdm import tqdm_notebook
#pbar = tqdm_notebook(total=final_year, desc=desc, bar_format=desc+"{bar}"+r_bar)
#pbar.update(final_year - years_left)
#pbar.close()
In [11]:
temp2_file = home +exp_id +'/' +exp_id +'_temp2_miroc5-iso_catted_yearmean_fldmean.nc'
temp2_dataset = xr.open_dataset(temp2_file, decode_times=False)
temp2_rm = temp2_dataset.rolling(time=15, center=True).mean()
d18Op_file = home +exp_id +'/' +exp_id +'_d18Op_miroc5-iso_catted_yearmean_fldmean.nc'
d18Op_dataset = xr.open_dataset(d18Op_file, decode_times=False)
d18Op_rm = d18Op_dataset.rolling(time=15, center=True).mean()

Monitoring of temperature and $\delta$$^{18}O$ in precipitation

In [12]:
f, (ax1, ax2) = plt.subplots(1, 2, dpi=150, figsize=(8,4))

ax1.plot(temp2_dataset.tas.data.squeeze())
ax1.plot(temp2_rm.tas.data.squeeze(), 'r')
ax1.set_xlabel("Simulation Time (Years 1981-2010)")
ax1.set_ylabel('Global $\mathsf{T_{2m}}$' + ' [\u00b0C]')

ax2.plot(d18Op_dataset.d18O_precip.data.squeeze())
ax2.plot(d18Op_rm.d18O_precip.data.squeeze(), 'r')
ax2.set_xlabel("Simulation Time (Years 1981-2010)")
ax2.set_ylabel('Global $\mathsf{\delta^{18}O_{p}}$' + ' [\u2030]')

f.suptitle(exp_id +" - miroc5-iso", fontsize=14)

plt.tight_layout()
plt.subplots_adjust(top=0.85)
In [13]:
# Iso-Miroc5
ifile_PI_timmean_atm = home +exp_id +'/' +exp_id +'_miroc5-iso_last10yrs_timmean.iso.nc'
ifile_PI_yseasmean_atm = home +exp_id +'/' +exp_id +'_miroc5-iso_last10yrs_yseasmean.iso.nc'
nc_PI_timmean_atm = NetCDFFile(ifile_PI_timmean_atm)
nc_PI_yseasmean_atm = NetCDFFile(ifile_PI_yseasmean_atm)

# Coordinates
lon_miroc = nc_PI_timmean_atm.variables['lon']
lat_miroc = nc_PI_timmean_atm.variables['lat']
# Time-mean values
d18Op_PI = nc_PI_timmean_atm.variables['d18O_precip'][0,:,:]
temp2m_PI = nc_PI_timmean_atm.variables['tas'][0,:,:]
precip_PI = nc_PI_timmean_atm.variables['pr'][0,:,:]
# Multiyear seasonal mean values
d18Op_yseasmean_PI = nc_PI_yseasmean_atm.variables['d18O_precip'][:,:,:]
temp2m_yseasmean_PI = nc_PI_yseasmean_atm.variables['tas'][:,:,:]
precip_yseasmean_PI = nc_PI_yseasmean_atm.variables['pr'][:,:,:]

# Time-mean values
d18Op_PI_cyc, lon_miroc_cyc = add_cyclic_point(d18Op_PI, coord=lon_miroc)
temp2m_PI_cyc, lon_miroc_cyc = add_cyclic_point(temp2m_PI, coord=lon_miroc)
precip_PI_cyc, lon_miroc_cyc = add_cyclic_point(precip_PI, coord=lon_miroc)
# Multiyear seasonal mean values
d18Op_yseasmean_PI_cyc, lon_miroc_cyc = add_cyclic_point(d18Op_yseasmean_PI, coord=lon_miroc)
temp2m_yseasmean_PI_cyc, lon_miroc_cyc = add_cyclic_point(temp2m_yseasmean_PI, coord=lon_miroc)
precip_yseasmean_PI_cyc, lon_miroc_cyc = add_cyclic_point(precip_yseasmean_PI, coord=lon_miroc)
In [14]:
# Data
ignip = home +exp_test +'/' +'GNIP_world_2010_timemean_TPd18oD_12mon.csv'
iice = home +exp_test +'/' +'IceCores_6k_PI.csv'
ispelothem = home +exp_test +'/' +'SISALv1b_d18O_carb_MH_anom_20190402.csv'
era40 = home +exp_test +'/' +'ERA40.1961.1990.t2m.stl1.stl4.timmean.nc'

# GNIP data
lat_gnip, lon_gnip, nyr_gnip, d18O_gnip, dD_gnip, prec_gnip, temp_gnip = np.genfromtxt(ignip, unpack=True, skip_header=1, delimiter=';', usecols = (1,2,4,6,7,8,9))
for i in range(len(lat_gnip)):
    if nyr_gnip[i] < 5:
        lat_gnip[i]  = np.nan
        lon_gnip[i]  = np.nan
        d18O_gnip[i] = np.nan 
        dD_gnip[i]   = np.nan
        prec_gnip[i] = np.nan 
        temp_gnip[i] = np.nan
dex_gnip = dD_gnip - 8.*d18O_gnip

# Ice core data
lat_ice, lon_ice, pd_tsurf_ice, pd_accu_ice, pd_d18O_ice, pd_dex_ice, delta_6k_pd_d18O_ice, delta_6k_pd_dex_ice = np.genfromtxt(iice, unpack=True, skip_header=1, delimiter=';', usecols = (1,2,3,4,5,6,7,8))

# Speleothems data
lat_speleo, lon_speleo, d18O_pi_pdb_speleo, d18O_6k_pdb_speleo, delta_6k_pi_d18O_pdb_speleo = np.genfromtxt(ispelothem, unpack=True, skip_header=1, delimiter=',', usecols = (1,2,5,6,7))
d18O_pi_csmow_speleo = (d18O_pi_pdb_speleo + 29.98) / 0.97002 #PDB to SMOW scale
temp_era40 = np.empty(len(lat_speleo))
for i in range(len(lat_speleo)):
    t = cdo.selvar('stl1', input="-remapbil,lon="+np.str(lon_speleo[i])+"_lat="+np.str(lat_speleo[i])+" " + era40, options = "-f nc", returnArray  =  'stl1')
    temp_era40[i] = t[:,0,0]+273.15 # temperature needs to be in Kelvin
d18O_pi_vsmow_speleo = np.empty(len(lat_speleo))
for i in range(len(lat_speleo)):
    d18O_pi_vsmow_speleo[i] = d18O_pi_csmow_speleo[i] - ((16.1*1000.)/temp_era40[i] - 24.6) # Tremaine et al. 2011

Maps of 2m-temperature

We take 10 years (1996-2005) of simulation to make to the annual mean maps.

In [15]:
fig, axes = plt.subplots(1, 3, figsize=(13,3.5), subplot_kw=dict(projection=ccrs.Robinson()))
plt.rcParams.update({'font.family': 'sans-serif', 'text.usetex': False})

levels = [-50, -40, -30, -20, -15, -10, -5, -2, 0, 2, 5, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 30]
n_colors = len(levels) + 1
cmap = plt.get_cmap('RdYlBu_r', n_colors)
colors = cmap(range(cmap.N))
cmap, norm = mc.from_levels_and_colors(levels, colors, extend='both')

cs_world1 = axes[0].pcolormesh(lon_miroc_cyc, lat_miroc, temp2m_PI_cyc, norm=norm, cmap=cmap, transform=ccrs.PlateCarree())
axes[0].set_title('Annual mean', fontsize=16)
axes[0].coastlines()
cs_world2 = axes[1].pcolormesh(lon_miroc_cyc, lat_miroc, temp2m_yseasmean_PI_cyc[0,:,:], norm=norm, cmap=cmap, transform=ccrs.PlateCarree())
axes[1].set_title('DJF', fontsize=16)
axes[1].coastlines()
cs_world3 = axes[2].pcolormesh(lon_miroc_cyc, lat_miroc, temp2m_yseasmean_PI_cyc[2,:,:], norm=norm, cmap=cmap, transform=ccrs.PlateCarree())
axes[2].set_title('JJA', fontsize=16)
axes[2].coastlines()

cbar_ax = fig.add_axes([0.02, 0.2, 0.96, 0.05])
cbar=fig.colorbar(cs_world1, cax=cbar_ax, orientation="horizontal", ticks=levels)
cbar.ax.tick_params(labelsize=10)

cbar.set_label('$\mathsf{T_{2m}}$' + ' [\u00b0C]', size=18)
plt.subplots_adjust(top=0.92,
bottom=0.29,
left=0.02,
right=0.98,
hspace=0.2,
wspace=0.05)

plt.show()