# Read sea ice area in the Arctic
#  - Hindcast/forecast simulations
#  - sea ice area calculated from OISST
# Calibrate forecast
#  - Climatology (1993-2016)
#  - Compute the climatological differences between simulations and observations
# Plot observations and simulated sea ice area with forecast calibration 
#
#
# Script version 2.1
#
# Python version: 3.10
#
# Package version
#   numpy:  1.26.3
#   pandas:  2.2.0
#   matplotlib:  3.7.3
#
# Ver1.0.0: Created by Mariko Koseki, 20.06.2024
# Ver1.1.0: Updated by Mariko, 02.08.2024
# Ver2.0.0: Updated by Mariko, 15.08.2024
# ver2.1: updated by Mariko, 08.12.2025



import numpy as np
import pandas as pd
import datetime
from dateutil.relativedelta import relativedelta
import statistics
import math
import collections
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import matplotlib.ticker as ptick
import os
import glob
import sys
import platform
import calendar


print('--- Python version ---')
print(platform.python_version())

print('--- Package version ---')
print('numpy: ', np.__version__)
print('pandas: ', pd.__version__)
print('matplotlib: ', matplotlib.__version__)


##--- Set path 1 ----------------------------------##
## Set path to input files ##
'''
Set path to input directory ###
'''
### Workspace for sea ice forecast
seaice_dir = '/nird/datapeak/NS9873K/norcpm/validation/seaicepre/' # ver2.1

### Hindcast/forecast
out_hind = seaice_dir + 'hindcast/'
### Observation from OISST (merged)
out_oisst_post = seaice_dir + 'oisst/merge/'

### Climatology (1993-2016)
out_clm_mdl = seaice_dir + 'clm/mdl/clm_1993_2016/' #ver2.0.0'
out_clm_oisst = seaice_dir + 'clm/oisst/'


## Select NorCPM system ##
'''
Set system number 'system_num'
'''
system_num = '1'
#system_num = '2'

print('')
print('system: ', system_num)
print('')



##--- Input -----------------------------------------##
args = sys.argv
if len(args) == 2:
    input_date = args[1]
    if len(input_date) == 6:
        yyyy = int(input_date[0:4])
        mm = int(input_date[4:6])
        dd = 1

    else:
        print('')
        print('---How to use this script---')
        print('python plot_seaice_oisst_mdl.py <yyyymm>')
        print('<yyyymm> = Forecast start date')
        print('Example: python plot_seaice_oisst_mdl.py 202312')
        print('')
        sys.exit()

else:
    print('')
    print('---How to use this script---')
    print('python plot_seaice_oisst_mdl.py <yyyymm>')
    print('<yyyymm> = Forecast start date')
    print('Example: python plot_seaice_oisst_mdl.py 202312')
    print('')
    sys.exit()



if mm >= 8:
    #print(mm)
    start_month_mdl = mm
    start_year_mdl = yyyy
    end_month_mdl = mm - 7
    end_year_mdl = yyyy + 1
    init_month_mdl = mm - 1
    init_year_mdl = yyyy

    start_month_obs = mm - 6
    start_year_obs = yyyy
    end_month_obs = mm - 1
    end_year_obs = yyyy
else:
    start_month_mdl = mm
    start_year_mdl = yyyy
    end_month_mdl = mm + 5
    end_year_mdl = yyyy

    if start_month_mdl == 1:
        init_month_mdl = 12
        init_year_mdl = yyyy - 1
    
        start_month_obs = mm + 6
        start_year_obs = yyyy - 1
        end_month_obs = 12
        end_year_obs = yyyy - 1
    elif start_month_mdl == 7:
        init_month_mdl = mm - 1
        init_year_mdl = yyyy

        start_month_obs = mm - 6
        start_year_obs = yyyy
        end_month_obs = mm - 1
        end_year_obs = yyyy

    else:
        init_month_mdl = mm - 1
        init_year_mdl = yyyy
        
        start_month_obs = mm + 6
        start_year_obs = yyyy - 1
        end_month_obs = mm - 1
        end_year_obs = yyyy

recent_year_obs = str(yyyy); last_year_obs = str(yyyy - 1)
year_obs_1 = str(yyyy - 1); year_obs_2 = str(yyyy - 2)
year_obs_3 = str(yyyy - 3); year_obs_4 = str(yyyy - 4)
year_obs_5 = str(yyyy - 5); year_obs_6 = str(yyyy - 6)
year_obs_9 = str(yyyy - 9); year_obs_10 = str(yyyy - 10); year_obs_11 = str(yyyy - 11)
year_obs_19 = str(yyyy - 19); year_obs_20 = str(yyyy - 20); year_obs_21 = str(yyyy - 21)
year_obs_29 = str(yyyy - 29); year_obs_30 = str(yyyy - 30); year_obs_31 = str(yyyy - 31)



start_year = '{:0>4d}'.format(yyyy)
start_month = '{:0>2d}'.format(mm)
start_date = '{:0>2d}'.format(dd)

issue_date = start_year + start_month + start_date
print('')
print('Issued: ', issue_date)
print('')


forecast_start = '{:0>2d}/{:0>2d}/{:0>4d}'.format(dd,mm,yyyy)
print('Forecast start: ', forecast_start)

forecast_start_date = '{:0>4d}{:0>2d}'.format(start_year_mdl,start_month_mdl) #string
forecast_end_date = '{:0>4d}{:0>2d}'.format(end_year_mdl,end_month_mdl) #string
init_date = '{:0>4d}{:0>2d}'.format(init_year_mdl,init_month_mdl) #string

obs_start_date = '{:0>4d}{:0>2d}'.format(start_year_obs,start_month_obs) #string
obs_end_date = '{:0>4d}{:0>2d}'.format(end_year_obs,end_month_obs) #string

print('')
print('Forecast start date: ', forecast_start_date) #Forecast start date
print('Forecast end date: ', forecast_end_date) #Forecast end date

print('simulation start date: ', init_date) #Start date of simulation/forecast initialization date

print('Observation start date: ', obs_start_date) #Observation start date
print('Observation end date: ', obs_end_date) #Observation end date
print('')
#sys.exit()


##--- Set path 2 ----------------------------##
### Input: Set path to input files ###
'''
Set path to input files ###
'''
## Hindcast/forecast ##
file_hind = out_hind + 'norcpm-cf-system1_hindcast1_' + init_date + '15_mem01-60/seaice_hind_' + init_date + '.csv'
print(file_hind)


## Observation from OISST (merged) ##
### Latest data ###
file_obs_post_last = out_oisst_post + 'seaice_oisst_merge_' + last_year_obs + '.csv'
file_obs_post_recent = out_oisst_post + 'seaice_oisst_merge_' + recent_year_obs + '.csv'

### Historical data ###
file_obs_post_1 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_1 + '.csv'; file_obs_post_2 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_2 + '.csv'
file_obs_post_3 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_3 + '.csv'; file_obs_post_4 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_4 + '.csv'
file_obs_post_5 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_5 + '.csv'; file_obs_post_6 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_6 + '.csv'
file_obs_post_9 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_9 + '.csv'; file_obs_post_10 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_10 + '.csv'
file_obs_post_11 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_11 + '.csv'
file_obs_post_19 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_19 + '.csv'
file_obs_post_20 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_20 + '.csv'
file_obs_post_21 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_21 + '.csv'
file_obs_post_29 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_29 + '.csv'
file_obs_post_30 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_30 + '.csv'
file_obs_post_31 = out_oisst_post + 'seaice_oisst_merge_' + year_obs_31 + '.csv'


print(file_obs_post_last)
print(file_obs_post_recent)
#sys.exit()


## Climatology ##
file_clm_mdl = out_clm_mdl + 'mdl_area_' + str(start_month_mdl).zfill(2) + '_' + str(end_month_mdl).zfill(2) + '_clm.csv' #ver2.0.0
file_clm_oisst = out_clm_oisst + 'oisst_seaice_area_clm_' + str(start_month_mdl).zfill(2) + '_' + str(end_month_mdl).zfill(2) + '.csv'
print(file_clm_mdl)
print(file_clm_oisst)
#sys.exit()


### Output: Set path to output file ###
'''
Set path to output file
'''
#out_dir = seaice_dir + 'fig/' # <--- tmp forlder
out_dir  = '/nird/datapeak/NS9873K/www/norcpm/forecast/plots/version2/' # <--- folder for archiving figures # ver2.1

if not os.path.exists(out_dir):
    os.makedirs(out_dir)



##--- Read CSV files into pandas.DataFrame ------------------------------##
## Hindcast/forecast ##
ice_hind = pd.read_csv(file_hind)


## Observation from OISST (merged) ##
### Latest data ###
ice_obs_post_last = pd.read_csv(file_obs_post_last); ice_obs_post_recent = pd.read_csv(file_obs_post_recent)

### Historical data ###
if os.path.exists(file_obs_post_2):
    ice_obs_post_1 = pd.read_csv(file_obs_post_1); ice_obs_post_2 = pd.read_csv(file_obs_post_2)

if os.path.exists(file_obs_post_4):
    ice_obs_post_3 = pd.read_csv(file_obs_post_3); ice_obs_post_4 = pd.read_csv(file_obs_post_4)

if os.path.exists(file_obs_post_6):
    ice_obs_post_5 = pd.read_csv(file_obs_post_5); ice_obs_post_6 = pd.read_csv(file_obs_post_6)

if os.path.exists(file_obs_post_11):
    ice_obs_post_9 = pd.read_csv(file_obs_post_9)
    ice_obs_post_10 = pd.read_csv(file_obs_post_10)
    ice_obs_post_11 = pd.read_csv(file_obs_post_11)

if os.path.exists(file_obs_post_21):
    ice_obs_post_19 = pd.read_csv(file_obs_post_19)
    ice_obs_post_20 = pd.read_csv(file_obs_post_20)
    ice_obs_post_21 = pd.read_csv(file_obs_post_21)

if os.path.exists(file_obs_post_31):
    ice_obs_post_29 = pd.read_csv(file_obs_post_29)
    ice_obs_post_30 = pd.read_csv(file_obs_post_30)
    ice_obs_post_31 = pd.read_csv(file_obs_post_31)





## Climatology
ice_clm_mdl = pd.read_csv(file_clm_mdl)
ice_clm_oisst = pd.read_csv(file_clm_oisst)





##--- Forecast calibration -----------------------------------------------##
## Change 'object' to 'datetime' ##
ice_hind['dtime'] = pd.to_datetime(ice_hind['dtime'])

## Extract simulation period ##
ice_hind = ice_hind[(ice_hind["dtime"] >= datetime.datetime(start_year_mdl, start_month_mdl, 1)) \
        & (ice_hind["dtime"] <= datetime.datetime(end_year_mdl, end_month_mdl, calendar.monthrange(end_year_mdl, end_month_mdl)[1]))]


## Insert NaN if there is missing data ## #Ver1.1.0
all_date = pd.DataFrame(pd.date_range(start = datetime.datetime(start_year_mdl, start_month_mdl, 1), end = datetime.datetime(end_year_mdl, end_month_mdl, calendar.monthrange(end_year_mdl, end_month_mdl)[1]), freq="D"), columns=["dtime"])
#print(all_date)
ice_hind_fill = pd.merge(all_date, ice_hind, on = 'dtime', how = 'left')
#print(ice_hind_fill)


## Reset index ##
ice_hind_res = ice_hind_fill.reset_index()
#print(ice_hind)
#print(ice_hind_res)


## Extract variables ##
time_hind = ice_hind_res['dtime']

ensemble_mean = ice_hind_res['ensemble_mean']
#print(ensemble_mean)
std_ens_ice = ice_hind_res['std_ens_ice']
percentile_10 = ice_hind_res['percentile_10']
percentile_90 = ice_hind_res['percentile_90']
#print(percentile_10)


## Compute climatological difference ##
clm_mdl = ice_clm_mdl['mean']
clm_oisst = ice_clm_oisst['mean']
clm_diff = clm_mdl - clm_oisst

#print(clm_diff)


## Calibrate forecast ##
ensemble_mean_calib = ensemble_mean - clm_diff
percentile_10_calib = percentile_10 - clm_diff
percentile_90_calib = percentile_90 - clm_diff
#print(ensemble_mean_calib)
#print(percentile_10_calib)





##--- Read observation -----------------------------------------------------##
## Change 'object' to 'datetime' ##
### Latest data ###
ice_obs_post_recent['dtime'] = pd.to_datetime(ice_obs_post_recent['dtime'])
ice_obs_post_last['dtime'] = pd.to_datetime(ice_obs_post_last['dtime'])

### Historical data ###
if os.path.exists(file_obs_post_2):
    ice_obs_post_1['dtime'] = pd.to_datetime(ice_obs_post_1['dtime']); ice_obs_post_2['dtime'] = pd.to_datetime(ice_obs_post_2['dtime'])
if os.path.exists(file_obs_post_4):
    ice_obs_post_3['dtime'] = pd.to_datetime(ice_obs_post_3['dtime']); ice_obs_post_4['dtime'] = pd.to_datetime(ice_obs_post_4['dtime'])
if os.path.exists(file_obs_post_6):
    ice_obs_post_5['dtime'] = pd.to_datetime(ice_obs_post_5['dtime']); ice_obs_post_6['dtime'] = pd.to_datetime(ice_obs_post_6['dtime'])
if os.path.exists(file_obs_post_11):
    ice_obs_post_9['dtime'] = pd.to_datetime(ice_obs_post_9['dtime']); ice_obs_post_10['dtime'] = pd.to_datetime(ice_obs_post_10['dtime'])
    ice_obs_post_11['dtime'] = pd.to_datetime(ice_obs_post_11['dtime'])
if os.path.exists(file_obs_post_21):
    ice_obs_post_19['dtime'] = pd.to_datetime(ice_obs_post_19['dtime'])
    ice_obs_post_20['dtime'] = pd.to_datetime(ice_obs_post_20['dtime']); ice_obs_post_21['dtime'] = pd.to_datetime(ice_obs_post_21['dtime'])
if os.path.exists(file_obs_post_31):
    ice_obs_post_29['dtime'] = pd.to_datetime(ice_obs_post_29['dtime'])
    ice_obs_post_30['dtime'] = pd.to_datetime(ice_obs_post_30['dtime']); ice_obs_post_31['dtime'] = pd.to_datetime(ice_obs_post_31['dtime'])



## Combine two/three DataFrame objects ##
### Latest data ###
ice_obs_post_mer = pd.concat([ice_obs_post_last, ice_obs_post_recent], axis=0)

### Historical data ###
if os.path.exists(file_obs_post_2):
    ice_obs_post_mer1 = pd.concat([ice_obs_post_2, ice_obs_post_1, ice_obs_post_recent], axis=0)
if os.path.exists(file_obs_post_4):
    ice_obs_post_mer3 = pd.concat([ice_obs_post_4, ice_obs_post_3, ice_obs_post_2], axis=0)
if os.path.exists(file_obs_post_6):
    ice_obs_post_mer5 = pd.concat([ice_obs_post_6, ice_obs_post_5, ice_obs_post_4], axis=0)
if os.path.exists(file_obs_post_11):
    ice_obs_post_mer10 = pd.concat([ice_obs_post_11, ice_obs_post_10, ice_obs_post_9], axis=0)
if os.path.exists(file_obs_post_21):
    ice_obs_post_mer20 = pd.concat([ice_obs_post_21, ice_obs_post_20, ice_obs_post_19], axis=0)
if os.path.exists(file_obs_post_31):
    ice_obs_post_mer30 = pd.concat([ice_obs_post_31, ice_obs_post_30, ice_obs_post_29], axis=0)





## Extract simulation period ##
### Latest data ###
ice_obs_post_mer = ice_obs_post_mer[(ice_obs_post_mer["dtime"] >= datetime.datetime(start_year_obs, start_month_obs, 1)) \
        & (ice_obs_post_mer["dtime"] <= datetime.datetime(end_year_obs, end_month_obs, calendar.monthrange(end_year_obs, end_month_obs)[1]))]

### Historical data ###
if os.path.exists(file_obs_post_2):
    ice_obs_post_mer1 = ice_obs_post_mer1[(ice_obs_post_mer1["dtime"] >= datetime.datetime(start_year_obs - 1, start_month_obs, 1)) \
        & (ice_obs_post_mer1["dtime"] <= datetime.datetime(end_year_mdl - 1, end_month_mdl, calendar.monthrange(end_year_mdl - 1, end_month_mdl)[1]))]
if os.path.exists(file_obs_post_4):
    ice_obs_post_mer3 = ice_obs_post_mer3[(ice_obs_post_mer3["dtime"] >= datetime.datetime(start_year_obs - 3, start_month_obs, 1)) \
        & (ice_obs_post_mer3["dtime"] <= datetime.datetime(end_year_mdl - 3, end_month_mdl, calendar.monthrange(end_year_mdl - 3, end_month_mdl)[1]))]
if os.path.exists(file_obs_post_6):
    ice_obs_post_mer5 = ice_obs_post_mer5[(ice_obs_post_mer5["dtime"] >= datetime.datetime(start_year_obs - 5, start_month_obs, 1)) \
        & (ice_obs_post_mer5["dtime"] <= datetime.datetime(end_year_mdl - 5, end_month_mdl, calendar.monthrange(end_year_mdl - 5, end_month_mdl)[1]))]
if os.path.exists(file_obs_post_11):
    ice_obs_post_mer10 = ice_obs_post_mer10[(ice_obs_post_mer10["dtime"] >= datetime.datetime(start_year_obs - 10, start_month_obs, 1)) \
        & (ice_obs_post_mer10["dtime"] <= datetime.datetime(end_year_mdl - 10, end_month_mdl, calendar.monthrange(end_year_mdl - 10, end_month_mdl)[1]))]
if os.path.exists(file_obs_post_21):
    ice_obs_post_mer20 = ice_obs_post_mer20[(ice_obs_post_mer20["dtime"] >= datetime.datetime(start_year_obs - 20, start_month_obs, 1)) \
        & (ice_obs_post_mer20["dtime"] <= datetime.datetime(end_year_mdl - 20, end_month_mdl, calendar.monthrange(end_year_mdl - 20, end_month_mdl)[1]))]
if os.path.exists(file_obs_post_31):
    ice_obs_post_mer30 = ice_obs_post_mer30[(ice_obs_post_mer30["dtime"] >= datetime.datetime(start_year_obs - 30, start_month_obs, 1)) \
        & (ice_obs_post_mer30["dtime"] <= datetime.datetime(end_year_mdl - 30, end_month_mdl, calendar.monthrange(end_year_mdl - 30, end_month_mdl)[1]))]


## Reset index ##
### Latest data ###
ice_obs_post_mer_res = ice_obs_post_mer.reset_index()

### Historical data ###
if os.path.exists(file_obs_post_2):
    ice_obs_post_mer_res1 = ice_obs_post_mer1.reset_index()
if os.path.exists(file_obs_post_4):
    ice_obs_post_mer_res3 = ice_obs_post_mer3.reset_index()
if os.path.exists(file_obs_post_6):
    ice_obs_post_mer_res5 = ice_obs_post_mer5.reset_index()
if os.path.exists(file_obs_post_11):
    ice_obs_post_mer_res10 = ice_obs_post_mer10.reset_index()
if os.path.exists(file_obs_post_21):
    ice_obs_post_mer_res20 = ice_obs_post_mer20.reset_index()
if os.path.exists(file_obs_post_31):
    ice_obs_post_mer_res30 = ice_obs_post_mer30.reset_index()



## Extract variables ##
### Latest data ###
time_obs =ice_obs_post_mer_res['dtime']
area_inter_obs = ice_obs_post_mer_res['sea_ice_area_interpolated']

### Historical data ###
if os.path.exists(file_obs_post_2):
    ice_obs_post_mer_res1["DATE_SHIFT"] = ice_obs_post_mer_res1["dtime"].apply(lambda x: x + relativedelta(years=1))
    time_obs1 =ice_obs_post_mer_res1['DATE_SHIFT']
    area_inter_obs1 = ice_obs_post_mer_res1['sea_ice_area_interpolated']
if os.path.exists(file_obs_post_4):
    ice_obs_post_mer_res3["DATE_SHIFT"] = ice_obs_post_mer_res3["dtime"].apply(lambda x: x + relativedelta(years=3))
    time_obs3 =ice_obs_post_mer_res3['DATE_SHIFT']
    area_inter_obs3 = ice_obs_post_mer_res3['sea_ice_area_interpolated']
if os.path.exists(file_obs_post_6):
    ice_obs_post_mer_res5["DATE_SHIFT"] = ice_obs_post_mer_res5["dtime"].apply(lambda x: x + relativedelta(years=5))
    time_obs5 =ice_obs_post_mer_res5['DATE_SHIFT']
    area_inter_obs5 = ice_obs_post_mer_res5['sea_ice_area_interpolated']
if os.path.exists(file_obs_post_11):
    ice_obs_post_mer_res10["DATE_SHIFT"] = ice_obs_post_mer_res10["dtime"].apply(lambda x: x + relativedelta(years=10))
    time_obs10 =ice_obs_post_mer_res10['DATE_SHIFT']
    area_inter_obs10 = ice_obs_post_mer_res10['sea_ice_area_interpolated']
if os.path.exists(file_obs_post_21):
    ice_obs_post_mer_res20["DATE_SHIFT"] = ice_obs_post_mer_res20["dtime"].apply(lambda x: x + relativedelta(years=20))
    time_obs20 =ice_obs_post_mer_res20['DATE_SHIFT']
    area_inter_obs20 = ice_obs_post_mer_res20['sea_ice_area_interpolated']
if os.path.exists(file_obs_post_31):
    ice_obs_post_mer_res30["DATE_SHIFT"] = ice_obs_post_mer_res30["dtime"].apply(lambda x: x + relativedelta(years=30))
    time_obs30 =ice_obs_post_mer_res30['DATE_SHIFT']
    area_inter_obs30 = ice_obs_post_mer_res30['sea_ice_area_interpolated']





##--- Plot sea ice area -----------------------------------------------##
## Create a new figure and a set of subplots ##
fig1 = plt.figure(figsize=(10,6))
ax1 = fig1.add_subplot(111)

## Plot x, y ##
### Historical observation ###
hex_clr = ['#8e0152', '#c51b7d', '#de77ae', '#7fbc41', '#276419', '#734e30']
if start_month_mdl == 7:
    if os.path.exists(file_obs_post_2):
        ax1.plot(time_obs1, area_inter_obs1, color=hex_clr[0], linewidth = 1.5, label = 'Observation (' + str(start_year_obs -1) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_4):
        ax1.plot(time_obs3, area_inter_obs3, color=hex_clr[1], linewidth = 1.5, label = 'Observation (' + str(start_year_obs -3) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_6):
        ax1.plot(time_obs5, area_inter_obs5, color=hex_clr[2], linewidth = 1.5, label = 'Observation (' + str(start_year_obs -5) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_11):
        ax1.plot(time_obs10, area_inter_obs10, color=hex_clr[3], linewidth = 1.5, label = 'Observation (' + str(start_year_obs -10) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_21):
        ax1.plot(time_obs20, area_inter_obs20, color=hex_clr[4], linewidth = 1.5, label = 'Observation (' + str(start_year_obs -20) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_31):
        ax1.plot(time_obs30, area_inter_obs30, color=hex_clr[5], linewidth = 1.5, label = 'Observation (' + str(start_year_obs -30) + ')', alpha=0.7)
else:
    if os.path.exists(file_obs_post_2):
        ax1.plot(time_obs1, area_inter_obs1, color=hex_clr[0], linewidth = 1.5, label = 'Observation (' + str(start_year_obs-1) + '-' + str(end_year_mdl-1) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_4):
        ax1.plot(time_obs3, area_inter_obs3, color=hex_clr[1], linewidth = 1.5, label = 'Observation (' + str(start_year_obs-3) + '-' + str(end_year_mdl-3) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_6):
        ax1.plot(time_obs5, area_inter_obs5, color=hex_clr[2], linewidth = 1.5, label = 'Observation (' + str(start_year_obs-5) + '-' + str(end_year_mdl-5) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_11):
        ax1.plot(time_obs10, area_inter_obs10, color=hex_clr[3], linewidth = 1.5, label = 'Observation (' + str(start_year_obs-10) + '-' + str(end_year_mdl-10) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_21):
        ax1.plot(time_obs20, area_inter_obs20, color=hex_clr[4], linewidth = 1.5, label = 'Observation (' + str(start_year_obs-20) + '-' + str(end_year_mdl-20) + ')', alpha=0.7)
    if os.path.exists(file_obs_post_31):
        ax1.plot(time_obs30, area_inter_obs30, color=hex_clr[5], linewidth = 1.5, label = 'Observation (' + str(start_year_obs-30) + '-' + str(end_year_mdl-30) + ')', alpha=0.7)



### Latest observation ###
ax1.plot(time_obs, area_inter_obs, color='black', linewidth = 2.5, label= 'Observation (' + obs_start_date + '-' + obs_end_date + ')')


### Hindcast ###
ax1.fill_between(time_hind, percentile_10_calib, percentile_90_calib, facecolor='cyan', alpha=0.5, label = '10-90 percentile range')
ax1.plot(time_hind, ensemble_mean_calib, color='b', linewidth = 2.5, label = 'Ensemble mean (' + forecast_start_date + '-' + forecast_end_date + ')')


## Get legend, handles and labels ##
hans, labs = ax1.get_legend_handles_labels()



## Set xlim, ylim ##
span = [datetime.datetime(start_year_obs, start_month_obs, 1), datetime.datetime(end_year_mdl, end_month_mdl, calendar.monthrange(end_year_mdl, end_month_mdl)[1])]
ax1.set_xlim(span)
ax1.set_ylim([0, 16000000])


## Set X axis ##
#ax1.xaxis.set_major_locator(mdates.MonthLocator(interval=2)) # interval: 2 month
ax1.xaxis.set_major_locator(mdates.MonthLocator(interval=1)) # interval: 1 month
#datefmt = mdates.DateFormatter("%Y-%m") # format: yyyy-mm
datefmt = mdates.DateFormatter("%b") # format: Jan
ax1.xaxis.set_major_formatter(datefmt)


## Set Y axis ##
ax1.yaxis.set_major_formatter(ptick.ScalarFormatter(useMathText=True))
ax1.ticklabel_format(style="sci", axis="y", scilimits=(6,6)) # 1*10^6


## Set label size
ax1.tick_params(labelsize = 12, pad = 10)


## Add grid ##
ax1.minorticks_on()
ax1.grid(color = "gray", linestyle = "--")


## Add legend ##
### Set loc1 ###
if mm >= 8 or mm < 3:
    loc1 = 'lower right'
else:
    loc1 = 'upper right'
loc2 = 'lower left'

### Set legend number ###
ln_num = len(labs)-3
if ln_num > 0:
    l1 = ax1.legend(handles=hans[:ln_num], labels=labs[:ln_num], frameon=False, loc=loc1, fontsize = 'small')
    l2 = ax1.legend(handles=hans[ln_num:], labels=labs[ln_num:], frameon=False, loc=loc2, fontsize = 'small')
    ax1.add_artist(l1)
else:
    l2 = ax1.legend(handles=hans[ln_num:], labels=labs[ln_num:], frameon=False, loc=loc2, fontsize = 'small')


## Add labels, title ##
ax1.set_ylabel('Sea Ice Area [$\mathrm{km^2}$]', fontsize = 12)

fig_title = 'Daily Arctic sea ice area'
ax_title = '(forecast start: ' + forecast_start +')'

fig1.suptitle(fig_title, fontsize = 18)
ax1.set_title(ax_title, loc="right", fontsize = 14)

#plt.show()


## Save figures ##
figs_dir = out_dir + '/' + start_year + start_month + '/'
if not os.path.exists(figs_dir):
    os.makedirs(figs_dir)

### File name ###
'''
<bccr system number>_<variable name>_s<forecast start date (yyyymm)>_<aggregation>_<lead time>
'''
#fig_name1 = figs_dir + 'bccr_' + system_num + '_seaice_area_forecast_oisst.png' # <--- temp forlder
fig_name1 = figs_dir + 'bccr_' + system_num + '_seaice_area_s' + issue_date + '.png' # <--- archiving figure
fig1.savefig(fig_name1, dpi=600, format='png')

plt.show()