# Read the daily sea ice concentration from hindcast/forecast simulations
# Calculate sea ice area in the arctic region (Northern Hemisphere, at least 15% sea ice concentration)
# Calculate ensemble mean , Standard deviation (SD), 10 percentile, 90 percentile
# Convert results into CSV file
#
# Script version 2.1
#
# Script works with:
#   Python version 3.10.13
#
#   Package version
#    numpy:  1.26.3
#    pandas:  2.2.0
#    netCDF4:  1.6.4
#
# Ver1.0.0: Created by Mariko Koseki, 30.08.2023
# Ver1.1.0: Updated by Mariko, 01.02.2024
# Ver1.1.1: Updated by Mariko, 13.03.2024
# Ver2.0: Updated by Mariko, 10.03 2025
# ver2.1: updated by mariko, 08.12.2025


import netCDF4
import numpy as np
import pandas as pd
#import datetime as dt
import statistics
import math
import collections
import os
import glob
import sys
import platform  # Ver1.1.0


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

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

#sys.exit()



##--- Input -----------------------------------------## #ver2.0.0
args = sys.argv
if len(args) == 2:
    date_ice = args[1]
    if len(date_ice) == 6:
        yyyy1 = int(date_ice[0:4])
        mm1 = int(date_ice[4:6])
        yyyy2 = yyyy1
        mm2 = mm1

    else:
        print('---How to use this script---')
        print('python read_netcdf_ice.py <yyyymm>')
        print('or')
        print('python read_netcdf_ice.py <yyyymm1> <yyyymm2>')
        print('<yyyymm1> = start month; <yyyymm2> = end month')
        print('')
        print('Example: python read_netcdf_ice.py 202312')
        print('Example: python read_netcdf_ice.py 199501 202312')
        print('')
        sys.exit()

elif len(args) == 3:
    date_ice_start = args[1]
    date_ice_end = args[2]
    if (len(date_ice_start) == 6) & (len(date_ice_end) == 6):
        yyyy1 = int(date_ice_start[0:4])
        mm1 = int(date_ice_start[4:6])
        yyyy2 = int(date_ice_end[0:4])
        mm2 = int(date_ice_end[4:6])

    else:
        print('---How to use this script---')
        print('python read_netcdf_ice.py <yyyymm>')
        print('or')
        print('python read_netcdf_ice.py <yyyymm1> <yyyymm2>')
        print('<yyyymm1> = start month; <yyyymm2> = end month')        
        print('')
        print('Example: python read_netcdf_ice.py 202312')
        print('Example: python read_netcdf_ice.py 199501 202312')
        print('')
        sys.exit()


else:
    print('---How to use this script---')
    print('python read_netcdf_ice.py <yyyymm>')
    print('or')
    print('python read_netcdf_ice.py <yyyymm1> <yyyymm2>')
    print('<yyyymm1> = start month; <yyyymm2> = end month')
    print('')
    print('Example: python read_netcdf_ice.py 202312')
    print('Example: python read_netcdf_ice.py 199501 202312')
    print('')
    sys.exit()


if mm2 == 12:
    start_month = mm1
    start_year = yyyy1
    end_month = 1
    end_year = yyyy2 + 1

else:
    start_month = mm1
    start_year = yyyy1
    end_month = mm2 + 1
    end_year = yyyy2


#print(start_month, start_year)
#print(end_month, end_year)
#sys.exit()



##--- Set path ------------------------------------------##
## Set path to input files ##
#path1 = '/nird/projects/NS9873K/norcpm/raw/norcpm-cf-system1/norcpm-cf-system1_hindcast/' # Ver1.1.1
path1 = '/nird/datalake/NS9873K/norcpm/raw/norcpm-cf-system1/norcpm-cf-system1_hindcast/' # Ver1.2.0
path10 = '/nird/datapeak/NS9034K/CMIP5/output1/NCC/NorESM1-M/piControl/fx/ocean/fx/r0i0p0/latest/areacello/areacello_fx_NorESM1-M_piControl_r0i0p0.nc' #ver2.1


## Set path to output file ## # Ver1.1.0
'''
Set 'out_hind' (path to csv file)
'''
seaice_dir = '/nird/datapeak/NS9873K/norcpm/validation/seaicepre/' # Ver1.1.1, ver2.1
out_hind = seaice_dir + 'hindcast/'

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

#sys.exit()



##--- Set year, month ------------------------------------## #ver2.0
'''
Set 'start_year' and 'end_year'
'''
start_date = '{:0>4d}-{:0>2d}'.format(start_year,start_month)
end_date = '{:0>4d}-{:0>2d}'.format(end_year,end_month)
ym_str_list = pd.date_range(start_date, end_date, freq='ME').strftime("%Y%m").tolist()

print(ym_str_list)

print('')
print('start month: ', ym_str_list[0])
print('end month: ', ym_str_list[-1])
print('')




for ym_path in ym_str_list:

    ## Search for directory ##
    path2 = 'norcpm-cf-system1_hindcast_' + ym_path + '15/norcpm-cf-system1_hindcast_' + ym_path + '15_mem01-60'  # Ver1.1.1

    path_dir_ice = path1 + path2
    print('')
    print('Search for dorectory:')
    print(path_dir_ice)
    print('')

    if not os.path.exists(path_dir_ice):
        print("no such a directory!!!")
        continue
    
    else:

        ice_files = glob.glob(path_dir_ice + "/ocn/hist/*hd*")
        print("file names: ", ice_files)


    ## Loop over files in the directory and read variable from each file ##
    all_sea_ice_mem01 = []; all_sea_ice_mem02 = []; all_sea_ice_mem03 = []; all_sea_ice_mem04 = []; all_sea_ice_mem05 = []
    all_sea_ice_mem06 = []; all_sea_ice_mem07 = []; all_sea_ice_mem08 = []; all_sea_ice_mem09 = []; all_sea_ice_mem10 = []
    all_sea_ice_mem11 = []; all_sea_ice_mem12 = []; all_sea_ice_mem13 = []; all_sea_ice_mem14 = []; all_sea_ice_mem15 = []
    all_sea_ice_mem16 = []; all_sea_ice_mem17 = []; all_sea_ice_mem18 = []; all_sea_ice_mem19 = []; all_sea_ice_mem20 = []
    all_sea_ice_mem21 = []; all_sea_ice_mem22 = []; all_sea_ice_mem23 = []; all_sea_ice_mem24 = []; all_sea_ice_mem25 = []
    all_sea_ice_mem26 = []; all_sea_ice_mem27 = []; all_sea_ice_mem28 = []; all_sea_ice_mem29 = []; all_sea_ice_mem30 = []
    all_sea_ice_mem31 = []; all_sea_ice_mem32 = []; all_sea_ice_mem33 = []; all_sea_ice_mem34 = []; all_sea_ice_mem35 = []
    all_sea_ice_mem36 = []; all_sea_ice_mem37 = []; all_sea_ice_mem38 = []; all_sea_ice_mem39 = []; all_sea_ice_mem40 = []
    all_sea_ice_mem41 = []; all_sea_ice_mem42 = []; all_sea_ice_mem43 = []; all_sea_ice_mem44 = []; all_sea_ice_mem45 = []
    all_sea_ice_mem46 = []; all_sea_ice_mem47 = []; all_sea_ice_mem48 = []; all_sea_ice_mem49 = []; all_sea_ice_mem50 = []
    all_sea_ice_mem51 = []; all_sea_ice_mem52 = []; all_sea_ice_mem53 = []; all_sea_ice_mem54 = []; all_sea_ice_mem55 = []
    all_sea_ice_mem56 = []; all_sea_ice_mem57 = []; all_sea_ice_mem58 = []; all_sea_ice_mem59 = []; all_sea_ice_mem60 = []
    
    all_ensemble_mean = []
    all_std_ens_ice = []
    all_percentile_10 = []
    all_percentile_90 = []
    all_dtime = []
    for nc_file in ice_files:
        print("file name: ", nc_file)
        #sys.exit()

        ## Read netCDF files ##
        nc_sea_ice = nc_file # sea ice area
        nc_grid_cell_area = path10 # grid cell area

        nc = netCDF4.Dataset(nc_sea_ice, 'r')
        gc = netCDF4.Dataset(nc_grid_cell_area, 'r')



        ## Check length of variables ##
        lon_len = len(nc.dimensions['x']) # Longitude: 320
        lat_len = len(nc.dimensions['y']) # Latitude:  384
        t_len = len(nc.dimensions['time']) # Time
        mem_len = len(nc.dimensions['ensemble']) # Ensemble member: 60



        #----- Read variables 1 -----

        ## Read time, lon, lat ##
        time = nc.variables['time'][:] # daily, units: days since 1800-01-01 00:00, NO LEAP
        timeunits = nc.variables['time'].getncattr('units')
        timecalendar = nc.variables['time'].getncattr('calendar')
        #dtime = [dt.date(1800,1,1) + dt.timedelta(days=t) for t in time] # does not work well!!
        dtime = netCDF4.num2date(time,timeunits,timecalendar)
        dtime = [ f'{i.year}-{i.month:02}-{i.day:02}' for i in dtime]

        #sys.exit()

        lon = gc.variables['lon'][:][:] # Longitude
        lat = gc.variables['lat'][:][:] # Latitude

        areacello = gc.variables['areacello'][:][:] # Ocean Grid-Cell Area (m2)
        area = areacello/1000000 # Ocean Grid-Cell Area (km2)


        ## Convert variables into Numpy ndarray 1 ##
        lon_ar = np.array(lon); lat_ar = np.array(lat)
        area_ar = np.array(area)


        ## Handle missing values and replace them to NaN ##
        '''
        Missing value:
            areacello:missing_value = 1.e+20f ;
            areacello:_FillValue = 1.e+20f
        '''
        area_ar[np.where(area_ar >= 1.0e+20)] = np.nan
        
        #print('Number of NaN in grid area: ', np.count_nonzero(np.isnan(area_ar)))




        #----- Read variables 2 -----

        list_sea_ice_mem01 = []; list_sea_ice_mem02 = []; list_sea_ice_mem03 = []; list_sea_ice_mem04 = []; list_sea_ice_mem05 = [];
        list_sea_ice_mem06 = []; list_sea_ice_mem07 = []; list_sea_ice_mem08 = []; list_sea_ice_mem09 = []; list_sea_ice_mem10 = [];
        list_sea_ice_mem11 = []; list_sea_ice_mem12 = []; list_sea_ice_mem13 = []; list_sea_ice_mem14 = []; list_sea_ice_mem15 = [];
        list_sea_ice_mem16 = []; list_sea_ice_mem17 = []; list_sea_ice_mem18 = []; list_sea_ice_mem19 = []; list_sea_ice_mem20 = [];
        list_sea_ice_mem21 = []; list_sea_ice_mem22 = []; list_sea_ice_mem23 = []; list_sea_ice_mem24 = []; list_sea_ice_mem25 = [];
        list_sea_ice_mem26 = []; list_sea_ice_mem27 = []; list_sea_ice_mem28 = []; list_sea_ice_mem29 = []; list_sea_ice_mem30 = [];
        list_sea_ice_mem31 = []; list_sea_ice_mem32 = []; list_sea_ice_mem33 = []; list_sea_ice_mem34 = []; list_sea_ice_mem35 = [];
        list_sea_ice_mem36 = []; list_sea_ice_mem37 = []; list_sea_ice_mem38 = []; list_sea_ice_mem39 = []; list_sea_ice_mem40 = [];
        list_sea_ice_mem41 = []; list_sea_ice_mem42 = []; list_sea_ice_mem43 = []; list_sea_ice_mem44 = []; list_sea_ice_mem45 = [];
        list_sea_ice_mem46 = []; list_sea_ice_mem47 = []; list_sea_ice_mem48 = []; list_sea_ice_mem49 = []; list_sea_ice_mem50 = [];
        list_sea_ice_mem51 = []; list_sea_ice_mem52 = []; list_sea_ice_mem53 = []; list_sea_ice_mem54 = []; list_sea_ice_mem55 = [];
        list_sea_ice_mem56 = []; list_sea_ice_mem57 = []; list_sea_ice_mem58 = []; list_sea_ice_mem59 = []; list_sea_ice_mem60 = [];

        list_ensemble_mean = []
        list_std_ens_ice = []
        list_percentile_10 = []
        list_percentile_90 = []
        for t in range(0, t_len):
            sea_ice_60mem = []
            for mem in range(0, mem_len):
                #print(mem)
                var = np.squeeze(nc.variables['fice'][mem][t][:][:]) # Ice concentration , units: %, fice(ensemble, time, y (lat), x (lon))


                ## Convert variables into Numpy ndarray 2 ##
                var_ar = np.array(var)




                #----- Calculate sea ice extent -----
                '''
                The cumulative area of all polar grid cells of the Northern Hemisphere that have at least 15% sea ice concentration
                '''

                '''
                Missing value:
                    fice:_FillValue = 9.96921e+36
                '''

                ## Handle missing values and replace them to NaN ##
                var_ar[np.where(var_ar >= 9.96921e+36)] = np.nan

                #print('Number of NaN in sea ice: ', np.count_nonzero(np.isnan(var_ar)))



                ## Check index of the Northern Hemispher ##
                lat_nor_ind = np.where(lat_ar > 0)[0][0]
                lat_sou_ind = max(np.where(lat_ar <= 0)[:][0])

                #print(lat_nor_ind)
                #print(lat_sou_ind)



                ## Extract sea ice and grid area only in the Northern Hemispher ##
                var_nor = var_ar[lat_nor_ind:lat_len, :]
                var_sou = var_ar[0:lat_sou_ind+1, :]
                area_nor = area_ar[lat_nor_ind:lat_len, :]



                ## Extract grid cells of the Northern Hemisphere that have at least 15% sea ice concentration ##
                
                #print('Number of grid cells (>= 15%): ', np.count_nonzero(var_nor >= 15))
                #print('Number of grid cells (< 15%): ', np.count_nonzero(var_nor < 15))

                var_nor[np.where(var_nor < 15)] = np.nan
                
                #print('Number of NaN in sea ice (Northern Hemisphere): ', np.count_nonzero(np.isnan(var_nor)))



                ## Calculate sea ice extent ##
                sea_ice_extent = var_nor * area_nor/100 # unit of sea ice: %

                #print('Number of NaN in sea ice extent: ', np.count_nonzero(np.isnan(sea_ice_extent)))


                ## Calculate sum of sea ice extent ##
                sum_sea_ice_extent = np.nansum(sea_ice_extent)


                #print(sum_sea_ice_extent)


                ## Add results into list ##
                sea_ice_60mem.append(sum_sea_ice_extent)



            # ----- Calculate ensemble mean, SD -----

            sea_ice_mem01 = sea_ice_60mem[0]; sea_ice_mem02 = sea_ice_60mem[1]; sea_ice_mem03 = sea_ice_60mem[2]
            sea_ice_mem04 = sea_ice_60mem[3]; sea_ice_mem05 = sea_ice_60mem[4]; sea_ice_mem06 = sea_ice_60mem[5]
            sea_ice_mem07 = sea_ice_60mem[6]; sea_ice_mem08 = sea_ice_60mem[7] ;sea_ice_mem09 = sea_ice_60mem[8]
            sea_ice_mem10 = sea_ice_60mem[9]; sea_ice_mem11 = sea_ice_60mem[10]; sea_ice_mem12 = sea_ice_60mem[11]
            sea_ice_mem13 = sea_ice_60mem[12]; sea_ice_mem14 = sea_ice_60mem[13]; sea_ice_mem15 = sea_ice_60mem[14]
            sea_ice_mem16 = sea_ice_60mem[15]; sea_ice_mem17 = sea_ice_60mem[16]; sea_ice_mem18 = sea_ice_60mem[17]
            sea_ice_mem19 = sea_ice_60mem[18]; sea_ice_mem20 = sea_ice_60mem[19]; sea_ice_mem21 = sea_ice_60mem[20]
            sea_ice_mem22 = sea_ice_60mem[21]; sea_ice_mem23 = sea_ice_60mem[22]; sea_ice_mem24 = sea_ice_60mem[23]
            sea_ice_mem25 = sea_ice_60mem[24]; sea_ice_mem26 = sea_ice_60mem[25]; sea_ice_mem27 = sea_ice_60mem[26]
            sea_ice_mem28 = sea_ice_60mem[27]; sea_ice_mem29 = sea_ice_60mem[28]; sea_ice_mem30 = sea_ice_60mem[29]
            sea_ice_mem31 = sea_ice_60mem[30]; sea_ice_mem32 = sea_ice_60mem[31]; sea_ice_mem33 = sea_ice_60mem[32]
            sea_ice_mem34 = sea_ice_60mem[33]; sea_ice_mem35 = sea_ice_60mem[34]; sea_ice_mem36 = sea_ice_60mem[35]
            sea_ice_mem37 = sea_ice_60mem[36]; sea_ice_mem38 = sea_ice_60mem[37]; sea_ice_mem39 = sea_ice_60mem[38]
            sea_ice_mem40 = sea_ice_60mem[39]; sea_ice_mem41 = sea_ice_60mem[40]; sea_ice_mem42 = sea_ice_60mem[41]
            sea_ice_mem43 = sea_ice_60mem[42]; sea_ice_mem44 = sea_ice_60mem[43]; sea_ice_mem45 = sea_ice_60mem[44]
            sea_ice_mem46 = sea_ice_60mem[45]; sea_ice_mem47 = sea_ice_60mem[46]; sea_ice_mem48 = sea_ice_60mem[47]
            sea_ice_mem49 = sea_ice_60mem[48]; sea_ice_mem50 = sea_ice_60mem[49]; sea_ice_mem51 = sea_ice_60mem[50]
            sea_ice_mem52 = sea_ice_60mem[51]; sea_ice_mem53 = sea_ice_60mem[52]; sea_ice_mem54 = sea_ice_60mem[53]
            sea_ice_mem55 = sea_ice_60mem[54]; sea_ice_mem56 = sea_ice_60mem[55]; sea_ice_mem57 = sea_ice_60mem[56]
            sea_ice_mem58 = sea_ice_60mem[57]; sea_ice_mem59 = sea_ice_60mem[58]; sea_ice_mem60 = sea_ice_60mem[59]


            #sys.exit()
            ens_mean_ice2 = statistics.mean(sea_ice_60mem)
            std_ens_ice = statistics.stdev(sea_ice_60mem)



            # ----- Calculate percentiles -----
        
            def calculate_percentile(arry, percentile):
                size = len(arry)
                return sorted(arry)[int(math.ceil((size * percentile) / 100)) - 1]

            percentile_10 = calculate_percentile(sea_ice_60mem, 10)
            percentile_90 = calculate_percentile(sea_ice_60mem, 90)

            #print(percentile_10)
            #print(percentile_90)



            # ----- Add results into list -----
            
            list_sea_ice_mem01.append(sea_ice_mem01); list_sea_ice_mem02.append(sea_ice_mem02); list_sea_ice_mem03.append(sea_ice_mem03);
            list_sea_ice_mem04.append(sea_ice_mem04); list_sea_ice_mem05.append(sea_ice_mem05); list_sea_ice_mem06.append(sea_ice_mem06);
            list_sea_ice_mem07.append(sea_ice_mem07); list_sea_ice_mem08.append(sea_ice_mem08); list_sea_ice_mem09.append(sea_ice_mem09);
            list_sea_ice_mem10.append(sea_ice_mem10); list_sea_ice_mem11.append(sea_ice_mem11); list_sea_ice_mem12.append(sea_ice_mem12);
            list_sea_ice_mem13.append(sea_ice_mem13); list_sea_ice_mem14.append(sea_ice_mem14); list_sea_ice_mem15.append(sea_ice_mem15);
            list_sea_ice_mem16.append(sea_ice_mem16); list_sea_ice_mem17.append(sea_ice_mem17); list_sea_ice_mem18.append(sea_ice_mem18);
            list_sea_ice_mem19.append(sea_ice_mem19); list_sea_ice_mem20.append(sea_ice_mem20); list_sea_ice_mem21.append(sea_ice_mem21);
            list_sea_ice_mem22.append(sea_ice_mem22); list_sea_ice_mem23.append(sea_ice_mem23); list_sea_ice_mem24.append(sea_ice_mem24);
            list_sea_ice_mem25.append(sea_ice_mem25); list_sea_ice_mem26.append(sea_ice_mem26); list_sea_ice_mem27.append(sea_ice_mem27);
            list_sea_ice_mem28.append(sea_ice_mem28); list_sea_ice_mem29.append(sea_ice_mem29); list_sea_ice_mem30.append(sea_ice_mem30);
            list_sea_ice_mem31.append(sea_ice_mem31); list_sea_ice_mem32.append(sea_ice_mem32); list_sea_ice_mem33.append(sea_ice_mem33);
            list_sea_ice_mem34.append(sea_ice_mem34); list_sea_ice_mem35.append(sea_ice_mem35); list_sea_ice_mem36.append(sea_ice_mem36);
            list_sea_ice_mem37.append(sea_ice_mem37); list_sea_ice_mem38.append(sea_ice_mem38); list_sea_ice_mem39.append(sea_ice_mem39);
            list_sea_ice_mem40.append(sea_ice_mem40); list_sea_ice_mem41.append(sea_ice_mem41); list_sea_ice_mem42.append(sea_ice_mem42);
            list_sea_ice_mem43.append(sea_ice_mem43); list_sea_ice_mem44.append(sea_ice_mem44); list_sea_ice_mem45.append(sea_ice_mem45);
            list_sea_ice_mem46.append(sea_ice_mem46); list_sea_ice_mem47.append(sea_ice_mem47); list_sea_ice_mem48.append(sea_ice_mem48);
            list_sea_ice_mem49.append(sea_ice_mem49); list_sea_ice_mem50.append(sea_ice_mem50); list_sea_ice_mem51.append(sea_ice_mem51);
            list_sea_ice_mem52.append(sea_ice_mem52); list_sea_ice_mem53.append(sea_ice_mem53); list_sea_ice_mem54.append(sea_ice_mem54);
            list_sea_ice_mem55.append(sea_ice_mem55); list_sea_ice_mem56.append(sea_ice_mem56); list_sea_ice_mem57.append(sea_ice_mem57);
            list_sea_ice_mem58.append(sea_ice_mem58); list_sea_ice_mem59.append(sea_ice_mem59); list_sea_ice_mem60.append(sea_ice_mem60);

            list_ensemble_mean.append(ens_mean_ice2)
            list_std_ens_ice.append(std_ens_ice)
            list_percentile_10.append(percentile_10)
            list_percentile_90.append(percentile_90)

        # Close netCDF files
        nc.close()
        gc.close()

 
        #print(list_ensemble_mean)
        all_sea_ice_mem01.append(list_sea_ice_mem01); all_sea_ice_mem02.append(list_sea_ice_mem02); all_sea_ice_mem03.append(list_sea_ice_mem03);
        all_sea_ice_mem04.append(list_sea_ice_mem04); all_sea_ice_mem05.append(list_sea_ice_mem05); all_sea_ice_mem06.append(list_sea_ice_mem06);
        all_sea_ice_mem07.append(list_sea_ice_mem07); all_sea_ice_mem08.append(list_sea_ice_mem08); all_sea_ice_mem09.append(list_sea_ice_mem09);
        all_sea_ice_mem10.append(list_sea_ice_mem10); all_sea_ice_mem11.append(list_sea_ice_mem11); all_sea_ice_mem12.append(list_sea_ice_mem12);
        all_sea_ice_mem13.append(list_sea_ice_mem13); all_sea_ice_mem14.append(list_sea_ice_mem14); all_sea_ice_mem15.append(list_sea_ice_mem15);
        all_sea_ice_mem16.append(list_sea_ice_mem16); all_sea_ice_mem17.append(list_sea_ice_mem17); all_sea_ice_mem18.append(list_sea_ice_mem18);
        all_sea_ice_mem19.append(list_sea_ice_mem19); all_sea_ice_mem20.append(list_sea_ice_mem20); all_sea_ice_mem21.append(list_sea_ice_mem21);
        all_sea_ice_mem22.append(list_sea_ice_mem22); all_sea_ice_mem23.append(list_sea_ice_mem23); all_sea_ice_mem24.append(list_sea_ice_mem24);
        all_sea_ice_mem25.append(list_sea_ice_mem25); all_sea_ice_mem26.append(list_sea_ice_mem26); all_sea_ice_mem27.append(list_sea_ice_mem27);
        all_sea_ice_mem28.append(list_sea_ice_mem28); all_sea_ice_mem29.append(list_sea_ice_mem29); all_sea_ice_mem30.append(list_sea_ice_mem30);
        all_sea_ice_mem31.append(list_sea_ice_mem31); all_sea_ice_mem32.append(list_sea_ice_mem32); all_sea_ice_mem33.append(list_sea_ice_mem33);
        all_sea_ice_mem34.append(list_sea_ice_mem34); all_sea_ice_mem35.append(list_sea_ice_mem35); all_sea_ice_mem36.append(list_sea_ice_mem36);
        all_sea_ice_mem37.append(list_sea_ice_mem37); all_sea_ice_mem38.append(list_sea_ice_mem38); all_sea_ice_mem39.append(list_sea_ice_mem39);
        all_sea_ice_mem40.append(list_sea_ice_mem40); all_sea_ice_mem41.append(list_sea_ice_mem41); all_sea_ice_mem42.append(list_sea_ice_mem42);
        all_sea_ice_mem43.append(list_sea_ice_mem43); all_sea_ice_mem44.append(list_sea_ice_mem44); all_sea_ice_mem45.append(list_sea_ice_mem45);
        all_sea_ice_mem46.append(list_sea_ice_mem46); all_sea_ice_mem47.append(list_sea_ice_mem47); all_sea_ice_mem48.append(list_sea_ice_mem48);
        all_sea_ice_mem49.append(list_sea_ice_mem49); all_sea_ice_mem50.append(list_sea_ice_mem50); all_sea_ice_mem51.append(list_sea_ice_mem51);
        all_sea_ice_mem52.append(list_sea_ice_mem52); all_sea_ice_mem53.append(list_sea_ice_mem53); all_sea_ice_mem54.append(list_sea_ice_mem54);
        all_sea_ice_mem55.append(list_sea_ice_mem55); all_sea_ice_mem56.append(list_sea_ice_mem56); all_sea_ice_mem57.append(list_sea_ice_mem57);
        all_sea_ice_mem58.append(list_sea_ice_mem58); all_sea_ice_mem59.append(list_sea_ice_mem59); all_sea_ice_mem60.append(list_sea_ice_mem60);

        all_ensemble_mean.append(list_ensemble_mean)
        all_std_ens_ice.append(list_std_ens_ice)
        all_percentile_10.append(list_percentile_10)
        all_percentile_90.append(list_percentile_90)
        all_dtime.append(dtime)




    ## Multidimensional list -> 1D list ##
    def flatten(lst):
        for item in lst:
            if isinstance(item, collections.abc.Iterable) and not isinstance(item, (str, bytes)):
                yield from flatten(item)
            else:
                yield item

    fla01 = list(flatten(all_sea_ice_mem01)); fla02 = list(flatten(all_sea_ice_mem02)); fla03 = list(flatten(all_sea_ice_mem03));
    fla04 = list(flatten(all_sea_ice_mem04)); fla05 = list(flatten(all_sea_ice_mem05)); fla06 = list(flatten(all_sea_ice_mem06));
    fla07 = list(flatten(all_sea_ice_mem07)); fla08 = list(flatten(all_sea_ice_mem08)); fla09 = list(flatten(all_sea_ice_mem09));
    fla10 = list(flatten(all_sea_ice_mem10)); fla11 = list(flatten(all_sea_ice_mem11)); fla12 = list(flatten(all_sea_ice_mem12));
    fla13 = list(flatten(all_sea_ice_mem13)); fla14 = list(flatten(all_sea_ice_mem14)); fla15 = list(flatten(all_sea_ice_mem15));
    fla16 = list(flatten(all_sea_ice_mem16)); fla17 = list(flatten(all_sea_ice_mem17)); fla18 = list(flatten(all_sea_ice_mem18));
    fla19 = list(flatten(all_sea_ice_mem19)); fla20 = list(flatten(all_sea_ice_mem20)); fla21 = list(flatten(all_sea_ice_mem21));
    fla22 = list(flatten(all_sea_ice_mem22)); fla23 = list(flatten(all_sea_ice_mem23)); fla24 = list(flatten(all_sea_ice_mem24));
    fla25 = list(flatten(all_sea_ice_mem25)); fla26 = list(flatten(all_sea_ice_mem26)); fla27 = list(flatten(all_sea_ice_mem27));
    fla28 = list(flatten(all_sea_ice_mem28)); fla29 = list(flatten(all_sea_ice_mem29)); fla30 = list(flatten(all_sea_ice_mem30));
    fla31 = list(flatten(all_sea_ice_mem31)); fla32 = list(flatten(all_sea_ice_mem32)); fla33 = list(flatten(all_sea_ice_mem33));
    fla34 = list(flatten(all_sea_ice_mem34)); fla35 = list(flatten(all_sea_ice_mem35)); fla36 = list(flatten(all_sea_ice_mem36));
    fla37 = list(flatten(all_sea_ice_mem37)); fla38 = list(flatten(all_sea_ice_mem38)); fla39 = list(flatten(all_sea_ice_mem39));
    fla40 = list(flatten(all_sea_ice_mem40)); fla41 = list(flatten(all_sea_ice_mem41)); fla42 = list(flatten(all_sea_ice_mem42));
    fla43 = list(flatten(all_sea_ice_mem43)); fla44 = list(flatten(all_sea_ice_mem44)); fla45 = list(flatten(all_sea_ice_mem45));
    fla46 = list(flatten(all_sea_ice_mem46)); fla47 = list(flatten(all_sea_ice_mem47)); fla48 = list(flatten(all_sea_ice_mem48));
    fla49 = list(flatten(all_sea_ice_mem49)); fla50 = list(flatten(all_sea_ice_mem50)); fla51 = list(flatten(all_sea_ice_mem51));
    fla52 = list(flatten(all_sea_ice_mem52)); fla53 = list(flatten(all_sea_ice_mem53)); fla54 = list(flatten(all_sea_ice_mem54));
    fla55 = list(flatten(all_sea_ice_mem55)); fla56 = list(flatten(all_sea_ice_mem56)); fla57 = list(flatten(all_sea_ice_mem57));
    fla58 = list(flatten(all_sea_ice_mem58)); fla59 = list(flatten(all_sea_ice_mem59)); fla60 = list(flatten(all_sea_ice_mem60));


    fla61 = list(flatten(all_ensemble_mean))
    fla62 = list(flatten(all_dtime))
    fla63 = list(flatten(all_percentile_10))
    fla64 = list(flatten(all_percentile_90))
    fla65 = list(flatten(all_std_ens_ice))


    ## List -> Numpy array ##
    all_sea_ice_mem01_ar = np.array(fla01); all_sea_ice_mem02_ar = np.array(fla02); all_sea_ice_mem03_ar = np.array(fla03);
    all_sea_ice_mem04_ar = np.array(fla04); all_sea_ice_mem05_ar = np.array(fla05); all_sea_ice_mem06_ar = np.array(fla06);
    all_sea_ice_mem07_ar = np.array(fla07); all_sea_ice_mem08_ar = np.array(fla08); all_sea_ice_mem09_ar = np.array(fla09);
    all_sea_ice_mem10_ar = np.array(fla10); all_sea_ice_mem11_ar = np.array(fla11); all_sea_ice_mem12_ar = np.array(fla12);
    all_sea_ice_mem13_ar = np.array(fla13); all_sea_ice_mem14_ar = np.array(fla14); all_sea_ice_mem15_ar = np.array(fla15);
    all_sea_ice_mem16_ar = np.array(fla16); all_sea_ice_mem17_ar = np.array(fla17); all_sea_ice_mem18_ar = np.array(fla18);
    all_sea_ice_mem19_ar = np.array(fla19); all_sea_ice_mem20_ar = np.array(fla20); all_sea_ice_mem21_ar = np.array(fla21);
    all_sea_ice_mem22_ar = np.array(fla22); all_sea_ice_mem23_ar = np.array(fla23); all_sea_ice_mem24_ar = np.array(fla24);
    all_sea_ice_mem25_ar = np.array(fla25); all_sea_ice_mem26_ar = np.array(fla26); all_sea_ice_mem27_ar = np.array(fla27);
    all_sea_ice_mem28_ar = np.array(fla28); all_sea_ice_mem29_ar = np.array(fla29); all_sea_ice_mem30_ar = np.array(fla30);
    all_sea_ice_mem31_ar = np.array(fla31); all_sea_ice_mem32_ar = np.array(fla32); all_sea_ice_mem33_ar = np.array(fla33);
    all_sea_ice_mem34_ar = np.array(fla34); all_sea_ice_mem35_ar = np.array(fla35); all_sea_ice_mem36_ar = np.array(fla36);
    all_sea_ice_mem37_ar = np.array(fla37); all_sea_ice_mem38_ar = np.array(fla38); all_sea_ice_mem39_ar = np.array(fla39);
    all_sea_ice_mem40_ar = np.array(fla40); all_sea_ice_mem41_ar = np.array(fla41); all_sea_ice_mem42_ar = np.array(fla42);
    all_sea_ice_mem43_ar = np.array(fla43); all_sea_ice_mem44_ar = np.array(fla44); all_sea_ice_mem45_ar = np.array(fla45);
    all_sea_ice_mem46_ar = np.array(fla46); all_sea_ice_mem47_ar = np.array(fla47); all_sea_ice_mem48_ar = np.array(fla48);
    all_sea_ice_mem49_ar = np.array(fla49); all_sea_ice_mem50_ar = np.array(fla50); all_sea_ice_mem51_ar = np.array(fla51);
    all_sea_ice_mem52_ar = np.array(fla52); all_sea_ice_mem53_ar = np.array(fla53); all_sea_ice_mem54_ar = np.array(fla54);
    all_sea_ice_mem55_ar = np.array(fla55); all_sea_ice_mem56_ar = np.array(fla56); all_sea_ice_mem57_ar = np.array(fla57);
    all_sea_ice_mem58_ar = np.array(fla58); all_sea_ice_mem59_ar = np.array(fla59); all_sea_ice_mem60_ar = np.array(fla60);



    all_ensemble_mean_ar = np.array(fla61)
    all_dtime_ar = np.array(fla62)
    all_percentile_10_ar = np.array(fla63)
    all_percentile_90_ar = np.array(fla64)
    all_std_ens_ice_ar = np.array(fla65)


    ## Stack 1-D arrays as columns into a 2-D array ##
    np_arr = np.column_stack([all_dtime_ar, all_sea_ice_mem01_ar, all_sea_ice_mem02_ar, all_sea_ice_mem03_ar, all_sea_ice_mem04_ar, \
            all_sea_ice_mem05_ar, all_sea_ice_mem06_ar, all_sea_ice_mem07_ar, all_sea_ice_mem08_ar, all_sea_ice_mem09_ar, all_sea_ice_mem10_ar, \
            all_sea_ice_mem11_ar, all_sea_ice_mem12_ar, all_sea_ice_mem13_ar, all_sea_ice_mem14_ar, all_sea_ice_mem15_ar, \
            all_sea_ice_mem16_ar, all_sea_ice_mem17_ar, all_sea_ice_mem18_ar, all_sea_ice_mem19_ar, all_sea_ice_mem20_ar, \
            all_sea_ice_mem21_ar, all_sea_ice_mem22_ar, all_sea_ice_mem23_ar, all_sea_ice_mem24_ar, all_sea_ice_mem25_ar, \
            all_sea_ice_mem26_ar, all_sea_ice_mem27_ar, all_sea_ice_mem28_ar, all_sea_ice_mem29_ar, all_sea_ice_mem30_ar, \
            all_sea_ice_mem31_ar, all_sea_ice_mem32_ar, all_sea_ice_mem33_ar, all_sea_ice_mem34_ar, all_sea_ice_mem35_ar, \
            all_sea_ice_mem36_ar, all_sea_ice_mem37_ar, all_sea_ice_mem38_ar, all_sea_ice_mem39_ar, all_sea_ice_mem40_ar, \
            all_sea_ice_mem41_ar, all_sea_ice_mem42_ar, all_sea_ice_mem43_ar, all_sea_ice_mem44_ar, all_sea_ice_mem45_ar, \
            all_sea_ice_mem46_ar, all_sea_ice_mem47_ar, all_sea_ice_mem48_ar, all_sea_ice_mem49_ar, all_sea_ice_mem50_ar, \
            all_sea_ice_mem51_ar, all_sea_ice_mem52_ar, all_sea_ice_mem53_ar, all_sea_ice_mem54_ar, all_sea_ice_mem55_ar, \
            all_sea_ice_mem56_ar, all_sea_ice_mem57_ar, all_sea_ice_mem58_ar, all_sea_ice_mem59_ar, all_sea_ice_mem60_ar, \
            all_ensemble_mean_ar, all_std_ens_ice_ar, all_percentile_10_ar, all_percentile_90_ar])

    #print(np_arr)



    # -----Check if directory for output exists and create it if not -----

    outdir = out_hind + 'norcpm-cf-system1_hindcast1_' + ym_path + '15_mem01-60/'
    if not os.path.exists(outdir):
        os.makedirs(outdir)



    # ----- Save results as CSV -----

    ## Convert all results into Pandas DataFrame ##
    df = pd.DataFrame(np_arr, columns=["dtime", "ens_mem01", "ens_mem02", "ens_mem03", "ens_mem04", "ens_mem05", \
            "ens_mem06", "ens_mem07", "ens_mem08", "ens_mem09", "ens_mem10", \
            "ens_mem11", "ens_mem12", "ens_mem13", "ens_mem14", "ens_mem15", \
            "ens_mem16", "ens_mem17", "ens_mem18", "ens_mem19", "ens_mem20", \
            "ens_mem21", "ens_mem22", "ens_mem23", "ens_mem24", "ens_mem25", \
            "ens_mem26", "ens_mem27", "ens_mem28", "ens_mem29", "ens_mem30", \
            "ens_mem31", "ens_mem32", "ens_mem33", "ens_mem34", "ens_mem35", \
            "ens_mem36", "ens_mem37", "ens_mem38", "ens_mem39", "ens_mem40", \
            "ens_mem41", "ens_mem42", "ens_mem43", "ens_mem44", "ens_mem45", \
            "ens_mem46", "ens_mem47", "ens_mem48", "ens_mem49", "ens_mem50", \
            "ens_mem51", "ens_mem52", "ens_mem53", "ens_mem54", "ens_mem55", \
            "ens_mem56", "ens_mem57", "ens_mem58", "ens_mem59", "ens_mem60", \
            "ensemble_mean", "std_ens_ice", "percentile_10", "percentile_90"])


    ## Sort dataset by dtime ##
    df = df.sort_values('dtime');


    ## Save DataFrame as CSV ##
    df.to_csv(outdir + 'seaice_hind_' + ym_path + '.csv', index=False)


print('')
print('Completed!')