Call for experiences: CICE4 restart conversion

[phil-blain]

Hi everyone,

We (ECCC) are moving forward with our transition to CICE6 and are having difficulty with starting CICE6 from CICE4 restarts.
At this point, we would just like to hear if others in the Consortium did that and what difficulty was encountered.

Our CICE4 restarts are in an in-house format (😢), so we can't directly use the included conversion program (https://github.com/CICE-Consortium/CICE/blob/main/configuration/tools/cice4_restart_conversion/convert_restarts.f90), but I had a student last term basically re-writing a Python program doing exactly the same logic to convert from our in-house format, used for CICE4, to CICE6 NetCDF.

Apart from the snow zapping code which is commented in the program:

CICE/configuration/tools/cice4_restart_conversion/convert_restarts.f90

Lines 225 to 237 in b946a95

	!!!!!!! modify as needed (begin)
	! if your restarts do not work, try uncommenting this section
	! ! snow temperature
	! zTsn = (Lfresh + trcrn(i,j,nt_qsno+k-1,n)/rhos)/cp_ice
	! ! zap entire snow volume if temperature is out of bounds
	! if (zTsn < Tmin .or. zTsn > Tmax) then
	! print*, 'zapping snow volume ', i,j,n,vsnon(i,j,n)
	! vsnon(i,j,n) = c0
	! do m = 1, nslyr
	! trcrn(i,j,nt_qsno+m-1,n) = c0
	! enddo
	! endif
	!!!!!!! modify as needed (begin)

which we indeed had to implement, do you remember having to manually adjust the converted restarts in anyway ?

/cc @JFLemieux73 @dupontf

[dabail10]

Hi Phillipe,

We have had similar experiences. I have tools in the NCAR Command Language, but one could imagine a similar tool in Fortran. They key issue is the energy versus enthalpy conversion. Snow enthalpy is very tricky. We ended up setting the snow enthalpy based on simple a simple function of the surface temperature. We had troubles just converting esnon to qsno. Probably one could just set a maximum value here. Also, salinity was added to the restart because of the mushy-layer thermodynamics (ktherm=2).

Here is some NCL code we used for this purpose. Note that we also went from 4 levels to 8 in the ice and 1 to 3 in the snow. If your levels do not change, you won't need the extra level code here.

qice = new((/nilyr,ncat,nj,ni/),double)
sice = new((/nilyr,ncat,nj,ni/),double)
qsno = new((/nslyr,ncat,nj,ni/),double)

printVarSummary(qice)

; Convert energy to enthalpy
do n=0,ncat-1
do k=0,nilyr1-1
   nt = n*nilyr1+k
   tmp = where(vicen(n,:,:).ne.0.,vicen(n,:,:)/int2dble(nilyr1),default_fillvalue("double"))
   tmp@_FillValue = default_fillvalue("double")
   if (phys .eq. "cice5") then
      qice(2*k,n,:,:) = eicen(nt,:,:)*0.5/tmp
      qice(2*k,n,:,:) = where(ismissing(qice(2*k,n,:,:)),0.,qice(2*k,n,:,:))
      qice(2*k+1,n,:,:) = eicen(nt,:,:)*0.5/tmp
      qice(2*k+1,n,:,:) = where(ismissing(qice(2*k+1,n,:,:)),0.,qice(2*k+1,n,:,:))
   end if
   if (phys .eq. "cice4") then
      qice(k,n,:,:) = eicen(nt,:,:)/tmp
      qice(k,n,:,:) = where(ismissing(qice(k,n,:,:)),0.,qice(k,n,:,:))
   end if
   print(n+1)
   print(k+1)
   tmp1 = ndtooned(qice(k,n,:,:))
   ii = minind(tmp1)
   tmp2 = ndtooned(eicen(nt,:,:))
   tmp3 = ndtooned(vicen(n,:,:))
   print(tmp1(ii))
   print(tmp2(ii))
   print(tmp3(ii))
end do
end do

rhos = 330.d0
cp_ice = 2.11727d3
Lfresh = 3.337d5

do n=0,ncat-1
do k=0,nslyr-1
   nt = n*nslyr1+k
   tmp = where(vsnon(n,:,:).ne.0.,vsnon(n,:,:)/int2dble(nslyr1),default_fillvalue("double"))
   tmp@_FillValue = default_fillvalue("double")
;  qsno(k,n,:,:) = esnon(nt,:,:)/tmp
   qsno(k,n,:,:) = -rhos*(Lfresh - cp_ice*Tsfcn(n,:,:))
   qsno(k,n,:,:) = where(ismissing(qsno(k,n,:,:)),0.,qsno(k,n,:,:))
end do
end do

; Salinity
saltmax = 3.2d0
nsal = 0.407d0
msal = 0.573d0
pi = atan(1.0d0)*4.0d0

salinz = new((/nilyr/),double)
do k=0,nilyr-1
   zn = (int2dble(k+1)-0.5d0)/int2dble(nilyr)
   salinz(k) = (saltmax/2.d0)*(1.d0-cos(pi*zn^(nsal/(msal+zn))))
   sice(k,:,:,:) = salinz(k)
end do

[daveh150]

Hi Philippe,

Snow enthalpy was the worst problem for us. I ended up adding this to convert_restarts to always recompute the qsno

           ! snow temperature
           zTsn = (Lfresh + trcrn(i,j,nt_qsno+k-1,n)/rhos)/cp_ice


           ! dah: added method to specify qsno based on Tsfc.
           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 ! the usual case for NRL
                 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
           endif ! zTsn out of bounds

Also, for completeness, since we used the CICE4 format, our typical methodology is:

1. Run the convert_restarts.f90 program (with above to re-compute qsno).
2. In CICE_InitMod.F90, make sure to uncomment the line to 'call restartfile_v4 (ice_ic)'. (Also need to add restartfile_v4 in the top include statement of CICE_InitMod.F90::ice_restart)
3. Must set runtype == 'initial' and specify the converted file in ice_ic for step 2 to work (that has caused me problems in the past)

Hope that helps. I have not tried to add the methodology of restartfile_v4 to a separte routine to write a NetCDF, that might be helpful as I think about it more.

Thanks,
Dave

[phil-blain]

Hi @dabail10 and @daveh150 , thanks a lot for your insights. I'll definitely keep those in mind if we hit more problems later on.

In the end, what we realized was that one of our in-house tools was writing the netCDF _FillValue attribute with value 1e30 for land points for all fields in the restart, and this lead to the model misbehaving(!). This was "fixed" by setting land points to zero instead, which is what the model does for its restarts.

From what @dupontf was able to deduce, it seems the advection code somehow uses those values on land, which is a little unsettling. I'll try to investigate that separately.

[dabail10]

The advection needs to run at points outside the ice pack. So, yes the restart should have zeroes for land. We had the problem where we had land block elimination that changed from the initial file to a new run and some land FillValues became 1.0e30. We do a check in the io_pio2 restart module for PIO_FILL_DOUBLE. However, this could be generalized to look for values greater than 1.0e29 or something.

[phil-blain]

Interesting, I was not aware of that . Note that we tried with upwind instead of remap, and this worked without zero-ing out land points.