diff --git a/CFM_main/CFMoutput_example/csv/CFMresults.hdf5 b/CFM_main/CFMoutput_example/csv/CFMresults.hdf5 index fe56d56..e58b6fd 100644 Binary files a/CFM_main/CFMoutput_example/csv/CFMresults.hdf5 and b/CFM_main/CFMoutput_example/csv/CFMresults.hdf5 differ diff --git a/CFM_main/CFMoutput_example/csv/CFMspin.hdf5 b/CFM_main/CFMoutput_example/csv/CFMspin.hdf5 index 6161a22..b107e1c 100644 Binary files a/CFM_main/CFMoutput_example/csv/CFMspin.hdf5 and b/CFM_main/CFMoutput_example/csv/CFMspin.hdf5 differ diff --git a/CFM_main/CFMoutput_example/csv/example_csv.json b/CFM_main/CFMoutput_example/csv/example_csv.json deleted file mode 100644 index c985ea6..0000000 --- a/CFM_main/CFMoutput_example/csv/example_csv.json +++ /dev/null @@ -1,112 +0,0 @@ -{ -"InputFileFolder": "CFMinput_example", -"InputFileNameTemp": "example_TSKIN.csv", -"InputFileNamebdot": "example_BDOT.csv", -"InputFileNameIso": "example_ISOTOPE.csv", -"InputFileNamerho": "example_RHOS.csv", -"InputFileNamemelt": "example_SMELT.csv", -"InputFileNameStrain": "example_STRAIN.csv", -"InputFileNameSublim": "example_SUBLIM.csv", -"InputFileNameRain": "example_RAIN.csv", -"resultsFolder": "CFMoutput_example/csv", -"initfirnFile": "example_firndata.csv", -"initprofile": false, -"input_type": "csv", -"input_type_options": ["csv","dataframe"], -"DFresample": "1M", -"DFfile": "CFM_example_input_df.pkl", -"physRho": "KuipersMunneke2015", -"physRho_options":["HLdynamic","HLSigfus","Li2004","Li2011","Helsen2008","Arthern2010S","Arthern2010T","Li2015","Goujon2003","Barnola1991","Morris2014","KuipersMunneke2015","Crocus","Ligtenberg2011"], -"MELT": true, -"RAIN": true, -"ReehCorrectedT": false, -"FirnAir": false, -"AirConfigName": "AirConfig.json", -"TWriteInt": 1, -"TWriteStart": 1980.0, -"int_type": "nearest", -"int_type_options": ["nearest","linear"], -"SeasonalTcycle": false, -"SeasonalThemi": "north", -"coreless": true, -"TAmp":10.0, -"physGrain": true, -"calcGrainSize": false, -"GrGrowPhysics": "Arthern", -"GrGrowPhysics_options": ["Arthern", "Katsushima"], -"heatDiff": true, -"conductivity": "Calonne2019", -"conductivity_options": ["Schwander","Yen_fixed","Yen_var","Anderson","Yen_b","Sturm","VanDusen","Schwerdtfeger","Riche","Jiawen","mix","Calonne2011","Calonne2019"], -"variable_srho": false, -"srho_type": "userinput", -"srho_type_options": ["userinput","param","noise"], -"rhos0": 350.0, -"r2s0": 1.0e-8, -"AutoSpinUpTime": false, -"yearSpin": 1, -"H": 3000, -"HbaseSpin": 2880.0, -"stpsPerYear": 12.0, -"D_surf": 1.0, -"bdot_type": "mean", -"bdot_type_options": ["instant","mean","stress"], -"grid_outputs": true, -"grid_output_res": 0.25, -"isoDiff": true, -"iso": ["d18O","dD"], -"isoOptions":["d18O","dD","NoDiffusion"], -"spacewriteint": 1, -"horizontal_divergence": false, -"strain_softening": false, -"tuning_bias_correction": false, -"residual_strain": 0.0002, -"outputs": ["density", "depth", "temperature", "age", "DIP","LWC","meltoutputs","climate"], -"outputs_options": ["density", "depth", "temperature", "age", "Dcon", "bdot_mean", "climate", "compaction", "grainsize", "temp_Hx", "isotopes", "BCO", "DIPc", "DIP", "LWC","gasses", "PLWC_mem", "viscosity", "runoff", "refrozen"], -"resultsFileName": "CFMresults.hdf5", -"spinFileName": "CFMspin.hdf5", -"doublegrid": true, -"nodestocombine": 3, -"multnodestocombine": 12, -"grid1bottom": 5.0, -"grid2bottom": 15.0, -"spinup_climate_type": "mean", -"spinup_climate_type_options": ["mean", "initial"], -"manual_climate": false, -"deepT": 255.88, -"bdot_long": 0.49073, -"manual_iceout": false, -"iceout": 0.23, -"QMorris": 110.0e3, -"timesetup": "exact", -"timesetup_options": ["exact","interp","retmip"], -"liquid": "bucket", -"liquid_options": ["bucket","darcy","resingledomain","prefsnowpack"], -"merging": true, -"merge_min": 1e-9, -"LWCcorrect": false, -"manualT": false, -"no_densification": false, -"rad_pen": false, -"site_pressure": 1013.25, -"output_bits": "float32", -"spinUpdate": true, -"spinUpdateDate": 1980.0, -"DIPhorizon": 100.0, -"NewSpin": false, -"ColeouLesaffre": true, -"IrrVal": 0.02, -"RhoImp": 830, -"DownToIce": false, -"ThickImp": 0.1, -"Ponding": false, -"DirectRunoff": 0.0, -"RunoffZuoOerlemans": false, -"Slope": 0.1, -"SUBLIM": true -} - - - - - - diff --git a/CFM_main/RCMpkl_to_spin.py b/CFM_main/RCMpkl_to_spin.py index c4b0334..e8276cd 100644 --- a/CFM_main/RCMpkl_to_spin.py +++ b/CFM_main/RCMpkl_to_spin.py @@ -125,63 +125,8 @@ def makeSpinFiles(CLIM_name,timeres='1D',Tinterp='mean',spin_date_st = 1980.0, s else: #CLIM_name is not a pickle, it is the dataframe being passed df_CLIM = CLIM_name -# <<<<<<< HEAD if not SEB: - # drn = {'TS':'TSKIN'} #customize this to change your dataframe column names to match the required inputs - # try: - # df_CLIM['RAIN'] = df_CLIM['PRECTOT'] - df_CLIM['PRECSNO'] - # df_CLIM['BDOT'] = df_CLIM['PRECSNO'] + df_CLIM['EVAP'] - # except: - # pass - # df_CLIM.rename(mapper=drn,axis=1,inplace=True) - # try: - # df_CLIM.drop(['EVAP','PRECTOT','PRECSNO'],axis=1,inplace=True) - # except: - # pass - # l1 = df_CLIM.columns.values.tolist() - # l2 = ['SMELT','BDOT','RAIN','TSKIN','SRHO'] - # notin = list(np.setdiff1d(l1,l2)) - # df_CLIM.drop(notin,axis=1,inplace=True) - # # df_BDOT = pd.DataFrame(df_CLIM.BDOT) - # df_TS = pd.DataFrame(df_CLIM.TSKIN) - # res_dict_all = {'SMELT':'sum','BDOT':'sum','RAIN':'sum','TSKIN':'mean','SRHO':'mean'} # resample type for all possible variables - # res_dict = {key:res_dict_all[key] for key in df_CLIM.columns} # resample type for just the data types in df_CLIM - - # # df_BDOT_re = df_BDOT.resample(timeres).sum() - # if Tinterp == 'mean': - # df_TS_re = df_TS.resample(timeres).mean() - # elif Tinterp == 'effective': - # df_TS_re = df_TS.resample(timeres).apply(effectiveT) - # elif Tinterp == 'weighted': - # df_TS_re = pd.DataFrame(data=(df_BDOT.BDOT*df_TS.TSKIN).resample(timeres).sum()/(df_BDOT.BDOT.resample(timeres).sum()),columns=['TSKIN']) - # # pass - - # df_CLIM_re = df_CLIM.resample(timeres).agg(res_dict) - # df_CLIM_re.TSKIN = df_TS_re.TSKIN - # df_CLIM_ids = list(df_CLIM_re.columns) - - # df_CLIM_re['decdate'] = [toYearFraction(qq) for qq in df_CLIM_re.index] - # df_CLIM_re = df_CLIM_re.fillna(method='pad') - - # # df_TS_re['decdate'] = [toYearFraction(qq) for qq in df_TS_re.index] - # # df_BDOT_re['decdate'] = [toYearFraction(qq) for qq in df_BDOT_re.index] - # # df_TS_re = df_TS_re.fillna(method='pad') - - # stepsperyear = 1/(df_CLIM_re.decdate.diff().mean()) - - # BDOT_mean_IE = (df_CLIM_re['BDOT']*stepsperyear/917).mean() - # T_mean = (df_TS_re['TSKIN']).mean() - # print(BDOT_mean_IE) - # print(T_mean) - - # hh = np.arange(0,501) - # age, rho = hla.hl_analytic(350,hh,T_mean,BDOT_mean_IE) - # if not desired_depth: - # desired_depth = hh[np.where(rho>=rho_bottom)[0][0]] - # depth_S1 = hh[np.where(rho>=550)[0][0]] - # depth_S2 = hh[np.where(rho>=750)[0][0]] - # ======= drn = {'TS':'TSKIN','EVAP':'SUBLIM'} #customize this to change your dataframe column names to match the required inputs try: df_CLIM['RAIN'] = df_CLIM['PRECTOT'] - df_CLIM['PRECSNO'] @@ -249,51 +194,6 @@ def makeSpinFiles(CLIM_name,timeres='1D',Tinterp='mean',spin_date_st = 1980.0, s depth_S1 = desired_depth * 0.5 depth_S2 = desired_depth * 0.75 - # #### Make spin up series ### - # RCI_length = spin_date_end-spin_date_st+1 - # num_reps = int(np.round(desired_depth/BDOT_mean_IE/RCI_length)) - # years = num_reps*RCI_length - # sub = np.arange(-1*years,0,RCI_length) - # startyear = int(df_CLIM_re.index[0].year + sub[0]) - # startmonth = df_CLIM_re.index[0].month - # startday = df_CLIM_re.index[0].day - # startstring = '{}/{}/{}'.format(startday,startmonth,startyear) - - # msk = df_CLIM_re.decdate.values>>>>>> master - # else: - # desired_depth = desired_depth - # depth_S1 = desired_depth * 0.5 - # depth_S2 = desired_depth * 0.75 - #### Make spin up series ### RCI_length = spin_date_end-spin_date_st+1 num_reps = int(np.round(desired_depth/BDOT_mean_IE/RCI_length)) @@ -380,7 +280,7 @@ def makeSpinFiles(CLIM_name,timeres='1D',Tinterp='mean',spin_date_st = 1980.0, s hh = np.arange(0,501) age, rho = hla.hl_analytic(350,hh,T_mean,BDOT_mean_IE) if not desired_depth: - desired_depth = hh[np.where(rho>=916)[0][0]] + desired_depth = hh[np.where(rho>=rho_bottom)[0][0]] depth_S1 = hh[np.where(rho>=550)[0][0]] depth_S2 = hh[np.where(rho>=750)[0][0]] else: diff --git a/CFM_main/SEB.py b/CFM_main/SEB.py index 08fbf89..044d144 100644 --- a/CFM_main/SEB.py +++ b/CFM_main/SEB.py @@ -60,6 +60,7 @@ def __init__(self,config,climateTS,start_ind): self.SW_d = climateTS['SW_d'][start_ind:] self.LW_d = climateTS['LW_d'][start_ind:] self.ALBEDO = climateTS['ALBEDO'][start_ind:] + self.ALBEDO[np.isnan(self.ALBEDO)] = np.nanmean(self.ALBEDO) self.T2m = climateTS['T2m'][start_ind:] self.TSKIN = climateTS['TSKIN'][start_ind:] if (('EVAP' in climateTS.keys()) and ('SUBLIM' in climateTS.keys())): @@ -87,8 +88,170 @@ def __init__(self,config,climateTS,start_ind): self.emissivity_air = 1 self.emissivity_snow = 0.98 - # self.D_sh = 15 # Sensible heat flux coefficient, Born et al. 2019 [W m^-2 K^-1] + # self.D_sh = 15 # Sensible heat flux coefficient, Born et al. 2019 [W m^-2 K^-1] + def SEB_fqs(self,PhysParams,iii,T_old): + ''' + Calculate surface energy using fqs solver + Positive fluxes are into the surface, negative are out + SEBparams: mass, Tz, dt + ''' + + Tz = PhysParams['Tz'] + mass = PhysParams['mass'] + dt = PhysParams['dt'] # [s] + Tguess = self.T2m[iii] + dz = PhysParams['dz'] + z = PhysParams['z'] + mtime = PhysParams['mtime'] + + T_rain = np.max((self.T2m[iii],T_MELT)) + # Qrain_i = 0 + + rain_mass = self.RAIN[iii] * RHO_I / S_PER_YEAR * dt #[kg] of rain at this timestep + Qrain_i = CP_W * rain_mass * (T_rain - T_MELT) # Assume rain temperature is air temp, Hock 2005, eq 19 + #latent heat for rain falling on top of cold snow should be handled in melt.py + + Q_SW_net = self.SW_d[iii] * (1-self.ALBEDO[iii]) + # Q_LW_d = self.SBC * (self.emissivity_air * self.T2m[iii]**4) + Q_LW_d = self.emissivity_air * self.LW_d[iii] + + i_GL = np.where(z>=1)[0][0] + z_GL = z[i_GL] + m_GL = np.cumsum(mass)[i_GL] + T_GL = np.cumsum(mass*Tz)[i_GL]/m_GL + rho_GL = m_GL/z_GL + K_ice = 9.828 * np.exp(-0.0057 * T_GL) #[W/m/K] + + K_GL = K_ice * (rho_GL/RHO_I) ** (2 - 0.5 * (rho_GL/RHO_I)) + + G = (K_GL * (Tz[i_GL] - Tz[0])/z_GL) # estimated temperature flux in firn due to temperature gradient + # G = 0 + + + + TL_thick = 0.01 # thickness of snow/firn "Top Layer" that energy goes into. Reducing results in higher melt. + iTL = np.where(z>=TL_thick)[0][0] + + for kk in range(10): + dTL = np.cumsum(dz)[iTL] + + m = np.cumsum(mass)[iTL] #mass of the TL + + TTL = np.cumsum(mass*Tz)[iTL]/m # mean temperature of top X cm (weighted mean) + cold_content_TL = CP_I * m * (T_MELT - TTL) # cold content [J], positive quantity if T0))].real) + + if Tsurface>=273.15: + Tsurface = 273.15 + meltmass = (Qnet - self.SBC*273.15**4) * dt / LF_I #multiply by dt to put in units per time step + # do not need to subtract cold content to calculate cold content b/c Q_melt = sum(energies), Q_melt=0 if energy can be balanced, i.e. sum(energies)=0 + # melt_mass = (Qnet - self.SBC*273.15**4) * dt / LF_I + ### meltmass has units [kg/m2/s] + else: + meltmass = 0 + + if meltmass<=m: #if the melt mass is greater than the mass of the TL layer, we need a thicker TL because the next layer could be below freezing, and it needs to warm before melting + break + else: + iTL = np.where(np.cumsum(mass)>=meltmass)[0][0] + + Tz[0:iTL+1] = Tsurface + + return Tsurface, Tz, meltmass, self.TSKIN[iii] + ############################ + ### end SEB_fqs + ############################ + + def SEB_loop(self,PhysParams,iii,T_old): + ''' + Calculate surface energy using fqs solver + Positive fluxes are into the surface, negative are out + SEBparams: mass, Tz, dt + ''' + + Tz = PhysParams['Tz'] + mass = PhysParams['mass'] + dt = PhysParams['dt'] # [s] + Tguess = self.T2m[iii] + dz = PhysParams['dz'] + z = PhysParams['z'] + + T_rain = np.max((self.T2m[iii],T_MELT)) + # Qrain_i = 0 + + rain_mass = self.RAIN[iii] * RHO_I / S_PER_YEAR * dt #[kg] of rain at this timestep + Qrain_i = CP_W * rain_mass * (T_rain - T_MELT) # Assume rain temperature is air temp, Hock 2005, eq 19 + #latent heat for rain falling on top of cold snow should be handled in melt.py + + Q_SW_net = self.SW_d[iii] * (1-self.ALBEDO[iii]) + # Q_LW_d = self.SBC * (self.emissivity_air * self.T2m[iii]**4) + Q_LW_d = self.emissivity_air * self.LW_d[iii] + + TL_thick = 0.1 # thickness of snow/firn "Top Layer" that energy goes into. Reducing results in higher melt. + iTL = np.where(z>=TL_thick)[0][0] + dTL = z[iTL] + + i_GL = np.where(z>=1)[0][0] + z_GL = z[i_GL] + m_GL = np.cumsum(mass)[i_GL] + T_GL = np.cumsum(mass*Tz)[i_GL]/m_GL + rho_GL = m_GL/z_GL + K_ice = 9.828 * np.exp(-0.0057 * T_GL) #[W/m/K] + + K_GL = K_ice * (rho_GL/RHO_I) ** (2 - 0.5 * (rho_GL/RHO_I)) + G = (K_GL * (Tz[i_GL] - Tz[0])/z_GL) # estimated temperature flux in firn due to temperature gradient + + m = np.cumsum(mass)[iTL] #mass of the top layer + TTL = np.cumsum(mass*Tz)[iTL]/m # mean temperature of top layer (weighted mean) + cold_content_TL = CP_I * m * (T_MELT - TTL) # cold content [J], positive quantity if T=273.15: + Tsurface = 273.15 + meltmass = (Q_sum - self.SBC*273.15**4) * dt / LF_I #*dt #multiply by dt to put in units per day + # do not need to subtract cold content to calculate cold content b/c Q_melt = sum(energies), Q_melt=0 if energy can be balanced, i.e. sum(energies)=0 + # melt_mass = (Qnet - self.SBC*273.15**4) * dt / LF_I + ### meltmass has units [kg/m2/s] + else: + meltmass = 0 + + Tz[0:iTL+1] = Tsurface + + return Tsurface, Tz, meltmass + ############################ + ### end SEB_loop + ############################ def SEB(self, PhysParams,iii,T_old): ''' @@ -132,10 +295,27 @@ def SEB(self, PhysParams,iii,T_old): # Q_LW_d = self.SBC * (self.emissivity_air * self.T2m[iii]**4) Q_LW_d = self.emissivity_air * self.LW_d[iii] - i10cm = np.where(z>=0.1)[0][0] - d10cm = z[i10cm] - G = (0.3*(Tz[i10cm] - Tz[0])/d10cm) - m = mass[i10cm] + + TL_thick = 0.1 # thickness of snow/firn "Top Layer" that energy goes into. Reducing results in higher melt. + iTL = np.where(z>=TL_thick)[0][0] + dTL = z[iTL] + + i_GL = np.where(z>=1)[0][0] + z_GL = z[i_GL] + m_GL = np.cumsum(mass)[i_GL] + T_GL = np.cumsum(mass*Tz)[i_GL]/m_GL + rho_GL = m_GL/z_GL + K_ice = 9.828 * np.exp(-0.0057 * T_GL) #[W/m/K] + + K_GL = K_ice * (rho_GL/RHO_I) ** (2 - 0.5 * (rho_GL/RHO_I)) + + G = (K_GL * (Tz[i_GL] - Tz[0])/z_GL) # estimated temperature flux in firn due to temperature gradient + + iTL = np.where(z>=0.1)[0][0] + dTL = z[iTL] + + m = np.cumsum(mass)[iTL] # this was on staging + # G=0 Qnet = Q_SW_net + Q_LW_d + self.QH[iii] + self.QL[iii] + Qrain_i + G @@ -175,138 +355,344 @@ def SEB(self, PhysParams,iii,T_old): else: Tsurface = 273.15 melt_mass = (Qnet - self.SBC*273.15**4) * dt / LF_I - - # df12_daily['meltvol'] = meltvold - # df12_daily['Tcalc'] = Tcalcd - - # if Tsurface>=T_MELT: - # Q_LW_me_old = (self.SBC * (self.emissivity_air * self.T2m[iii]**4 - self.emissivity_snow * (T_MELT)**4)) - # Q_LW_up = (self.SBC * self.emissivity_snow * (T_MELT)**4) - # Q_LW_me = Q_LW_d - Q_LW_up - - # E_melt = (Q_SW_net + Q_LW_me + self.QH[iii] + self.QL[iii] + Qrain_i + G) #[W/m2] - - # Tsurface = T_MELT - # melt_mass = E_melt/LF_I * dt - - ################### - # if (E_melt*dt) < cold_content[0]: - # Tsurface = Tz[0] + E_melt/(CP_I * mass[0]) * dt - # melt_mass = 0 - # else: - # Tsurface=T_MELT - # energy_excess = E_melt*dt - (CP_I * mass[0] * (T_MELT - Tz[0])) # [J] - # melt_mass = energy_excess/LF_I - ################### - - ################### - # E_iter = E_melt.copy() - # E_iter = E_melt*dt - # kk = 0 - # melt_mass = 0 - - # while E_iter > 0: - # if E_iter <= cold_content[kk]: # all energy warms/cools this firn; no melt in it - # Tz[kk] = Tz[kk] + E_iter / (CP_I * mass[kk]) #* dt - # E_iter = 0 - - # elif E_iter > cold_content[kk]: - # Tz[kk] = T_MELT - # E_melt1 = E_iter - cold_content[kk] #energy available for melting - # melt_mass_pot = E_melt1/LF_I#*dt # mass that can be melted with the energy available - - # if melt_mass_pot <= mass[kk]: #less than the entire node melts - # melt_mass += melt_mass_pot - # E_iter = 0 - - # else: #entire node melts - # melt_mass += mass[kk] - # E_iter -= mass[kk]*LF_I#*dt - - # kk+=1 - - # Tsurface = T_MELT - ###################### - - - # else: - # melt_mass = 0 Tz[0] = Tsurface if melt_mass<0: melt_mass = 0 - # print(iii,Tsurface,melt_mass) - return Tsurface, Tz, melt_mass + +# Fast Quartic Solver: analytically solves quartic equations (needed to calculate melt) +# Takes methods from fqs package (@author: NKrvavica) +# full documentation: https://github.com/NKrvavica/fqs/blob/master/fqs.py - # def SEB_direct(self, PhysParams,iii): - # ''' - # Calculate the surface energy budget - # Positive fluxes are into the surface, negative are out - # SEBparams: mass, Tz, dt - # ''' +def single_quadratic(a0, b0, c0): + ''' + Analytical solver for a single quadratic equation + ''' + a, b = b0 / a0, c0 / a0 - # # max_iter = 20 - # # ijk = 0 - # Tz = PhysParams['Tz'] - # mass = PhysParams['mass'] - # dt = PhysParams['dt'] - # print('dt',dt) - - # T_rain = np.max((self.T2m[iii],T_MELT)) - # Qrain_i = RHO_W_KGM * CP_W * self.RAIN[iii] * (T_rain - T_MELT) # Assume rain temperature is air temp, Hock 2005, eq 19 - # # heat for rain falling on top of cold snow should be handled in melt.py + # Some repating variables + a0 = -0.5*a + delta = a0*a0 - b + sqrt_delta = np.sqrt(delta) - # Q_SW_i = self.SW_d[iii] * (1-self.ALBEDO[iii]) + # Roots + r1 = a0 - sqrt_delta + r2 = a0 + sqrt_delta - # emissivity_air = emissivity_snow = 1 - # Q_LW_i = self.SBC * (emissivity_air * self.T2m[iii]**4 - emissivity_snow * (Tz[0])**4) + return r1, r2 - # E_net = Q_SW_i + Q_LW_i + self.QH[iii] + self.QL[iii] + Qrain_i + self.G[iii] - # print('E_net',E_net) - # # time dimension for conversion? - # # dTz[0] = E_net * dt[iii] / CP_I / mass - # cold_content = CP_I * mass * (T_MELT - Tz) # cold content [J] - # print('cold_content',cold_content) - # print('mass',mass[0:5]) +def single_cubic(a0, b0, c0, d0): + ''' + Analytical closed-form solver for a single cubic equation + ''' + a, b, c = b0 / a0, c0 / a0, d0 / a0 + + # Some repeating constants and variables + third = 1./3. + a13 = a*third + a2 = a13*a13 + sqr3 = np.sqrt(3) + + # Additional intermediate variables + f = third*b - a2 + g = a13 * (2*a2 - b) + c + h = 0.25*g*g + f*f*f + + def cubic_root(x): + ''' Compute cubic root of a number while maintaining its sign''' + if x.real >= 0: + return x**third + else: + return -(-x)**third + + if f == g == h == 0: + r1 = -cubic_root(c) + return r1, r1, r1 + + elif h <= 0: + j = np.sqrt(-f) + k = np.arccos(-0.5*g / (j*j*j)) + m = np.cos(third*k) + n = sqr3 * np.sin(third*k) + r1 = 2*j*m - a13 + r2 = -j * (m + n) - a13 + r3 = -j * (m - n) - a13 + return r1, r2, r3 + + else: + sqrt_h = np.sqrt(h) + S = cubic_root(-0.5*g + sqrt_h) + U = cubic_root(-0.5*g - sqrt_h) + S_plus_U = S + U + S_minus_U = S - U + r1 = S_plus_U - a13 + r2 = -0.5*S_plus_U - a13 + S_minus_U*sqr3*0.5j + r3 = -0.5*S_plus_U - a13 - S_minus_U*sqr3*0.5j + return r1, r2, r3 + + + +def single_cubic_one(a0, b0, c0, d0): + ''' + Analytical closed-form solver for a single cubic equation + ''' + a, b, c = b0 / a0, c0 / a0, d0 / a0 - # E_iter = E_net.copy() - # kk = 0 - # melt_mass = 0 - - # while E_iter > 0: - # if E_iter <= cold_content[kk]: # all energy warms/cools the firn; no melt - # Tz[kk] = Tz[kk] + E_net / (CP_I * mass[kk]) * dt - # E_iter = 0 + # Some repeating constants and variables + third = 1./3. + a13 = a*third + a2 = a13*a13 - # elif E_net > cold_content[kk]: - # Tz[kk] = T_MELT - # E_melt = E_net - cold_content[kk] #energy available for melting - # melt_mass_pot = E_melt/LF_I # mass that can be melted + # Additional intermediate variables + f = third*b - a2 + g = a13 * (2*a2 - b) + c + h = 0.25*g*g + f*f*f - # if melt_mass_pot <= mass[kk]: #less than the entire node melts - # melt_mass += melt_mass_pot - # E_iter = 0 + def cubic_root(x): + ''' Compute cubic root of a number while maintaining its sign + ''' + if x.real >= 0: + return x**third + else: + return -(-x)**third + + if f == g == h == 0: + return -cubic_root(c) + + elif h <= 0: + j = np.sqrt(-f) + k = np.arccos(-0.5*g / (j*j*j)) + m = np.cos(third*k) + return 2*j*m - a13 + + else: + sqrt_h = np.sqrt(h) + S = cubic_root(-0.5*g + sqrt_h) + U = cubic_root(-0.5*g - sqrt_h) + S_plus_U = S + U + return S_plus_U - a13 - # else: #entire node melts - # melt_mass += mass[kk] - # E_iter -= mass[kk]*LF_I +def single_quartic(a0, b0, c0, d0, e0): + ''' + Analytical closed-form solver for a single quartic equation + ''' + a, b, c, d = b0/a0, c0/a0, d0/a0, e0/a0 - # kk+=1 + # Some repeating variables + a0 = 0.25*a + a02 = a0*a0 - # print(Tz[0:5]) - # input('SEB') - # return Tz, melt_mass + # Coefficients of subsidiary cubic euqtion + p = 3*a02 - 0.5*b + q = a*a02 - b*a0 + 0.5*c + r = 3*a02*a02 - b*a02 + c*a0 - d + # One root of the cubic equation + z0 = single_cubic_one(1, p, r, p*r - 0.5*q*q) + # Additional variables + s = np.sqrt(2*p + 2*z0.real + 0j) + if s == 0: + t = z0*z0 + r + else: + t = -q / s + # Compute roots by quadratic equations + r0, r1 = single_quadratic(1, s, z0 + t) + r2, r3 = single_quadratic(1, -s, z0 - t) + return r0 - a0, r1 - a0, r2 - a0, r3 - a0 +def multi_quadratic(a0, b0, c0): + ''' + Analytical solver for multiple quadratic equations + ''' + a, b = b0 / a0, c0 / a0 + + # Some repating variables + a0 = -0.5*a + delta = a0*a0 - b + sqrt_delta = np.sqrt(delta + 0j) + + # Roots + r1 = a0 - sqrt_delta + r2 = a0 + sqrt_delta + + return r1, r2 + + +def multi_cubic(a0, b0, c0, d0, all_roots=True): + ''' + Analytical closed-form solver for multiple cubic equations + ''' + a, b, c = b0 / a0, c0 / a0, d0 / a0 + + # Some repeating constants and variables + third = 1./3. + a13 = a*third + a2 = a13*a13 + sqr3 = np.sqrt(3) + + # Additional intermediate variables + f = third*b - a2 + g = a13 * (2*a2 - b) + c + h = 0.25*g*g + f*f*f + + # Masks for different combinations of roots + m1 = (f == 0) & (g == 0) & (h == 0) # roots are real and equal + m2 = (~m1) & (h <= 0) # roots are real and distinct + m3 = (~m1) & (~m2) # one real root and two complex + + def cubic_root(x): + ''' Compute cubic root of a number while maintaining its sign + ''' + root = np.zeros_like(x) + positive = (x >= 0) + negative = ~positive + root[positive] = x[positive]**third + root[negative] = -(-x[negative])**third + return root + + def roots_all_real_equal(c): + ''' Compute cubic roots if all roots are real and equal + ''' + r1 = -cubic_root(c) + if all_roots: + return r1, r1, r1 + else: + return r1 + + def roots_all_real_distinct(a13, f, g, h): + ''' Compute cubic roots if all roots are real and distinct + ''' + j = np.sqrt(-f) + k = np.arccos(-0.5*g / (j*j*j)) + m = np.cos(third*k) + r1 = 2*j*m - a13 + if all_roots: + n = sqr3 * np.sin(third*k) + r2 = -j * (m + n) - a13 + r3 = -j * (m - n) - a13 + return r1, r2, r3 + else: + return r1 + + def roots_one_real(a13, g, h): + ''' Compute cubic roots if one root is real and other two are complex + ''' + sqrt_h = np.sqrt(h) + S = cubic_root(-0.5*g + sqrt_h) + U = cubic_root(-0.5*g - sqrt_h) + S_plus_U = S + U + r1 = S_plus_U - a13 + if all_roots: + S_minus_U = S - U + r2 = -0.5*S_plus_U - a13 + S_minus_U*sqr3*0.5j + r3 = -0.5*S_plus_U - a13 - S_minus_U*sqr3*0.5j + return r1, r2, r3 + else: + return r1 + + # Compute roots + if all_roots: + roots = np.zeros((3, len(a))).astype(complex) + roots[:, m1] = roots_all_real_equal(c[m1]) + roots[:, m2] = roots_all_real_distinct(a13[m2], f[m2], g[m2], h[m2]) + roots[:, m3] = roots_one_real(a13[m3], g[m3], h[m3]) + else: + roots = np.zeros(len(a)) # .astype(complex) + roots[m1] = roots_all_real_equal(c[m1]) + roots[m2] = roots_all_real_distinct(a13[m2], f[m2], g[m2], h[m2]) + roots[m3] = roots_one_real(a13[m3], g[m3], h[m3]) + + return roots + + +def multi_quartic(a0, b0, c0, d0, e0): + ''' + Analytical closed-form solver for multiple quartic equations + ''' + a, b, c, d = b0/a0, c0/a0, d0/a0, e0/a0 + + # Some repeating variables + a0 = 0.25*a + a02 = a0*a0 + + # Coefficients of subsidiary cubic euqtion + p = 3*a02 - 0.5*b + q = a*a02 - b*a0 + 0.5*c + r = 3*a02*a02 - b*a02 + c*a0 - d + + # One root of the cubic equation + z0 = multi_cubic(1, p, r, p*r - 0.5*q*q, all_roots=False) + + # Additional variables + s = np.sqrt(2*p + 2*z0.real + 0j) + t = np.zeros_like(s) + mask = (s == 0) + t[mask] = z0[mask]*z0[mask] + r[mask] + t[~mask] = -q[~mask] / s[~mask] + + # Compute roots by quadratic equations + r0, r1 = multi_quadratic(1, s, z0 + t) - a0 + r2, r3 = multi_quadratic(1, -s, z0 - t) - a0 + + return r0, r1, r2, r3 + + +def cubic_roots(p): + ''' + A caller function for a fast cubic root solver (3rd order polynomial). + ''' + # Convert input to array (if input is a list or tuple) + p = np.asarray(p) + + # If only one set of coefficients is given, add axis + if p.ndim < 2: + p = p[np.newaxis, :] + + # Check if four coefficients are given + if p.shape[1] != 4: + raise ValueError('Expected 3rd order polynomial with 4 ' + 'coefficients, got {:d}.'.format(p.shape[1])) + + if p.shape[0] < 100: + roots = [single_cubic(*pi) for pi in p] + return np.array(roots) + else: + roots = multi_cubic(*p.T) + return np.array(roots).T + + +def quartic_roots(p): + ''' + A caller function for a fast quartic root solver (4th order polynomial). + p[0]*x^4 + p[1]*x^3 + p[2]*x^2 + p[3]*x + p[4] = 0 + ''' + # Convert input to an array (if input is a list or tuple) + p = np.asarray(p) + + # If only one set of coefficients is given, add axis + if p.ndim < 2: + p = p[np.newaxis, :] + + # Check if all five coefficients are given + if p.shape[1] != 5: + raise ValueError('Expected 4th order polynomial with 5 ' + 'coefficients, got {:d}.'.format(p.shape[1])) + + if p.shape[0] < 100: + roots = [single_quartic(*pi) for pi in p] + return np.array(roots) + else: + roots = multi_quartic(*p.T) + return np.array(roots).T ''' diff --git a/CFM_main/diffusion.py b/CFM_main/diffusion.py index a50dc97..db02efa 100644 --- a/CFM_main/diffusion.py +++ b/CFM_main/diffusion.py @@ -125,6 +125,7 @@ def heatDiff(self,iii): Gamma_P = K_firn + tot_rho = self.rho c_vol = self.rho * c_firn @@ -207,9 +208,11 @@ def enthalpyDiff(self,iii): lwc_old = self.LWC.copy() - phi_ret, g_liq, count, iterdiff, g_sol = transient_solve_EN(z_edges_vec, z_P_vec, nt, self.dt[iii], K_eff, phi_0, nz_P, nz_fv, phi_s, tot_rho, c_vol, self.LWC, self.mass, self.dz,ICT,self.rho,iii) - + phi_ret, g_liq, count, iterdiff,g_sol = transient_solve_EN(z_edges_vec, z_P_vec, nt, self.dt[iii], K_eff, phi_0, nz_P, nz_fv, phi_s, tot_rho, c_vol, self.LWC, self.mass, self.dz,ICT,self.rho,iii) + LWC_ret = g_liq * self.dz + # self.LWC = g_liq * vol_tot + delta_mass_liq = mass_liq - (LWC_ret * RHO_W_KGM) dml_sum = 0.0 @@ -238,8 +241,10 @@ def enthalpyDiff(self,iii): if np.any(delta_mass_liq<0): if np.any(np.abs(delta_mass_liq[delta_mass_liq<0])>1e-7): print('------') + print('If you are seeing this message there was a liquid mass gain in diffusion.') print('Please email maxstev@umd.edu so I can fix it.') + dml_sum = np.sum(delta_mass_liq[delta_mass_liq<0]) delta_mass_liq = np.maximum(delta_mass_liq,0) # fix for numerical instabilities with small time steps. @@ -255,8 +260,8 @@ def enthalpyDiff(self,iii): ############################## def heatDiff_highC(self,iii): - ''' + ''' IN DEVELOPMENT One way of dealing with liquid water in the firn @@ -319,11 +324,8 @@ def heatDiff_highC(self,iii): ############################## def heatDiff_Teff(self,iii): - ''' - + ''' IN DEVELOPMENT - - artificially set the temperature of volumes with liquid water to be higher than T_melt ''' @@ -395,6 +397,7 @@ def heatDiff_Teff(self,iii): def heatDiff_LWCcorr(self,iii, iters,correct_therm_prop): ''' + IN DEVELOPMENT just run the heat diffusion as normal and then balance energy. diff --git a/CFM_main/example.json b/CFM_main/example.json index 0bafb60..af346a9 100644 --- a/CFM_main/example.json +++ b/CFM_main/example.json @@ -8,7 +8,7 @@ "InputFileNameStrain": "example_STRAIN.csv", "InputFileNameSublim": "example_SUBLIM.csv", "InputFileNameRain": "example_RAIN.csv", -"resultsFolder": "CFMoutput_example/df", +"resultsFolder": "CFMoutput_example/csv", "initfirnFile": "example_firndata.csv", "initprofile": false, "input_type": "csv", @@ -102,7 +102,8 @@ "DirectRunoff": 0.0, "RunoffZuoOerlemans": false, "Slope": 0.1, -"SUBLIM": true +"SUBLIM": true, +"keep_firnthickness": false } diff --git a/CFM_main/firn_density_nospin.py b/CFM_main/firn_density_nospin.py index 804a554..a26637f 100755 --- a/CFM_main/firn_density_nospin.py +++ b/CFM_main/firn_density_nospin.py @@ -84,6 +84,7 @@ class FirnDensityNoSpin: (unit: kg m^-3, type: array of floats) : bdot_mean: mean accumulation over the lifetime of each parcel (units are m I.E. per year) + : sublim: sublimation/deposition. Negative means sublimation, positive means depostion :returns D_surf: diffusivity tracker (unit: ???, type: array of floats) @@ -306,8 +307,10 @@ def __init__(self, configName, climateTS = None, NewSpin = False): self.MELT = True try: self.LWC = initLWC[1:] + self.LWC_init = np.sum(self.LWC) except: self.LWC = np.zeros_like(self.z) + self.LWC_init = 0 self.PLWC_mem = np.zeros_like(self.z) #VV keep track of water content that was in PFdom if 'RAIN' not in self.c: @@ -337,6 +340,11 @@ def __init__(self, configName, climateTS = None, NewSpin = False): self.LWC = np.zeros_like(self.z) self.PLWC_mem = np.zeros_like(self.z) self.c['RAIN'] = False + + if 'keep_firnthickness' not in self.c: + self.c['keep_firnthickness'] = False + print('Note: keep_firnthickness is not in your .json (see changelog for v2.2.0)') + print('Defaulting to false (old CFM behavior)') ##################### ##################### @@ -450,14 +458,14 @@ def __init__(self, configName, climateTS = None, NewSpin = False): susf = interpolate.interp1d(input_year_sublim,input_sublim,int_type,fill_value='extrapolate') # interpolation function self.sublim = susf(self.modeltime) # [m ice eq per year] if np.any(self.sublim>0): - ipSUB = np.where(self.sublim>0) + ipSUB = np.where(self.sublim>0)[0] self.bdot[ipSUB] = self.bdot[ipSUB] + self.sublim[ipSUB] self.sublim[ipSUB] = 0 self.sublimSec = self.sublim / S_PER_YEAR / (S_PER_YEAR/self.dt) # sublimation at each time step (meters i.e. per second). gets multiplied by S_PER_YEAR later. (sort of hacky, I know) self.bdotSec = self.bdot / S_PER_YEAR / (S_PER_YEAR/self.dt) # accumulation at each time step (meters i.e. per second). gets multiplied by S_PER_YEAR later. (sort of hacky, I know) - try: #Rolling mean average surface temperature and accumulation rate (vector) - # (i.e. the long-term average climate) + try: # Rolling mean average surface temperature and accumulation rate (vector) + # (i.e. the long-term average climate) Nyears = 10 #number of years to average for T_mean NN = int(np.mean(S_PER_YEAR/self.dt)*Nyears) # NN = int(self.c['stpsPerYear']*Nyears) @@ -620,7 +628,7 @@ def __init__(self, configName, climateTS = None, NewSpin = False): ### transform mass in meters ice equiv -> divide by age(in sec) [m/s] -> multiply by years per step and by steps per year (cancels) -> multiply by secperyear -> [mIE/yr] ### for surf layer -> mass in mIE is only multiplied by steps per year: if 1 stp/yr,mean acc is the mass of surf layer; if 2 stps/yr,mean acc is 2* what has been accumulated over the last step, etc. ####################### - + ####################### if self.c['manualT']: self.Tz = np.interp(self.z, self.manualT_dep, init_Tz) @@ -653,7 +661,7 @@ def __init__(self, configName, climateTS = None, NewSpin = False): # Load strain rate input if self.c['horizontal_divergence'] or self.c['strain_softening']: self.eps_eff_hor_2, self.eps_divergence, self.c = load_strain(self, spin=False) - + ### initial grain growth (if specified in config file) if self.c['physGrain']: initr2 = read_init(self.c['resultsFolder'], self.c['spinFileName'], 'r2Spin') @@ -856,6 +864,13 @@ def __init__(self, configName, climateTS = None, NewSpin = False): self.MOutputs = ModelOutputs(self.c,MOd,TWlen, init_time, len(self.dz)) ### self.MOutputs is a class instance + + self.na_count = 0 + self.na_sum = 0 + self.melt_sum = 0 + self.acc_sum = 0 + self.dsdz_sum = 0 + self.ddz_bdot = 0 ###################################### ###################################### @@ -882,7 +897,7 @@ def time_evolve(self): runoff2check = np.zeros(self.stp) meltvol2check = np.zeros(self.stp) dml2check = np.zeros(self.stp) - + self.mismatch = 0 #################################### ##### START TIME-STEPPING LOOP ##### #################################### @@ -890,6 +905,9 @@ def time_evolve(self): print('modeltime',self.modeltime[0],self.modeltime[-1]) for iii in range(self.stp): mtime = self.modeltime[iii] + + zbot_old = self.z[-1] + self.D_surf[iii] = iii # This gives each layer a tracking number that is just iteration number. if iii==1000: if ((self.MELT) and (self.c['liquid']=='darcy')): @@ -902,12 +920,10 @@ def time_evolve(self): if self.c['merging']: # merging may be deprecated (check with VV) if ((self.dz[1] < self.c['merge_min']) or (self.dz[0] < 1e-10)): # Start with surface merging self.dz,self.z,self.gridLen,self.dx,self.rho,self.age,self.LWC,self.PLWC_mem,self.mass,self.mass_sum,self.sigma,self.bdot_mean,\ - self.Dcon,self.T_mean,self.T10m,self.r2 = mergesurf(self,self.c['merge_min'],iii) - # print('merging1 at ', iii) + self.Dcon,self.T_mean,self.T10m,self.r2 = mergesurf(self,self.c['merge_min'],iii) if (np.any(self.dz[2:] < self.c['merge_min'])): # Then merge rest of the firn column self.dz,self.z,self.gridLen,self.dx,self.rho,self.age,self.LWC,self.PLWC_mem,self.mass,self.mass_sum,self.sigma,self.bdot_mean,\ self.Dcon,self.T_mean,self.T10m,self.r2 = mergenotsurf(self,self.c['merge_min'],iii) - # print('merging2 at ', iii) ### dictionary of the parameters that get passed to physics PhysParams = { @@ -1005,23 +1021,37 @@ def time_evolve(self): self.sdz_old = np.sum(self.dz) # old total column thickness (s for sum) self.z_old = np.copy(self.z) self.dz = self.mass / self.rho * self.dx # new dz after compaction + self.z = self.dz.cumsum(axis = 0) + # znew = np.copy(self.z) + self.z = np.concatenate(([0], self.z[:-1])) + + sdz_new = np.sum(self.dz) + dsdz = sdz_new - self.sdz_old + self.dsdz_sum = self.dsdz_sum + dsdz ### Surface energy balance ##### if self.c['SEB']: if iii==0: print('CAUTION: SEB still in beta') - PhysParams.update(dz=self.dz,rho=self.rho) # update dict, will be used for SEB + PhysParams.update(dz=self.dz,rho=self.rho,mtime=mtime) # update dict, will be used for SEB + if iii == 0: T_old = self.Tz[0] else: T_old = self.Ts[iii-1] - self.Ts[iii], self.Tz, melt_mass = self.SEB.SEB(PhysParams,iii,T_old) + + self.Ts[iii], self.Tz, melt_mass, M2TS = self.SEB.SEB_fqs(PhysParams,iii,T_old) + # self.Ts[iii] = self.Tz[0] # set the surface temp to the skin temp calclated by SEB (needed for diffusion module) - self.snowmeltSec[iii] = melt_mass / RHO_I / S_PER_YEAR - self.snowmelt[iii] = self.snowmeltSec[iii] * S_PER_YEAR * (S_PER_YEAR/self.dt[iii]) + self.snowmelt[iii] = melt_mass / RHO_I / self.dt[iii] * S_PER_YEAR # m i.e. per year ([kg/m2/timestep] / [kg/m3] / [s/timestep] * [s/year]) + self.snowmeltSec[iii] = self.snowmelt[iii] / S_PER_YEAR / (S_PER_YEAR/self.dt[iii]) # melt at this time step (mIE/s) + + # self.snowmeltSec[iii] = melt_mass / RHO_I / S_PER_YEAR + # self.snowmelt[iii] = self.snowmeltSec[iii] * S_PER_YEAR * (S_PER_YEAR/self.dt[iii]) # self.snowmeltSec = self.snowmelt / S_PER_YEAR / (S_PER_YEAR/self.dt) # melt for each time step (meters i.e. per second) self.forcing_dict['SMELT'][self.start_ind+iii] = self.snowmelt[iii] + ############################ if self.THist: @@ -1035,13 +1065,27 @@ def time_evolve(self): if self.c['liquid'] == 'bucket': if (self.snowmeltSec[iii]>0) or (np.any(self.LWC > 0.)) or (self.rainSec[iii] > 0.): #i.e. there is water - self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, meltgridtrack, self.refreeze, self.runoff, self.dh_melt = bucket(self,iii) + LWCpre = np.sum(self.LWC * RHO_W_KGM) #mass of liquid before bucket + self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, \ + self.dzn, self.LWC, meltgridtrack, self.refreeze, self.runoff, self.dh_melt = bucket(self,iii) + + # refreeze is [m we] + m2X = self.refreeze + self.runoff + np.sum(self.LWC * RHO_W_KGM) + if ((self.rainSec[iii] == 0) & (LWCpre!=m2X)): + self.mismatch = self.mismatch + (m2X - LWCpre) / RHO_W_KGM + # print('mismatch') + # print(m2X) + # print(LWCpre) + if self.doublegrid==True: # if we use doublegrid -> use the gridtrack corrected for melting self.gridtrack = np.copy(meltgridtrack) + self.meltvol = self.snowmeltSec[iii]*S_PER_YEAR*0.917 #[m w.e.] + else: # Dry firn column and no input of meltwater self.dzn = self.dz[0:self.compboxes] # Not sure this is necessary self.refreeze, self.runoff, self.meltvol, self.dh_melt = 0.,0.,0.,0. + ### end bucket ################## elif self.c['liquid'] == 'darcy': @@ -1099,8 +1143,8 @@ def time_evolve(self): ### end prefsnowpack ################## if self.LWC[-1] > 0.: #VV we don't want to lose water - pass - # print('LWC in last layer that is going to be removed, amount is:',self.LWC[-1]) + # pass + print('LWC in last layer that is going to be removed, amount is:',self.LWC[-1]) # self.LWC[-2] += self.LWC[-1] #VV, flow routine will deal with saturation exceeding 1 # This should never happen if bottom of modelled firn column is at rho >= 830 # self.LWC = np.concatenate(([0], self.LWC[:-1])) @@ -1148,7 +1192,7 @@ def time_evolve(self): for isotope in self.c['iso']: self.Isoz[isotope], self.Iso_sig2_z[isotope] = self.Isotopes[isotope].isoDiff(IsoParams,iii) ### new box gets added on within isoDiff function - #################### + #################### self.sdz_new = np.sum(self.dz) #total column thickness after densification, melt, horizontal strain, before new snow added @@ -1161,25 +1205,39 @@ def time_evolve(self): ### update model grid, mass, stress, and mean accumulation rate ### If SEB, Ts was set in the SEB module - do not update if there is not snow that has temperature T2m. - - if self.bdotSec[iii]>0: # there is accumulation at this time step + + zz1 = self.z.copy() + bd_flag = None + + if self.bdotSec[iii]>0: # there is accumulation at this time step self.age = np.concatenate(([0], self.age[:-1])) + self.dt[iii] self.dzNew = self.bdotSec[iii] * RHO_I / self.rhos0[iii] * S_PER_YEAR self.dh_acc = self.dzNew + + dz_bot_old = self.dz[-1] + z_bot_old = self.z[-1] + self.dz = np.concatenate(([self.dzNew], self.dz[:-1])) self.z = self.dz.cumsum(axis = 0) znew = np.copy(self.z) self.z = np.concatenate(([0], self.z[:-1])) self.rho = np.concatenate(([self.rhos0[iii]], self.rho[:-1])) + dz_bot_new = self.dz[-1] + z_bot_new = self.z[-1] + + dzb_diff = dz_bot_new - dz_bot_old + # self.ddz_bdot = self.ddz_bdot + dzdiff + if self.c['physGrain']: # update grain radius r2surface = FirnPhysics(PhysParams).surfacegrain() #grain size for new surface layer self.r2 = np.concatenate(([r2surface], self.r2[:-1])) if not self.c['manualT']: # If SEB, the new snow layer will be T2m if self.c['SEB']: - newSnowT = np.max((self.T2m[iii],T_MELT)) + newSnowT = np.min((self.T2m[iii],T_MELT)) else: newSnowT = self.Ts[iii] + self.Tz = np.concatenate((np.array([newSnowT],float), self.Tz[:-1])) self.Dcon = np.concatenate(([self.D_surf[iii]], self.Dcon[:-1])) @@ -1191,17 +1249,26 @@ def time_evolve(self): self.LWC = np.concatenate(([0], self.LWC[:-1])) else: # no accumulation during this time step + bd_flag = 'no accumulation' self.age = self.age + self.dt[iii] - self.z = self.dz.cumsum(axis=0) - self.z = np.concatenate(([0],self.z[:-1])) + ddz = self.dz_old - self.dz + self.z[1:] = self.z[1:] - np.cumsum(ddz[0:-1]) + # self.z = self.dz.cumsum(axis=0) + # self.z = np.concatenate(([0],self.z[:-1])) self.dzNew = 0 self.dh_acc = 0 znew = np.copy(self.z) self.compaction = (self.dz_old[0:self.compboxes]-self.dzn) + if ((not self.c['SEB']) or (not self.c['manualT'])): - self.Tz = np.concatenate(([self.Ts[iii]], self.Tz[1:])) + # self.Tz = np.concatenate(([self.Ts[iii]], self.Tz[1:])) + self.Tz[0] = self.Ts[iii] + + self.na_count += 1 + self.na_sum = self.na_sum + (self.z[-1]-zz1[-1]) if self.c['SUBLIM']: + zz1 = self.z.copy() if self.sublimSec[iii]<0: self.mass_sum = self.mass.cumsum(axis = 0) self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, self.PLWC_mem, self.totwatersublim, sublgridtrack, dh_sub = sublim(self,iii) # keeps track of sublimated water for mass conservation @@ -1211,9 +1278,12 @@ def time_evolve(self): if self.doublegrid == True: # gridtrack corrected for sublimation self.gridtrack = np.copy(sublgridtrack) znew = np.copy(self.z) - if not self.c['manualT']: - self.Tz = np.concatenate(([self.Ts[iii]], self.Tz[1:])) + # if ((not self.c['SEB']) or (not self.c['manualT'])): + # self.Tz = np.concatenate(([self.Ts[iii]], self.Tz[1:])) + # self.Tz[0] = self.Ts[iii] + d_zbot = (self.z[-1]-zz1[-1]) + self.melt_sum = self.melt_sum + d_zbot self.w_firn = (znew - self.z_old) / self.dt[iii] # advection rate of the firn, m/s @@ -1226,6 +1296,7 @@ def time_evolve(self): ### NOTE: sigma = bdot_mean*GRAVITY*age/S_PER_YEAR*917.0) (or, sigma = bdot*g*tau, steady state conversion.) ### All temperature work happens at the end of the time loop + Ts_old = self.Tz[0] if self.c['manualT']: # manual temperature measurements are fed into CFM tif = interpolate.interp1d(self.manualT_dep, self.manualT_temp[:,iii],kind='cubic') @@ -1243,7 +1314,7 @@ def time_evolve(self): self.Tz, self.T10m = heatDiff(self,iii) dml_sum = 0 - elif np.any(self.LWC>0.): #VV enthalpy diffusion if water in column + elif np.any(self.LWC>0.): # enthalpy diffusion if water in column LWC0e = sum(self.LWC) tot_heat_pre = np.sum(CP_I_kJ*self.mass*self.Tz + T_MELT*CP_W/1000*self.LWC*RHO_W_KGM + LF_I_kJ*self.LWC*RHO_W_KGM) if "LWC_heat" not in self.c: @@ -1260,20 +1331,12 @@ def time_evolve(self): tot_heat_post = np.sum(CP_I_kJ*self.mass*self.Tz + T_MELT*CP_W/1000*self.LWC*RHO_W_KGM + LF_I_kJ*self.LWC*RHO_W_KGM) - ### debugging: check if total energy is the same pre and post diffusion solver. - # if (np.abs(tot_heat_post-tot_heat_pre)/tot_heat_pre)>1e-3: - # print(f'change in enthalpy at iteration {iii}!') - # print('pre:', tot_heat_pre) - # print('post:', tot_heat_post) - # ediff = (tot_heat_post-tot_heat_pre) - # print('difference (kJ):', (tot_heat_post-tot_heat_pre)) - # print('difference %:', ediff/tot_heat_pre) - self.refreeze += LWC0e-sum(self.LWC) self.T50 = np.mean(self.Tz[self.z<50]) # Temperature at 50 self.rho[self.rho>RHO_I] = RHO_I + ######### ############################################################# @@ -1286,7 +1349,11 @@ def time_evolve(self): if 'viscosity' in self.output_list: self.viscosity = RD['viscosity'] - self.climate = np.array([self.bdot[iii],self.Ts[iii]]) + if not self.c['SEB']: + self.climate = np.array([self.bdot[iii],self.Ts[iii]]) + else: #if SEB true + SMBiii = self.bdot[iii] + self.sublim[iii] - self.snowmelt[iii] #sublim negative means mass loss, snowmelt positive is amount lost + self.climate = np.array([SMBiii,self.Ts[iii]]) bcoAgeMart, bcoDepMart, bcoAge830, bcoDep830, LIZAgeMart, LIZDepMart, bcoAge815, bcoDep815 = self.update_BCO(iii) @@ -1355,11 +1422,13 @@ def time_evolve(self): if self.MELT: print(f'Totals (m w.e.)\n' - f'Melt+Rain: {sum(self.snowmeltSec+self.rainSec)*S_PER_YEAR*RHO_I_MGM}\n' - f'meltvol: {sum(meltvol2check)}\n' + f'Melt+Rain: {sum(self.snowmeltSec + self.rainSec)*S_PER_YEAR*RHO_I_MGM}\n' + f'meltvol: {sum(meltvol2check)}\n' #m w.e. f'Refreezing: {sum(refreezing2check)}\n' f'Runoff: {sum(runoff2check)}\n' + f'mismatch: {self.mismatch}\n' f'LWC (current): {sum(self.LWC)}\n' + f'LWC (init): {self.LWC_init}\n' f'Refrz + Rnff +LWC: {sum(runoff2check)+sum(refreezing2check)+sum(self.LWC)}\n' f'DML: {sum(dml2check)}') write_nospin_hdf5(self,self.MOutputs.Mout_dict,self.forcing_dict) diff --git a/CFM_main/firn_density_spin.py b/CFM_main/firn_density_spin.py index 4baf267..0de586e 100755 --- a/CFM_main/firn_density_spin.py +++ b/CFM_main/firn_density_spin.py @@ -135,7 +135,10 @@ def __init__(self, config, climateTS = None): ### temperature ### if climateTS != None: input_temp = climateTS['TSKIN'] - input_bdot = climateTS['BDOT'] + try: + input_bdot = climateTS['BDOT'] + climateTS['SUBLIM'] + except: + input_bdot = climateTS['BDOT'] input_year_temp = input_year_bdot = climateTS['time'] else: @@ -164,6 +167,12 @@ def __init__(self, config, climateTS = None): except: self.bdot0 = np.mean(input_bdot) + # try: + # if self.c['MELT']: + # self.bdot0 = self.bdot0 - np.mean(climateTS['SMELT']) + # except: + # pass + if 'manual_climate' in self.c: pass else: @@ -210,10 +219,19 @@ def __init__(self, config, climateTS = None): THL = self.temp0 AHL = self.bdot0 + try: #VV use Reeh corrected T if self.c['ReehCorrectedT'] and self.c['MELT']: if climateTS != None: - input_snowmelt = climateTS['SMELT'] + if not self.c['SEB']: + input_snowmelt = climateTS['SMELT'] + else: # estimate melt using simple PDD + PDDs = climateTS['T2m'] - T_MELT + PDDs[PDDs<0] = 0 + fPDD = 5 # positive degree day factor, mm/C/day + days_per_timestep = 365.25/self.c['stpsPerYear'] + input_snowmelt_mm = PDDs * fPDD * days_per_timestep # mm w.e. per time step + input_snowmelt = input_snowmelt_mm / 1000 / 0.917 * self.c['stpsPerYear'] # m i.e./year else: input_snowmelt, input_year_snowmelt, input_snowmelt_full, input_year_snowmelt_full = read_input(os.path.join(self.c['InputFileFolder'],self.c['InputFileNamemelt'])) #VV meanmelt = np.mean(input_snowmelt) # mean melt per year [mIE/yr] (units are specified in Reeh 2008) @@ -399,6 +417,9 @@ def __init__(self, config, climateTS = None): self.c['no_densification']=False ####################### + if self.c['ReehCorrectedT']: + self.rho = 900 * np.ones_like(self.rho) # use to just set a solid ice column to initialize + ############################ ##### END INIT ############# ############################ @@ -444,7 +465,7 @@ def time_evolve(self): 'LWC': self.LWC, 'MELT': self.MELT, 'FirnAir': False, - 'bdot_av': self.bdot_av + 'bdot_av': self.bdot_av, } if self.c['physRho']=='MaxSP': diff --git a/CFM_main/main.py b/CFM_main/main.py index 7cc26dd..5d46833 100755 --- a/CFM_main/main.py +++ b/CFM_main/main.py @@ -6,6 +6,11 @@ for the spin up and main runs, and it runs the model. This can be bypassed if you write your own script to call firn_density_spin and firn_density_nospin.py +Note that this file is becoming deprecated (though still works, and is useful for +general runs). I have started using all-in-one scripts that set up model forcings, +configure model run details, and then launch the run. Please let me know if you +need a script like this, and I'll be happy to send you one. + Copyright © 2021 C. Max Stevens Distributed under terms of the MIT license. @@ -22,7 +27,7 @@ __author__ = "C. Max Stevens, Vincent Verjans, Brita Horlings, Annika Horlings, Jessica Lundin" __license__ = "MIT" -__version__ = "2.0.0" +__version__ = "2.2.0" __maintainer__ = "Max Stevens" __email__ = "maxstev@umd.edu" __status__ = "Production" diff --git a/CFM_main/melt.py b/CFM_main/melt.py index 0abdaba..634bac1 100755 --- a/CFM_main/melt.py +++ b/CFM_main/melt.py @@ -88,18 +88,22 @@ def bucket(self,iii): n_melted = ind1+1 # number of nodes melted ### Partially melted node properties ### - pm_mass = self.mass_sum[ind1]-melt_mass #remaining mass + pm_mass = self.mass_sum[ind1] - melt_mass #remaining mass pm_dz = pm_mass/self.rho[ind1] #remaining thickness pm_rho = self.rho[ind1] #density of the pm node pm_lwc = self.LWC[ind1]/self.dz[ind1]*pm_dz #LWC of the pm node pm_Tz = T_MELT dzo = self.dz.copy() + if ind1>0: - dh_melt = -1 * (np.sum(dzo[0:ind1]) + (dzo[ind1]-pm_dz)) + dh_melt = -1 * (np.sum(dzo[0:ind1]) + (dzo[ind1]-pm_dz)) #thickness of melted nodes else: dh_melt = -1 * (dzo[ind1]-pm_dz) + avg_dh_melted = -1 * dh_melt/n_melted # average thickness of melted nodes + + ### Liquid water input at the surface ### liq_in_mass = max(melt_mass + (np.sum(self.LWC[0:ind1+1]) - pm_lwc) * RHO_W_KGM, 0) #avoid negative lwcinput due to numerical round-off errors liq_in_vol = liq_in_mass/RHO_W_KGM @@ -112,28 +116,43 @@ def bucket(self,iii): liqmcinit = pm_lwc+sum(self.LWC[ind1+1:])+liq_in_vol #mass conservation checks ### Regridding ### - if ind1>0: - self.rho = np.concatenate((self.rho[ind1:-1],self.rho[-1]*np.ones(n_melted))) - # self.Tz = np.concatenate((self.Tz[ind1:-1],self.Tz[-1]*np.ones(n_melted))) - self.r2 = np.concatenate((self.r2[ind1:-1],self.r2[-1]*np.ones(n_melted))) - self.bdot_mean = np.concatenate((self.bdot_mean[ind1:-1],self.bdot_mean[-1]*np.ones(n_melted))) - self.age = np.concatenate((self.age[ind1:-1],self.age[-1]*np.ones(n_melted))) - self.Dcon = np.concatenate((self.Dcon[ind1:-1],self.Dcon[-1]*np.ones(n_melted))) - self.dzn = np.concatenate((np.zeros(n_melted),self.dz[1:])) - self.dzn = self.dzn[0:self.compboxes] + if melt_mass>0: + + if ind1>0: + self.rho = np.concatenate((self.rho[ind1:-1],self.rho[-1]*np.ones(n_melted))) + # self.Tz = np.concatenate((self.Tz[ind1:-1],self.Tz[-1]*np.ones(n_melted))) + self.r2 = np.concatenate((self.r2[ind1:-1],self.r2[-1]*np.ones(n_melted))) + self.bdot_mean = np.concatenate((self.bdot_mean[ind1:-1],self.bdot_mean[-1]*np.ones(n_melted))) + self.age = np.concatenate((self.age[ind1:-1],self.age[-1]*np.ones(n_melted))) + self.Dcon = np.concatenate((self.Dcon[ind1:-1],self.Dcon[-1]*np.ones(n_melted))) + self.dzn = np.concatenate((np.zeros(n_melted),self.dz[1:])) + self.dzn = self.dzn[0:self.compboxes] + else: + self.dzn = self.dz[0:self.compboxes] #VV avoids bug due to undefined self.dzn + + self.LWC = np.concatenate(([pm_lwc],self.LWC[ind1+1:-1],self.LWC[-1]*np.ones(n_melted))) + + keep_firnthickness = self.c['keep_firnthickness'] + + if keep_firnthickness: + nb_th = np.maximum(avg_dh_melted,self.dz[-1]) + self.dz = np.concatenate(([pm_dz],self.dz[ind1+1:-1],nb_th*np.ones(n_melted))) + zbot_old = self.z[-1] + else: + self.dz = np.concatenate(([pm_dz],self.dz[ind1+1:-1],self.dz[-1]*np.ones(n_melted))) + + self.Tz = np.concatenate(([pm_Tz],self.Tz[ind1+1:-1],self.Tz[-1]*np.ones(n_melted))) # PM layer should have temp=T_MELT + self.z = self.dz.cumsum(axis=0) + self.z = np.concatenate(([0],self.z[:-1])) + + self.mass = self.rho*self.dz + if self.doublegrid: # if we have doublegrid: need to adjust gridtrack + meltgridtrack = np.concatenate((self.gridtrack[ind1:-1],self.gridtrack[-1]*np.ones(n_melted))) + elif self.doublegrid==False: + meltgridtrack = np.zeros(nnd) #just return a zero array else: - self.dzn = self.dz[0:self.compboxes] #VV avoids bug due to undefined self.dzn - - self.LWC = np.concatenate(([pm_lwc],self.LWC[ind1+1:-1],self.LWC[-1]*np.ones(n_melted))) - self.dz = np.concatenate(([pm_dz],self.dz[ind1+1:-1],self.dz[-1]*np.ones(n_melted))) - self.Tz = np.concatenate(([pm_Tz],self.Tz[ind1+1:-1],self.Tz[-1]*np.ones(n_melted))) # PM layer should have temp=T_MELT - self.z = self.dz.cumsum(axis=0) - self.z = np.concatenate(([0],self.z[:-1])) - self.mass = self.rho*self.dz - if self.doublegrid: # if we have doublegrid: need to adjust gridtrack - meltgridtrack = np.concatenate((self.gridtrack[ind1:-1],self.gridtrack[-1]*np.ones(n_melted))) - elif self.doublegrid==False: - meltgridtrack = np.zeros(nnd) #just return a zero array + meltgridtrack = self.gridtrack + self.dzn = self.dz[0:self.compboxes] ### end regridding ### ### Calculate excessive LWC (above irreducible holding capacity) ### @@ -528,6 +547,11 @@ def bucket(self,iii): self.rho[self.rho>RHO_I] = RHO_I + ### Mass conservation check 2 ### + liqmcfinal = sum(self.LWC) + refrozentot + runofftot + if abs(liqmcfinal - liqmcinit) > 1e-3: + print(f'Mass conservation error (2) (melt.py) at step {iii}\n Init: {liqmcinit} m\n Final: {liqmcfinal} m') + total_liquid_mass_end = np.sum(self.LWC*RHO_W_KGM) mass_runoff = runofftot*RHO_W_KGM mass_refreeze = refrozentot*RHO_W_KGM @@ -535,7 +559,8 @@ def bucket(self,iii): tot_mass_end = total_liquid_mass_end+mass_refreeze+mass_runoff mass_diff = tot_mass_end-total_liquid_mass_start - return self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, meltgridtrack, refrozentot, runofftot, dh_melt + return self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, \ + self.dzn, self.LWC, meltgridtrack, refrozentot, runofftot, dh_melt ############# @@ -948,3 +973,5 @@ def darcyscheme(self,iii): + + diff --git a/CFM_main/plotter.py b/CFM_main/plotter.py deleted file mode 100644 index 8018ecc..0000000 --- a/CFM_main/plotter.py +++ /dev/null @@ -1,47 +0,0 @@ -#!/usr/bin/env python -''' -Example script to load and plot results from the CFM. -''' - -# Copyright (C) 2019 C. Max Stevens -# Distributed under terms of the MIT license. - -import matplotlib.pyplot as plt -import h5py as h5 -import os -import sys -import numpy as np -import scipy as sp -import pickle -import seaborn as sns -sns.set() - -def plotter(rfolder,saver): - ''' - Function that plots. - ''' - - rfile = 'CFMresults.hdf5' - fn = os.path.join(rfolder,rfile) - f = h5.File(fn,'r') - - timesteps = f['depth'][1:,0] - stps = len(timesteps) - depth = f['depth'][1:,1:] - density = f['density'][1:,1:] - temperature = f['temperature'][1:,1:] - dip_all = f['DIP'][:,:] - f.close() - - f1,a1 = plt.subplots() - a1.plot(density[-1,:],depth[-1,:]) - a1.invert_yaxis() - a1.grid(True) - if saver: - f1.savefig('Example_DepthDensity.png') - -if __name__ == '__main__': - - saver = True - rfolder = '/PATH/TO/RESULTS/FOLDER' # alter this to point to the results folder. - diff --git a/CFM_main/solver.py b/CFM_main/solver.py index 318c8ea..f7ba6ff 100755 --- a/CFM_main/solver.py +++ b/CFM_main/solver.py @@ -60,9 +60,10 @@ def transient_solve_TR(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, ''' phi_t = phi_0 + phi_t_old = phi_t.copy() for i_time in range(nt): - + dZ = np.diff(z_edges) #width of nodes deltaZ_u = np.diff(Z_P) @@ -177,7 +178,9 @@ def transient_solve_TR(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, bc_u = np.concatenate(([ bc_u_0], [bc_type_u])) bc_d_0 = 0 - bc_type_d = 2 + bc_type_d = 2 # 2 is gradient + # bc_d_0 = 273.149 + # bc_type_d = 1 bc_d = np.concatenate(([ bc_d_0 ], [ bc_type_d ])) ######################################### @@ -200,6 +203,7 @@ def transient_solve_TR(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, b[-1] = bc_d[0] phi_t = solver(a_U, a_D, a_P, b) + a_P = a_U + a_D + a_P_0 if airdict!=None: @@ -262,7 +266,9 @@ def transient_solve_EN(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, itercheck = 0.9 count = 0 + ### Big H stands for latent enthalpy, little h is sensible enthalpy + ### H_tot is the sum of latent and sensible enthalpy H_lat = H_L_liq*g_liq # Latent enthalpy for each layer H_lat_old = H_lat.copy() @@ -310,11 +316,10 @@ def transient_solve_EN(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, a_U = D_u #* dt # [W/m2/K] a_D = D_d #* dt # [W/m2/K] - # c_vol1 = RHO_I*CP_I*g_sol_old - c_vol1 = RHO_I * CP_I - - a_P_0 = c_vol1 * dZ / dt # [W/m2/K] (new) Patankar eq. 4.41c, this is b_p in Voller (1990; Eq. 30) + c_vol1 = RHO_I * CP_I # gets multiplied by g_sol below + a_P_0 = c_vol1 * dZ / dt # [W/m2/K] (new) Patankar eq. 4.41c, this is b_p in Voller (1990; Eq. 30) + if update_gsol: a_P = a_U + a_D + a_P_0 * g_sol_iter - S_P * dZ #* dt # check the multiply on the S_P else: @@ -322,6 +327,7 @@ def transient_solve_EN(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, b_0 = S_C * dZ/dt #* dt # [W/m2] b = b_0 + a_P_0 * g_sol_old * phi_t_old # By this phi_t_old has to be in K + ############### ### Boundary conditions: @@ -330,7 +336,6 @@ def transient_solve_EN(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, # bc_u_0 = phi_s bc_u_0 = phi_t_old[0] # bc_u_0 = phi_iter[0] - bc_type_u = 1 bc_u = np.concatenate(([ bc_u_0], [bc_type_u])) @@ -347,7 +352,6 @@ def transient_solve_EN(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, #Down boundary a_P[-1] = 1 a_D[-1] = 0 - if bc_type_d==2: # Gradient (flux) a_U[-1] = 1 b[-1] = deltaZ_u[-1] * bc_d[0] @@ -359,56 +363,23 @@ def transient_solve_EN(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, phi_t = solver(a_U, a_D, a_P, b) #sensible enthalpy. 0 for layers at freezing (have LWC), negative for dry layers ##### - ### Adjust liquid fraction and temperture field based on solution - ### Calculations are (partially) based on the fact that the freezing temp is 0, - ### which means that if H_tot<0 there is no liquid and you can calculate temperature, and if H_tot>0 there is liquid and T is 0. - # if update_gsol: - # h_updated = phi_t * CP_I * RHO_I * g_sol_iter # updated sensible enthalpy after solver - # else: - - ### from the Neumann notebook - # h_updated = phi_t * CP_I * RHO_I * g_sol_old # updated sensible enthalpy after solver - - # delta_h = (h_updated - h_old) * 0.6 # overshoot correction. delta_h should always be negative for layers with LWC - # h_updated_trun = h_old + delta_h # sensible enthalpy after overshoot - - # cond0 = ((delta_h>0) & (LWC_old>0)) # Layers where there was sensible enthalpy increased and there is liquid water - # delta_h[cond0] = 0 # sensible enthalpy should not increase if LWC present (should either stay at 0C or cool down) - - # ndh = -1*delta_h #negative delta_h (which makes it positive), makes corrections below easy - - # ### everything refreezes if the calculated change in enthalpy is greater than the latent enthalpy - # cond1 = ((ndh>=H_lat_old) & (g_liq_old>0)) #layers where energy change is larger than the needed to refreeze, and where there is water - # H_tot[cond1] = (delta_h[cond1] + H_tot_old[cond1]) # total enthalpy in those layers is change+total, should be net negative - # g_liq[cond1] = 0 #no liquid left - # g_sol[cond1] = (mass_tot[cond1]/RHO_I)/dz[cond1] - # phi_t[cond1] = H_tot[cond1]/(CP_I * RHO_I * g_sol[cond1]) - - # ### partial refreezing if the delta_h is less than the latent enthalpy - # cond2 = ((ndh0)) - # g_liq[cond2] = (H_lat_old[cond2] + delta_h[cond2])/H_L_liq #remaining liquid - # dgl = g_liq - g_liq_old - # H_tot[cond2] = H_L_liq * g_liq[cond2] - # phi_t[cond2] = 0 - # g_sol[cond2] = g_sol[cond2] + -1 * dgl[cond2]/0.917 # Update solid fraction - - # ### Make sure that there is no liquid in layers that did not have liquid at start - # cond3 = (g_liq_old<=0) - # g_liq[cond3] = 0 - # H_tot[cond3] = h_updated[cond3] - # g_sol[cond3] = (mass_tot[cond3]/RHO_I)/dz[cond3] - - # g_liq[g_liq<0]=0 - # g_liq[g_liq>1]=1 - ############### - - - ### Older way - h_updated = phi_t * CP_I * RHO_I * g_sol_old # updated sensible enthalpy after solver + ''' + #### + The crux is to adjust liquid fraction and temperture field based on solution + Calculations are (partially) based on the fact that the freezing temp is 0, + which means that if H_tot<0 there is no liquid and you can calculate temperature, and if H_tot>0 there is liquid and T is 0. + + Note previous ways of solving in dev branch and previous releases. + ### + ''' - delta_h = (h_updated - h_old) * 0.6 # overshoot correction. delta_h should always be negative for layers with LWC - h_updated_trun = h_old + delta_h # sensible enthalpy after overshoot + ### The best way to solve: calculate new g_liq, and apply overshoot correction on g_liq + ### Break loop if iteration (prior to overshoot) is the same as previous solution + h_updated = phi_t * CP_I * RHO_I * g_sol_old # updated sensible enthalpy after solver. g_sol is volume_solid/dz + delta_h = h_updated - h_old # change in sensible enthalpy, relative to initial (not iteration) + + ### Figure out what delta_h and g_liq should be based on different conditions cond0 = ((delta_h>0) & (LWC_old>0)) # Layers where there was sensible enthalpy increased and there is liquid water delta_h[cond0] = 0 # sensible enthalpy should not increase if LWC present (should either stay at 0C or cool down) @@ -427,29 +398,30 @@ def transient_solve_EN(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, ### Make sure that there is no liquid in layers that did not have liquid at start cond3 = (g_liq_old<=0) g_liq[cond3] = 0 - # H_tot[cond3] = h_updated_trun[cond3] - H_tot[cond3] = h_updated_trun[cond3] + H_tot[cond3] = h_updated[cond3] + + ### if this iteration gave the same solution as the last iteration, break + iterdiff = (np.sum(g_liq_iter) - np.sum(g_liq)) # Deprecated? Used to use to calulate time to break loop + if ((np.allclose(phi_iter,phi_t,rtol=1e-4,atol=1e-3))): + break + elif ((np.allclose(g_liq_iter,g_liq,rtol=1e-4,atol=1e-4))): + break + ### otherwise apply overshoot correction to g_liq and iterate again. + delta_g_liq = g_liq - g_liq_iter # change in liquid fraction. Should always be negative. + g_liq = g_liq_iter + 0.6*delta_g_liq g_liq[g_liq<0]=0 g_liq[g_liq>1]=1 ################ - delta_g_liq = g_liq - g_liq_old # change in liquid fraction. Should always be negative. - ###### - g_sol = g_sol + -1 * delta_g_liq/0.917 # Update solid fraction ## Now update temperatures after liquid corrections phi_t[H_tot>=0] = 0 # H_tot>0 means liquid present, T=0 phi_t[H_tot<0] = H_tot[H_tot<0] / (CP_I * RHO_I * g_sol_old[H_tot<0]) # H_tot<0 means no liquid; all enthalpy is sensible - - iterdiff = (np.sum(g_liq_iter) - np.sum(g_liq)) # Deprecated? Used to use to calulate time to break loop - - if ((np.allclose(phi_iter,phi_t,rtol=1e-4,atol=1e-3))): - break - - elif ((np.allclose(g_liq_iter,g_liq,rtol=1e-4,atol=1e-4))): - break + + phi_t[g_liq>0] = 0 + ############# count += 1 @@ -666,6 +638,7 @@ def apparent_heat(z_edges, Z_P, nt, dt, Gamma_P, phi_0, nz_P, nz_fv, phi_s, mix_ ### end apparent_heat ######## ################################### + ''' Functions below are for firn air Works, but consider to be in beta diff --git a/CFM_main/sublim.py b/CFM_main/sublim.py index 153ee46..d78005b 100755 --- a/CFM_main/sublim.py +++ b/CFM_main/sublim.py @@ -23,17 +23,8 @@ def sublim(self,iii): sublim_mass = sublim_volume_WE * 1000. # [kg] ind1a = np.where((np.cumsum(self.mass)+np.cumsum(1000*self.LWC)) <= sublim_mass)[0] # indices of boxes that will be sublimated away num_boxes_sublim = len(ind1a)+1 # number of boxes that sublimate away, include the box that is partially sublimated - try: - ind1 = np.where((np.cumsum(self.mass)+np.cumsum(1000*self.LWC)) > sublim_mass)[0][0] # index which will become the new surface - except: - print('error!') - print(iii) - print(self.modeltime[iii]) - print('sublim_mass',sublim_mass) - print('rho',self.rho[0:20]) - print('tz',self.Tz[0:20]) - print('sublimsec',(self.sublimSec[iii] * S_PER_YEAR)) - sys.exit() + ind1 = np.where((np.cumsum(self.mass)+np.cumsum(1000*self.LWC)) > sublim_mass)[0][0] # index which will become the new surface + # ps is the partial sublimation (the model volume that has a portion sublimated away) dzo = self.dz.copy() @@ -45,7 +36,7 @@ def sublim(self,iii): ps_sublimmass = ps_sublim - ps_sublimlwc # rest of sublimated mass will be ice mass ps_lwc = np.maximum(0,(self.LWC[ind1] - ps_sublimlwc/1000)) # new surface layer gets its lwc reduced but not below 0 ps_plwc = np.maximum(self.PLWC_mem[ind1] - ps_sublimlwc/1000, 0.) # assumes plwc is sublimated before mlwc - ps_mass = np.maximum(1.0e-6,(self.mass[ind1] - ps_sublimmass)) # finally ice is sublimated if there is still mass to sublimate (<-> if ps_sublim<1000*self.LWC[ind1]), avoid rounding errors to cause negative mass + ps_mass = float(np.maximum(1.0e-6,(self.mass[ind1] - ps_sublimmass))) # finally ice is sublimated if there is still mass to sublimate (<-> if ps_sublim<1000*self.LWC[ind1]), avoid rounding errors to cause negative mass ps_dz = ps_mass / self.rho[ind1] # remaining thickness [m] if ind1>0: dh_sub = -1 * (np.sum(dzo[0:ind1]) + (dzo[ind1]-ps_dz)) @@ -65,8 +56,24 @@ def sublim(self,iii): self.age = np.concatenate((self.age[ind1:-1] , self.age[-1]*np.ones(num_boxes_sublim))) #+ self.dt[iii] # age of each layer increases of dt # self.dz = np.concatenate((self.dz[ind1:-1] , self.dz[-1]/divider*np.ones(num_boxes_sublim))) # this splits the last box into many. - self.dz = np.concatenate((self.dz[ind1:-1] , self.dz[-1]*np.ones(num_boxes_sublim))) # this adds new boxes at the bottom. - self.dz[0] = ps_dz #VV dz calculated for the partially sublimated layer + + keep_firnthickness = self.c['keep_firnthickness'] + + if keep_firnthickness: + avg_dh_sub = -1 * dh_sub/num_boxes_sublim # average thickness of melted nodes + + nb_th = np.maximum(avg_dh_sub,self.dz[-1]) + self.dz = np.concatenate(([ps_dz],self.dz[ind1+1:-1],nb_th*np.ones(num_boxes_sublim))) + zbot_old = self.z[-1] + else: + # pass + self.dz = np.concatenate(([ps_dz],self.dz[ind1+1:-1],self.dz[-1]*np.ones(num_boxes_sublim))) + + + + # self.dz = np.concatenate((self.dz[ind1:-1] , self.dz[-1]*np.ones(num_boxes_sublim))) # this adds new boxes at the bottom. + # self.dz[0] = ps_dz #VV dz calculated for the partially sublimated layer + self.Dcon = np.concatenate((self.Dcon[ind1:-1] , self.Dcon[-1]*np.ones(num_boxes_sublim))) self.dzn = np.concatenate((np.zeros(num_boxes_sublim), self.dz[1:])) #this is not quite right because is assumes compaction for the pm box is zero. self.dzn = self.dzn[0:self.compboxes] diff --git a/changelog.md b/changelog.md index b203142..a836112 100644 --- a/changelog.md +++ b/changelog.md @@ -3,7 +3,12 @@ All notable changes to the Community Firn Model should be documented in this fil TL;DR: Write down the changes that you made to the the model in this document and update the version number here and in main.py, then update main branch on github. -To run the update: +General update protocol: +First, get the staging branch to match main branch. +Then, while on staging branch checkout files from whatever branch you want (e.g., git checkout dev melt.py) +While on staging, test functionality. At minimum, python main.py example.json should work. +Then, switch to main branch and merge from staging. +Then: git commit -a -m "updating to vX.Y.Z. Details in changelog." git push git tag -a vX.Y.Z -m "CFM version vX.Y.Z" @@ -12,7 +17,7 @@ git push origin vX.Y.Z Then, on github do a release, which will trigger an updated DOI. ## Current Version -2.1.0 +2.2.0 ## Full Documentation @@ -33,6 +38,37 @@ https://communityfirnmodel.readthedocs.io/en/latest/ - Goujon physics work, but could possibly be implemented more elegantly (it would be nice to avoid globals) - Not exactly in progress, but at some point adding a log file that gets saved in the results folder would be a good idea. +## [2.2.0] 2023-06-26 +### Notes +- This version fixes an issue in the SEB module where new snow that was added was too warm. It also updates the enthalpy solving routine. There are numerous small changes and fixes, detailed below. + +### Fixed +- *SEB.py, firn_density_nospin.py* There was an bug when running using the SEB module that new snowfall was always at the melting temperature (which led to very warm firn). It now correctly sets the new snow temperature to be the minium of the freezing temperature and the air temperature + +- *SEB.py* Fixed an issue to set the albedo to some value if it is NaN in the input (MERRA-2 albedo is assinged as NaN when there is no SW flux). Even with no SW flux the NaNs caused an issue in the SEB calculation. + +- *RCMpkl_to_spin.py* Fixed an issue that the domain depth was always being set to be the depth of the 916 density horizon, rather than being configurable. + +- *melt.py* Fixed an issue for timesteps where there was rain input but no melting, in which the grid was still being altered (removing top layers, adding new layers to the bottom to keep number of layers constant). Now rain-only timesteps do not alter the grid structure (just adds mass). + +- *sublim.py* Fixed an issue where the sublimated mass was an numpy array with a single value rather than a float, which caused issues elsewhere when using that value. + +- *firn_density_nospin.py* Fixed an issue where surface temperature was concatenated (effectively advecting the temperature profile) rather than reset during timesteps without accumulation. + +### Changed +- *SEB.py* There is a new solver built into the code (FQS, or fast quartic solver) to find the surface temperature based on the previous surface temperature. The code (presently just in SEB.py, not as an option in the config file) uses a defined thickness (presently 1 cm) as the top layer from an SEB perspective; i.e., the energy is assumed to warm/cool that layer, and that layer's temperature is set to the new temperature at the end of the SEB routine (regardless of the resolution of the CFM grid). (Note that there are several other SEB solving schemes, but FQS is the most tested code. In theory they should all give the same value.) + +- *firn_density_spin.py* The firn column is initialized as a solid-ice column is ReehCorrectedT is true in the .json. The effect of this should be removed during spin up, but in areas around the ELA the column was sometimes disappearing or becoming too thin. + +- *melt.py, sublim.py* Previously, the CFM kept the number of grid layers constant when there was melt or sublimation by dividing the bottom layer into thinner layers. In places with a lot of melt this led to the column becoming very thin. Now, the model instead adds layers to the bottom that have the same properties (temperature, density) as the current bottom layer and the same thickness as was melted/sublimated away. Note that if the bottom of the domain is not close to the ice density this will affect the DIP calculation. This is toggled in the .json as 'keep_firnthickness'. 'keep_firnthickness' = False is the old behavior. + +- *example.json* Addition of 'keep_firnthickness' (see above). + +- *solver.py* Another update to the enthalpy solver. Previously, the overshoot factor applied to each iteration was applied to the new enthalpy solution, and now it just applied to the liquid portion. Additionally, the iteration's break statement is now prior to the overshoot, which (seems to) solves an issue where the solution was bouncing between values on each iteration instead of converging. + +### Removed +- *plotter.py* Plotter.py was written long ago and not updated. There are better ways of plotting the model outputs so I removed this file from the repository. Please let me know if you need a script or jupyter notebook to help with output processing/plotting and I can send you something. + ## [2.1.0] 2023-05-31 ### Notes - This version fixes an issue with how DIP was saved and updates the isotope diffusion module to allow inputs from a dictionary (parallel structure to loading temperature, accumulation rate, etc.)