Open guillaumevernieres opened 8 months ago
There are a few more cases for which the CICE model forecast becomes unstable.
Something needs to be done about snow enthalpy and small sea-ice concentration values. Solution from Philip Pegion:
import xarray as xr import numpy as np import sys import shutil infile=sys.argv[1] ds_in=xr.open_dataset(infile) var_list=['vicen','vsnon','Tsfcn','iage','alvl','vlvl', 'apnd', 'hpnd', 'ipnd',\ 'dhs', 'ffrac', 'sice001', 'qice001', 'sice002', 'qice002', 'sice003', 'qice003',\ 'sice004', 'qice004', 'sice005', 'qice005', 'sice006', 'qice006', 'sice007', 'qice007','qsno001','aicen'] for var in var_list: tmparr=np.where(ds_in['aicen'].values >0.15,ds_in[var],0.0) ds_in[var][:]=tmparr new_frac=tmparr.sum(axis=0) new_mask=np.where(new_frac > 0, 1.0,0.0) ds_in['iceumask'][:]=new_mask ds_in.to_netcdf('trimmed_'+infile) ds_in.close() shutil.move(infile,'old_'+infile) shutil.move('trimmed_'+infile,infile)
import xarray as xr import numpy as np import sys
infile=sys.argv[1] anl=xr.opendataset('anl'+infile) rhos = 330.0 cp_ice = 2106. c1 = 1.0 Lsub = 2.835e6 Lvap = 2.501e6 Lfresh=Lsub - Lvap rnslyr=1.0 puny=1.0E-012
A = c1 / (rhos cp_ice) B = Lfresh / cp_ice zTsn = A anl['qsno001'][:].values + B
Tmax = -anl['qsno001'][:].valuespunyrnslyr /(rhoscp_iceanl['vsnon'][:].values)
Qmax=rhoscp_ice(Tmax-Lfresh/cp_ice)
newq=np.where(zTsn <= Tmax,anl['qsno001'][:].values,Qmax) newf=np.where(anl['vicen'] > 0.00001,anl['aicen'][:].values,0.0)
newq2=np.where(anl['vsnon'][:]==0.0,anl['qsno001'][:].values,newq) anl['qsno001'][:]=newq2 anl['aicen'][:]=newf
anl.to_netcdf(infile)
The above features need to be added to the `soca2cice` var change.
Thanks for sharing @guillaumevernieres
Description
There are a few more cases for which the CICE model forecast becomes unstable.
Solution
Something needs to be done about snow enthalpy and small sea-ice concentration values. Solution from Philip Pegion:
infile=sys.argv[1] anl=xr.opendataset('anl'+infile) rhos = 330.0 cp_ice = 2106. c1 = 1.0 Lsub = 2.835e6 Lvap = 2.501e6 Lfresh=Lsub - Lvap rnslyr=1.0 puny=1.0E-012
icepack formulate for snow temperature
A = c1 / (rhos cp_ice) B = Lfresh / cp_ice zTsn = A anl['qsno001'][:].values + B
icepack formula for max snow tempature
Tmax = -anl['qsno001'][:].valuespunyrnslyr /(rhoscp_iceanl['vsnon'][:].values)
enthlap at max now tempetarure
Qmax=rhoscp_ice(Tmax-Lfresh/cp_ice)
fill in new enthalpy where snow temperature is too high
newq=np.where(zTsn <= Tmax,anl['qsno001'][:].values,Qmax) newf=np.where(anl['vicen'] > 0.00001,anl['aicen'][:].values,0.0)
fill in snow enthalpy (0) where there is no snow
newq2=np.where(anl['vsnon'][:]==0.0,anl['qsno001'][:].values,newq) anl['qsno001'][:]=newq2 anl['aicen'][:]=newf
write out file
anl.to_netcdf(infile)