What if the ice block expedition 1959 happens in 2021?

def q_latent(air_pressure,air_temp_K,wind_2,es,ea): # downward positive

    # calculate latent heat # ASSUME: water is always available! satuation = 1
    # input:
    # - es, ea, air_pressure  in Pa
    # - air_temp_K in K
    # - wind_2 in m/s
    # output:
    # - q_lat in W/m-2
    
    rho_air = air_pressure/(287.058*air_temp_K)  # air density [kg /m^3], T in K
    lv = (2500.8-2.36*(air_temp_K-273.15))*1000; #latent heat of evaporation of water [J/kg], T in C
    A = 0.4*0.4 / np.log(2/1e-2) / np.log(2/1e-2) # k = 0.4
    
    q_lat = -0.622*rho_air*lv*A*wind_2*(es-ea)/air_pressure; # potential evaporation ASSUME: water is available

    return q_lat
    
def q_sensi(air_pressure,air_temp_K,wind_2,ta,ts): # downward positive
    # calculate sensible heat 
    # input:
    # - ta, ts in same unit, K or T
    # - air_temp_K in K
    # - wind_2 in m/s
    # output:
    # - q_sen in W/m-2
    cp_air = 1005  # heat capacity of air   J/(kg*K)

    rho_air = air_pressure/(287.058*air_temp_K)  # air density [kg /m^3], T in K
    lv = (2500.8-2.36*(air_temp_K-273.15))*1000; #latent heat of evaporation of water [J/kg], T in C
    A = 0.4*0.4 / np.log(2/1e-2) / np.log(2/1e-2) # k = 0.4
    
    q_sen = cp_air*rho_air*A*wind_2*(ta-ts)
    return q_sen    

## 2d heat equation

import copy
from numba import jit
import numpy as np
from scipy.ndimage import binary_dilation, binary_erosion

# 1 environment
box_t = np.loadtxt('1959_box_temp.csv')

# 2 settings
a_ice = 1.19e-6 # heat difussivity [m2/s]
a_air = 18.46e-6 # heat difussivity [m2/s]
a_gw = 1.785e-6 # heat difussivity [m2/s]

lf = 334000 #Latent heat of fusion L f = 334 000 J/kg -1 (phase change ice – water)
cp_ice = 2108 # J/(kgK)
rho_ice = 917 # Kg/m3

t_ice = -6 # initial temperature of ice
dx=0.01 # 1cm nodes
dt = 1 # seconds
hours = len(box_t) # hours here is 796
# 3 initial 
heatmap = np.zeros([ice_width,ice_width])+t_ice # inital temp
block = np.ones([ice_width,ice_width])   # initial musk   1:ice 0:air
tt = lf*rho_ice*0.01*0.01*0.01/(cp_ice*rho_ice*0.01*0.01*0.01)-t_ice  # energy musk: the energy required to melt in [K]
energy = block*tt

block[[0,-1],:]=0 # air boundary
block[:,[0,-1]]=0

energy[[0,-1],:]=0 # air boundary
energy[:,[0,-1]]=0

@jit(nopython=True, nogil=True)
def iteration(cs,inner_musk,energy,outer_musk,box_t):
    ns = cs.copy() # new state
    ns_energy = energy.copy() # new state
    ly,lx = ns.shape
    for i in range(0,lx):
        for j in range(0,ly):
            if outer_musk[j][i] == 1: #outer node
                #heat flux -> warm-up 
                #q = (box_t-cs[j][i])/0.25*0.04   # W/m2
                #dT = q/(917*0.01)/2108
                ns[j][i] = cs[j][i] +  (box_t-cs[j][i])/0.25*0.04/(917*0.01)/2108*dt
                ns_energy[j][i] = energy[j][i] - (box_t-cs[j][i])/0.25*0.04/(917*0.01)/2108*dt
                #ns_energy[j][i] = energy[j][i] - 30/3600*dt

                # ->add energy storage -> latent heat flux
                if ns[j][i]>0:
                    ns[j][i] = 0

            if inner_musk[j][i] == 1: #inner node
                ns[j][i] = cs[j][i] + a_ice*dt/dx**2 * (cs[j+1][i] + cs[j-1][i] +
                                                        cs[j][i+1]+cs[j][i-1]-
                                                        4*cs[j][i])
            else: #pass boundry conditon
                pass
    return ns,ns_energy
def find_outer_musk(a):
    a= a.astype(int)
    k = np.ones((3,3),dtype=int) # for 4-connected
    #k = np.zeros((3,3),dtype=int); k[1] = 1; k[:,1] = 1 # for 8-connected
    out = binary_dilation(a==0, k) & a
    inside = binary_erosion(a).astype(a.dtype)
    return out,inside

def solve_heat(heatmap,block,box_temp_list,energy):
    hm_list = []
    block_list = []
    energy_list = []
    inner_musk_list = []
    cs = heatmap.copy() # initial temp state
    energy_musk = energy.copy()
    block_musk = block.copy()

    for h in range(0,hours):     
        
        for t in range(1,3600): # iteration in seconds
            # HEAT diffusion
            outer_musk,inner_musk = find_outer_musk(block_musk)
            cs,energy_musk = iteration(cs,inner_musk,energy_musk,outer_musk,box_temp_list[h]) #update cs and energy_musk
            block_musk[energy_musk <0] = 0 # update blockmusk: TURN INTO AIR
        if h%1 ==0: # return every hour
            hm_list.append(cs)
            block_list.append(block_musk.copy())
            energy_list.append(energy_musk)
            inner_musk_list.append(inner_musk)
            print('Now: in hour',len(hm_list))
    return hm_list,block_list,energy_list,inner_musk_list

### Main

hm_list_2d,block_2d,energy_2d,inner_musk_2d = solve_heat(heatmap,block,box_t,energy)