ice:vertical thermo error after assimilation

[Marina]

Hi everyone,

I am having an ice:vertical thermo error when I run the ROMS/CICE model after updating the CICE initial file after assimilation. This is the message I get from CICE:

Starting thermo, zTsn > Tmax​

zTsn= 158.594491927825​

Tmax= 0​

Starting thermo, zTsn > Tmax​

zTsn= 158.594491927825​

Tmax= 0.000000000000000E+000​

istep1, my_task, i, j: 316801 34 7 2​

zqsn 0.000000000000000E+000 -110220000.000000 110220000.000000​

istep1, my_task, iblk = 316801 34 1​

category n = 1​

Global block: 35​

Global i and j: 6 919​

Lat, Lon: 87.5075809260643 64.0937710028015​

aice: 0.990060922517092​

n: 1 aicen: 2.887602290415366E-002​

I do a postprocessing before running the model with my new cice initial file provided by the assimilation. For instance, I set to zero several variables when there is no sea ice. Related to the previous error, I apply the following conditions:
1) Over sea ice, 'qsno001' maximum value is -1.2e8; over sea water this variable equals 0.
2) I have set 'Tsfcn' maximum value to -1.8 over sea ice, and over sea water the global maximum is 0. This variable can reach minimums of -20 over sea ice.
If I unstersdand correctly 'qsno001' corresponds to the snow enthalpy (maximum value -rhos*Lfresh = -110220000.0), and 'Tsfcn' is the temperature of ice/snow top surface (in category n).
I looked at the global i and j shown on the CICE message, but I didn't see any strange value on my variables at that location. I don't know what else to try...don't hesitate to share any ideas/possible solutions!
Thank you!

Best regards,

Marina

 

[dbailey]

I may have been adding some extra zeros, but you are right that it should be a maximum of -rhos*Lfresh (-330 * 3.34x10^5). Do you have missing values in your restart? All of the land values should be zero. Do you have a plus rhos*Lfresh here?

 

[david_hebert@nrlssc_navy_mil]

When I have run into this it is due to having zero for qsno, causing the large temperature. What I did was use a method in

configuration/tools/cice4_restart_conversion/convert_restarts.f90​

where around line 230, instead of removing snow recompute the qsno:

if (zTsn < Tmin .or. zTsn > Tmax) then
print*, 'recomputing qsno i,j,n,k', i,j,n,k
Ti = min(c0,trcrn(i,j,nt_Tsfc,n)) ! check temp not > 0
if (heat_capacity) then
trcrn(i,j,nt_qsno+k-1,n) = -rhos *(Lfresh - cp_ice*Ti)
else
trcrn(i,j,nt_qsno+k-1,n) = -rhos * Lfresh
endif

Hope that helps. David

 

[TillRasmussen]

I am not sure which version of CICE you are working with, however in v6.3 the error message originates from init_vertical_profile
For Tmax to be zero the enthalpy must be zero as written in
zqsn 0.000000000000000E+000 -110220000.000000 110220000.000000
(zqsn(k), -Lfresh*rhos,zqsn(k)+Lfresh*rhos)
There is snow in the cell otherwise Tmax would be set to a small number in the calculation.

 

[Marina]

Hi, thank you everyone for your input!!

I may have been adding some extra zeros, but you are right that it should be a maximum of -rhos*Lfresh (-330 * 3.34x10^5). Do you have missing values in your restart? All of the land values should be zero. Do you have a plus rhos*Lfresh here?

In my file, qsno001 maximum negative value is -128989327.94752154 (so lower than the limit = -1.1022x10^8). In sea water and land, its value is zero. If there is no snow over sea ice, I understand that qsno001 should be zero too? Is there a relation to keep between qsno001 and vsnon?

I am not sure which version of CICE you are working with, however in v6.3 the error message originates from init_vertical_profile
For Tmax to be zero the enthalpy must be zero as written in
zqsn 0.000000000000000E+000 -110220000.000000 110220000.000000
(zqsn(k), -Lfresh*rhos,zqsn(k)+Lfresh*rhos)
There is snow in the cell otherwise Tmax would be set to a small number in the calculation.

I am using CICE5.1.2 version, so probably the output error changes a little bit. So this means there is no snow in the cell where the computation fails?

When I have run into this it is due to having zero for qsno, causing the large temperature. What I did was use a method in

configuration/tools/cice4_restart_conversion/convert_restarts.f90​

where around line 230, instead of removing snow recompute the qsno:

if (zTsn < Tmin .or. zTsn > Tmax) then
print*, 'recomputing qsno i,j,n,k', i,j,n,k
Ti = min(c0,trcrn(i,j,nt_Tsfc,n)) ! check temp not > 0
if (heat_capacity) then
trcrn(i,j,nt_qsno+k-1,n) = -rhos *(Lfresh - cp_ice*Ti)
else
trcrn(i,j,nt_qsno+k-1,n) = -rhos * Lfresh
endif

Hope that helps. David

I cannot exactly find the same file convert_restarts.f90, I think my setup is different maybe due to a different cice version. Did you use this qsno computation when you run the model, or did you use it before to create your restart file?

Thank you!!

 

[TillRasmussen]

In CICE v5.1 the same output is available as in in 6.3.
The output is found in source/ice_therm_vertical.F90 line 1028.
Just above Tmax can only be zero if there is snow and the enthalpy (zqsn) is zero.
If zqsn is zero then zTsn = (Lfresh + zqsn(ij,k)/rhos)/cp_ice=Lfresh/cp_ice which result in ~158 degrees.

In your case you have some snow but enthalpy= zero

Instead of enthalpy = zero try enthalpy = -rhos * Lfresh as default (I think that this should result in Temperature=0
Is it binary or netcdf input you are using?

 

[Marina]

In CICE v5.1 the same output is available as in in 6.3.
The output is found in source/ice_therm_vertical.F90 line 1028.
Just above Tmax can only be zero if there is snow and the enthalpy (zqsn) is zero.
If zqsn is zero then zTsn = (Lfresh + zqsn(ij,k)/rhos)/cp_ice=Lfresh/cp_ice which result in ~158 degrees.

In your case you have some snow but enthalpy= zero

Instead of enthalpy = zero try enthalpy = -rhos * Lfresh as default (I think that this should result in Temperature=0
Is it binary or netcdf input you are using?

Thanks for that info! I will replace the null enthalpy values by -rhos * Lfresh and see what happens. I am using netcdf input

 

[Marina]

I put all qsno001 maximum values equal to -rhos * Lfresh, and remove all the zero values. However, the same error comes up.

Looking at source/ice_therm_vertical.F90, around line 983:

if (hslyr(ij) > hs_min/rnslyr .and. heat_capacity) then​

! zqsn < 0 ​

zqsn(ij,k) = trcrn(i,j,nt_qsno+k-1)​

Tmax = -zqsn(ij,k)*puny*rnslyr / &​

(rhos*cp_ice*vsnon(i,j))​

else​

zqsn (ij,k) = -rhos * Lfresh​

Tmax = puny​

endif​

​

!-----------------------------------------------------------------​

! Compute snow temperatures from enthalpies.​

! Note: zqsn <= -rhos*Lfresh, so zTsn <= 0.​

!-----------------------------------------------------------------​

zTsn(ij,k) = (Lfresh + zqsn(ij,k)/rhos)/cp_ice​

It seems that when it computes zqsn(ij,k) = trcrn(i,j,nt_qsno+k-1), qzsn becomes zero, and that is why Tmax = 0 and then zTsn = 158...but why does this happen? Let me know if you have any ideas! :) Thank you!

 

[dbailey]

I have a feeling that you have a land/ocean mask problem. Here is the difference between the post and nopost files. There is a big area near land where there are differences on the order of 1.3e8. So, it is likely you have valid ice points where the qsno001 is above -rhos*Lfresh.
Dave

 

[TillRasmussen]

trcrn is the global value that contains all kinds of tracers with dimension (i,j,number_tracers,category). In trcrn(i,j,nt_qsno+k-1) nt_qsno is the index for enthalpy snow. The line you refer to should copy the value of enthalpy for snow from the trcrn array to the local array zqsn.
Restart file should be read in on line 283 in source/ice_restart_driver.F90. I dont think that it should be updated before you initialize the vertical profile where it fails.

do k=1,nslyr
write(nchar,'(i3.3)') k
call read_restart_field(nu_restart,0,trcrn(:,:,nt_qsno+k-1,:,:),'ruf8', &
'qsno'//trim(nchar),ncat,diag,field_loc_center,field_type_scalar)
enddo
You could write ouput after this in order to see whether the values are as expected

I dont have the restart files that @dbailey refers to, thus it is a bit difficult to get closer to the answer but I hope that you will find the bug

 

[dbailey]

This assumes she is using the CICE Fortran code for the assimilation. :) She has a single snow layer, and so the field is just qsno001 in the restart file. This could also be a problem if you are changing a point with the assimilation that previously did not have sea ice.

 

[Marina]

Thank you both your answers. I have checked for all new ice points after assimilation and put qsno001 = -1.2e8 instead of 0 for those new ice points. However the error stays. I was wondering that I might need to udpate all variables in these new ice points, for instance, Tsfcn should not be equal to -1.8 in these new ice points, however this is the case, so I guess I should modify that, but what values should I put instead, -10 for example? It seems to be a little bit random the way of defining these values for new ice points after assimilation.....

 

[Marina]

This assumes she is using the CICE Fortran code for the assimilation. :) She has a single snow layer, and so the field is just qsno001 in the restart file. This could also be a problem if you are changing a point with the assimilation that previously did not have sea ice.

I am not using the CICE Fortran code for assimilation, I use another code separately from CICE

 

[TillRasmussen]

Could you write what is in your restart file both the initial and the one that is changed after the assimilation
Just point
Global i and j: 6 919
variables: ice concentration, ice thickness and snow thickness, lon and lat

 

[Marina]

Hi,

Thanks a lot for all your feedback!!

Finally I got to solve the problem. I did some changes in the ice_in file. I set tr_iage = .true. and trestore = 0 (instead of .false. and 90 which was before). And after that I could make it work and no more ice:vertical thermo error. After some tests, I figured out that trestore = 0 was the key parameter to change, as I see in the doc it is related to the boundary conditions, maybe that is why the global i and j points were mostly located near the boundary. If any of you have any info/experience on this parameter, I would be happy to hear it so I can understand better the role this parameter plays!

 

[Riley]

Hi, Marina. I met the similar problem. Variable name trestore is not found in the ice_in file. How to set the trestore = 0? Thanks!