'time_beg' sometimes not initialized before being written to 'time_bounds' NetCDF variable (NetCDF files not usable by Panoply) #551
Comments
|
Just to clarify, standalone CICE output is read fine by panopoly but the coupled CICE output is not? Is the time axis different in those two sets of files? Is anything else in the header different in the two files? If not, could there be some bad data in the file that is creating problems? Maybe try producing a CICE history file with one field in it and see if that still fails? |
Yes
The only difference I see is that calendar...: # Standalone
$ ncdump -h iceh_inst.2005-01-01-43200.nc
netcdf iceh_inst.2005-01-01-43200 {
dimensions:
d2 = 2 ;
ni = 100 ;
nj = 116 ;
nc = 5 ;
nkice = 7 ;
nksnow = 1 ;
nkbio = 9 ;
nkaer = 11 ;
time = UNLIMITED ; // (1 currently)
nvertices = 4 ;
nf = 1 ;
variables:
double time(time) ;
time:long_name = "model time" ;
time:units = "days since 2005-01-01 00:00:00" ;
time:calendar = "Gregorian" ;
time:bounds = "time_bounds" ;
float time_bounds(time, d2) ;
time_bounds:long_name = "boundaries for time-averaging interval" ;
time_bounds:units = "days since 2005-01-01 00:00:00" ;
# ...
# Coupled
$ ncdump -h iceh_inst.2001-10-05-07200.nc
netcdf iceh_inst.2001-10-05-07200 {
dimensions:
d2 = 2 ;
ni = 528 ;
nj = 735 ;
nc = 5 ;
nkice = 7 ;
nksnow = 1 ;
nkbio = 9 ;
nkaer = 11 ;
time = UNLIMITED ; // (1 currently)
nvertices = 4 ;
nf = 1 ;
variables:
double time(time) ;
time:long_name = "model time" ;
time:units = "days since 2001-01-01 00:00:00" ;
time:calendar = "NoLeap" ;
time:bounds = "time_bounds" ;
float time_bounds(time, d2) ;
time_bounds:long_name = "boundaries for time-averaging interval" ;
time_bounds:units = "days since 2001-01-01 00:00:00"
# ...
Yeah I'll try some things. I'll try running standalone with the same namelist as the coupled run. |
|
My first guess would be that netcdf is struggling with the "NoLeap" calendar name. I wonder if the caps are creating problems? I believe netcdf should understand the following, http://cfconventions.org/cf-conventions/cf-conventions.html#calendar. Trying changing NoLeap to noleap in the file and see if panopoly is happier. You should be able to do this with nco or a similar tool. |
|
I tried that with |
|
I think it's the I'm compiling with CICE/cicecore/cicedynB/analysis/ice_history.F90 Lines 1863 to 1869 in f5f487f somehow is not reached and so that NaN gets written to the history file... I will run the model in a debugger tomorrow to try to understand how that happens. |
|
The NaN appears also in a standalone run when compiling with |
phil-blain
changed the title
phil-blain
changed the title
|
OK, I understand more what happens. The code that initializes CICE/cicecore/cicedynB/analysis/ice_history.F90 Line 1817 in f5f487f but we use In any case I'll investigate more our in-house CICE4 code next week and compare. |
|
I'm bringing @dabail10 into the discussion since you worked on the multiple-stream history functionality (CICE-Consortium/CICE-svn-trunk@ab1c53e) for CICE 5. In that commit you removed that line in ice_init: @@ -327,7 +332,6 @@ subroutine input_data
ocn_data_format = 'bin'
#endif
- if (histfreq == '1') hist_avg = .false. ! potential conflict
if (days_per_year /= 365) shortwave = 'default' ! definite conflict
chartmp = advection(1:6)but there are still places in the code that make
Apart from these 3 places, I actually think that averaging over several time steps does work, because we have not made any changes to the averaging code in our version of CICE4. But I'd like others to confirm that my thinking is correct... |
|
I don't remember this very well, but it is possible I introduced a bug here. I guess I couldn't figure out how to average some streams but have others instantaneous. So, I guess I was thinking all or none. However, I would be happy to revisit this and see if there is a way to do this. I suppose making hist_avg into an array ... |
|
Well that would be useful but here I'm actually talking about a single stream, with CICE/cicecore/cicedynB/analysis/ice_history.F90 Lines 1817 to 1836 in f5f487f that explicitely reset the accumulation arrays to zero (if I read the code correctly...) |
|
Is that not the behaviour you want? If it is step output (1) then you have nothing to accumulate. |
|
Could it have been a confusion between histreq_n(is)=1 and histfreq(is)='1'? The second condition would make more sense indeed if it was [...] |
|
@dabail10 I think there is indeed some confusion. As Fred is hinting at, we are using |
|
Makes sense. There are a few things in the history output that should probably be cleaned up. |
The code at the beginning of subroutine 'accum_hist' is responsible for
resetting the accumulation arrays a2D, a3Dc, etc. to zero when writing
snapshots, but it also resets them if "histfreq(ns) == '1'".
This has been the case since ab1c53e (Many thanks to Dave Bailey (NCAR)
for getting this ball rolling., 2009-06-01) [1], which implemented the
multiple-stream history functionality in CICE 5.
This means that averaging over multiple time steps, i.e.
histfreq = '1'
histfreq_n = m ! m > 1
hist_avg = .true.
does not work: the accumulation arrays will always be reset to zero at
the beginning of 'accum_hist' and no accumulation will happen. This also
means that the variable 'time_beg' is not initialized in that
configuration, since this it is initialized in the 'else' branch of that
'if' statement.
Commit ab1c53e also *removed* the following line in ice_init::input_data:
if (histfreq == '1') hist_avg = .false. ! potential conflict
so it's unclear if that configuration was forgotten at the time.
Allow overaging over multiple time steps by only resetting the
accumulation arrays if 'hist_avg' is false.
As a side effect, this correctly initializes the 'time_beg' variable,
which avoids writing that variable uninitialized to the output files, which can
lead to analysis tools failing if compiling with '-init=snan,arrays' (or
similar [2].
While at it, remove a commented line in the 'else' branch since it's
implemented just below.
This commit does not modify neither
ice_history_shared::construct_filename nor the different
ice_history_write::ice_write_hist subroutines, that also assume that
'histfreq = 1' implies 'histfreq_n = 1'. They will be modified in
subsequent commits.
[1] CICE-Consortium/CICE-svn-trunk@ab1c53e?branch=ab1c53eaf0c4a0b0dbc65e582a5e1b35c2208a88&diff=unified#diff-aebd77751eba77ed6b6c194d4cd2e4b353258f7e903cc7877343911c9b25c283L1296-L1306
[2] Closes CICE-Consortium#551
The code at the beginning of subroutine 'accum_hist' is responsible for
resetting the accumulation arrays a2D, a3Dc, etc. to zero when writing
snapshots, but it also resets them if "histfreq(ns) == '1'".
This has been the case since ab1c53e (Many thanks to Dave Bailey (NCAR)
for getting this ball rolling., 2009-06-01) [1], which implemented the
multiple-stream history functionality in CICE 5.
This means that averaging over multiple time steps, i.e.
histfreq = '1'
histfreq_n = m ! m > 1
hist_avg = .true.
does not work: the accumulation arrays will always be reset to zero at
the beginning of 'accum_hist' and no accumulation will happen. This also
means that the variable 'time_beg' is not initialized in that
configuration, since this it is initialized in the 'else' branch of that
'if' statement.
Commit ab1c53e also *removed* the following line in ice_init::input_data:
if (histfreq == '1') hist_avg = .false. ! potential conflict
so it's unclear if that configuration was forgotten at the time.
Allow overaging over multiple time steps by only resetting the
accumulation arrays if 'hist_avg' is false.
As a side effect, this correctly initializes the 'time_beg' variable,
which avoids writing that variable uninitialized to the output files, which can
lead to analysis tools failing if compiling with '-init=snan,arrays' (or
similar [2].
While at it, remove a commented line in the 'else' branch since it's
implemented just below.
This commit does not modify neither
ice_history_shared::construct_filename nor the different
ice_history_write::ice_write_hist subroutines, that also assume that
'histfreq = 1' implies 'histfreq_n = 1'. They will be modified in
subsequent commits.
[1] CICE-Consortium/CICE-svn-trunk@ab1c53e?branch=ab1c53eaf0c4a0b0dbc65e582a5e1b35c2208a88&diff=unified#diff-aebd77751eba77ed6b6c194d4cd2e4b353258f7e903cc7877343911c9b25c283L1296-L1306
[2] Closes CICE-Consortium#551
Hi everyone!
I'm still working on coupling CICE with NEMO in our systems and the output files that are created can't be read by Panoply (4.10.12), I get:
I can open that file in Ncview just fine (2.1.8) , so I'm not sure where the problem is, but I thought I'd ask here since I'm sure some people here know more about NetCDF than me...
Here is the header of the file (
ncdump -c iceh_inst.2001-10-05-07200.nc):$ ncdump -c iceh_inst.2001-10-05-07200.nc netcdf iceh_inst.2001-10-05-07200 { dimensions: d2 = 2 ; ni = 528 ; nj = 735 ; nc = 5 ; nkice = 7 ; nksnow = 1 ; nkbio = 9 ; nkaer = 11 ; time = UNLIMITED ; // (1 currently) nvertices = 4 ; nf = 1 ; variables: double time(time) ; time:long_name = "model time" ; time:units = "days since 2001-01-01 00:00:00" ; time:calendar = "NoLeap" ; time:bounds = "time_bounds" ; float time_bounds(time, d2) ; time_bounds:long_name = "boundaries for time-averaging interval" ; time_bounds:units = "days since 2001-01-01 00:00:00" ; float TLON(nj, ni) ; TLON:long_name = "T grid center longitude" ; TLON:units = "degrees_east" ; TLON:missing_value = 1.e+30f ; TLON:_FillValue = 1.e+30f ; float TLAT(nj, ni) ; TLAT:long_name = "T grid center latitude" ; TLAT:units = "degrees_north" ; TLAT:missing_value = 1.e+30f ; TLAT:_FillValue = 1.e+30f ; float ULON(nj, ni) ; ULON:long_name = "U grid center longitude" ; ULON:units = "degrees_east" ; ULON:missing_value = 1.e+30f ; ULON:_FillValue = 1.e+30f ; float ULAT(nj, ni) ; ULAT:long_name = "U grid center latitude" ; ULAT:units = "degrees_north" ; ULAT:missing_value = 1.e+30f ; ULAT:_FillValue = 1.e+30f ; ULAT:comment = "Latitude of NE corner of T grid cell" ; float NCAT(nc) ; NCAT:long_name = "category maximum thickness" ; NCAT:units = "m" ; float VGRDa(nkaer) ; VGRDa:long_name = "vertical snow-ice-bio levels" ; VGRDa:units = "1" ; float tmask(nj, ni) ; tmask:long_name = "ocean grid mask" ; tmask:coordinates = "TLON TLAT" ; tmask:comment = "0 = land, 1 = ocean" ; tmask:missing_value = 1.e+30f ; tmask:_FillValue = 1.e+30f ; float blkmask(nj, ni) ; blkmask:long_name = "ice grid block mask" ; blkmask:coordinates = "TLON TLAT" ; blkmask:comment = "mytask + iblk/100" ; blkmask:missing_value = 1.e+30f ; blkmask:_FillValue = 1.e+30f ; float tarea(nj, ni) ; tarea:long_name = "area of T grid cells" ; tarea:units = "m^2" ; tarea:coordinates = "TLON TLAT" ; tarea:missing_value = 1.e+30f ; tarea:_FillValue = 1.e+30f ; float uarea(nj, ni) ; uarea:long_name = "area of U grid cells" ; uarea:units = "m^2" ; uarea:coordinates = "ULON ULAT" ; uarea:missing_value = 1.e+30f ; uarea:_FillValue = 1.e+30f ; float ANGLE(nj, ni) ; ANGLE:long_name = "angle grid makes with latitude line on U grid" ; ANGLE:units = "radians" ; ANGLE:coordinates = "ULON ULAT" ; ANGLE:missing_value = 1.e+30f ; ANGLE:_FillValue = 1.e+30f ; float ANGLET(nj, ni) ; ANGLET:long_name = "angle grid makes with latitude line on T grid" ; ANGLET:units = "radians" ; ANGLET:coordinates = "TLON TLAT" ; ANGLET:missing_value = 1.e+30f ; ANGLET:_FillValue = 1.e+30f ; float hi(time, nj, ni) ; hi:units = "m" ; hi:long_name = "grid cell mean ice thickness" ; hi:coordinates = "TLON TLAT time" ; hi:cell_measures = "area: tarea" ; hi:missing_value = 1.e+30f ; hi:_FillValue = 1.e+30f ; hi:cell_methods = "time: mean" ; hi:time_rep = "instantaneous" ; float hs(time, nj, ni) ; hs:units = "m" ; hs:long_name = "grid cell mean snow thickness" ; hs:coordinates = "TLON TLAT time" ; hs:cell_measures = "area: tarea" ; hs:missing_value = 1.e+30f ; hs:_FillValue = 1.e+30f ; hs:cell_methods = "time: mean" ; hs:time_rep = "instantaneous" ; float Tsfc(time, nj, ni) ; Tsfc:units = "C" ; Tsfc:long_name = "snow/ice surface temperature" ; Tsfc:coordinates = "TLON TLAT time" ; Tsfc:cell_measures = "area: tarea" ; Tsfc:missing_value = 1.e+30f ; Tsfc:_FillValue = 1.e+30f ; Tsfc:cell_methods = "time: mean" ; Tsfc:time_rep = "instantaneous" ; float aice(time, nj, ni) ; aice:units = "1" ; aice:long_name = "ice area (aggregate)" ; aice:coordinates = "TLON TLAT time" ; aice:cell_measures = "area: tarea" ; aice:missing_value = 1.e+30f ; aice:_FillValue = 1.e+30f ; aice:cell_methods = "time: mean" ; aice:time_rep = "instantaneous" ; float uvel(time, nj, ni) ; uvel:units = "m/s" ; uvel:long_name = "ice velocity (x)" ; uvel:coordinates = "ULON ULAT time" ; uvel:cell_measures = "area: uarea" ; uvel:missing_value = 1.e+30f ; uvel:_FillValue = 1.e+30f ; uvel:cell_methods = "time: mean" ; uvel:time_rep = "instantaneous" ; float vvel(time, nj, ni) ; vvel:units = "m/s" ; vvel:long_name = "ice velocity (y)" ; vvel:coordinates = "ULON ULAT time" ; vvel:cell_measures = "area: uarea" ; vvel:missing_value = 1.e+30f ; vvel:_FillValue = 1.e+30f ; vvel:cell_methods = "time: mean" ; vvel:time_rep = "instantaneous" ; float uocn(time, nj, ni) ; uocn:units = "m/s" ; uocn:long_name = "ocean current (x)" ; uocn:coordinates = "ULON ULAT time" ; uocn:cell_measures = "area: uarea" ; uocn:missing_value = 1.e+30f ; uocn:_FillValue = 1.e+30f ; uocn:cell_methods = "time: mean" ; uocn:time_rep = "instantaneous" ; float vocn(time, nj, ni) ; vocn:units = "m/s" ; vocn:long_name = "ocean current (y)" ; vocn:coordinates = "ULON ULAT time" ; vocn:cell_measures = "area: uarea" ; vocn:missing_value = 1.e+30f ; vocn:_FillValue = 1.e+30f ; vocn:cell_methods = "time: mean" ; vocn:time_rep = "instantaneous" ; float congel(time, nj, ni) ; congel:units = "cm/day" ; congel:long_name = "congelation ice growth" ; congel:coordinates = "TLON TLAT time" ; congel:cell_measures = "area: tarea" ; congel:missing_value = 1.e+30f ; congel:_FillValue = 1.e+30f ; congel:cell_methods = "time: mean" ; congel:time_rep = "instantaneous" ; float strairx(time, nj, ni) ; strairx:units = "N/m^2" ; strairx:long_name = "atm/ice stress (x)" ; strairx:coordinates = "ULON ULAT time" ; strairx:cell_measures = "area: uarea" ; strairx:missing_value = 1.e+30f ; strairx:_FillValue = 1.e+30f ; strairx:cell_methods = "time: mean" ; strairx:time_rep = "instantaneous" ; float strairy(time, nj, ni) ; strairy:units = "N/m^2" ; strairy:long_name = "atm/ice stress (y)" ; strairy:coordinates = "ULON ULAT time" ; strairy:cell_measures = "area: uarea" ; strairy:missing_value = 1.e+30f ; strairy:_FillValue = 1.e+30f ; strairy:cell_methods = "time: mean" ; strairy:time_rep = "instantaneous" ; float strocnx(time, nj, ni) ; strocnx:units = "N/m^2" ; strocnx:long_name = "ocean/ice stress (x)" ; strocnx:coordinates = "ULON ULAT time" ; strocnx:cell_measures = "area: uarea" ; strocnx:missing_value = 1.e+30f ; strocnx:_FillValue = 1.e+30f ; strocnx:cell_methods = "time: mean" ; strocnx:time_rep = "instantaneous" ; float strocny(time, nj, ni) ; strocny:units = "N/m^2" ; strocny:long_name = "ocean/ice stress (y)" ; strocny:coordinates = "ULON ULAT time" ; strocny:cell_measures = "area: uarea" ; strocny:missing_value = 1.e+30f ; strocny:_FillValue = 1.e+30f ; strocny:cell_methods = "time: mean" ; strocny:time_rep = "instantaneous" ; float strintx(time, nj, ni) ; strintx:units = "N/m^2" ; strintx:long_name = "internal ice stress (x)" ; strintx:coordinates = "ULON ULAT time" ; strintx:cell_measures = "area: uarea" ; strintx:missing_value = 1.e+30f ; strintx:_FillValue = 1.e+30f ; strintx:cell_methods = "time: mean" ; strintx:time_rep = "instantaneous" ; float strinty(time, nj, ni) ; strinty:units = "N/m^2" ; strinty:long_name = "internal ice stress (y)" ; strinty:coordinates = "ULON ULAT time" ; strinty:cell_measures = "area: uarea" ; strinty:missing_value = 1.e+30f ; strinty:_FillValue = 1.e+30f ; strinty:cell_methods = "time: mean" ; strinty:time_rep = "instantaneous" ; float taubx(time, nj, ni) ; taubx:units = "N/m^2" ; taubx:long_name = "basal (seabed) stress (x)" ; taubx:coordinates = "ULON ULAT time" ; taubx:cell_measures = "area: uarea" ; taubx:missing_value = 1.e+30f ; taubx:_FillValue = 1.e+30f ; taubx:cell_methods = "time: mean" ; taubx:time_rep = "instantaneous" ; float tauby(time, nj, ni) ; tauby:units = "N/m^2" ; tauby:long_name = "basal (seabed) stress (y)" ; tauby:coordinates = "ULON ULAT time" ; tauby:cell_measures = "area: uarea" ; tauby:missing_value = 1.e+30f ; tauby:_FillValue = 1.e+30f ; tauby:cell_methods = "time: mean" ; tauby:time_rep = "instantaneous" ; float strength(time, nj, ni) ; strength:units = "N/m" ; strength:long_name = "compressive ice strength" ; strength:coordinates = "TLON TLAT time" ; strength:cell_measures = "area: tarea" ; strength:missing_value = 1.e+30f ; strength:_FillValue = 1.e+30f ; strength:cell_methods = "time: mean" ; strength:time_rep = "instantaneous" ; float divu(time, nj, ni) ; divu:units = "%/day" ; divu:long_name = "strain rate (divergence)" ; divu:coordinates = "TLON TLAT time" ; divu:cell_measures = "area: tarea" ; divu:missing_value = 1.e+30f ; divu:_FillValue = 1.e+30f ; divu:cell_methods = "time: mean" ; divu:time_rep = "instantaneous" ; float shear(time, nj, ni) ; shear:units = "%/day" ; shear:long_name = "strain rate (shear)" ; shear:coordinates = "TLON TLAT time" ; shear:cell_measures = "area: tarea" ; shear:missing_value = 1.e+30f ; shear:_FillValue = 1.e+30f ; shear:cell_methods = "time: mean" ; shear:time_rep = "instantaneous" ; float sig1(time, nj, ni) ; sig1:units = "1" ; sig1:long_name = "norm. principal stress 1" ; sig1:coordinates = "ULON ULAT time" ; sig1:cell_measures = "area: uarea" ; sig1:missing_value = 1.e+30f ; sig1:_FillValue = 1.e+30f ; sig1:cell_methods = "time: mean" ; sig1:time_rep = "instantaneous" ; float sig2(time, nj, ni) ; sig2:units = "1" ; sig2:long_name = "norm. principal stress 2" ; sig2:coordinates = "ULON ULAT time" ; sig2:cell_measures = "area: uarea" ; sig2:missing_value = 1.e+30f ; sig2:_FillValue = 1.e+30f ; sig2:cell_methods = "time: mean" ; sig2:time_rep = "instantaneous" ; float sigP(time, nj, ni) ; sigP:units = "1" ; sigP:long_name = "ice pressure" ; sigP:coordinates = "ULON ULAT time" ; sigP:cell_measures = "area: uarea" ; sigP:missing_value = 1.e+30f ; sigP:_FillValue = 1.e+30f ; sigP:cell_methods = "time: mean" ; sigP:time_rep = "instantaneous" ; float dvidtt(time, nj, ni) ; dvidtt:units = "cm/day" ; dvidtt:long_name = "volume tendency thermo" ; dvidtt:coordinates = "TLON TLAT time" ; dvidtt:cell_measures = "area: tarea" ; dvidtt:missing_value = 1.e+30f ; dvidtt:_FillValue = 1.e+30f ; dvidtt:cell_methods = "time: mean" ; dvidtt:time_rep = "instantaneous" ; float dvidtd(time, nj, ni) ; dvidtd:units = "cm/day" ; dvidtd:long_name = "volume tendency dynamics" ; dvidtd:coordinates = "TLON TLAT time" ; dvidtd:cell_measures = "area: tarea" ; dvidtd:missing_value = 1.e+30f ; dvidtd:_FillValue = 1.e+30f ; dvidtd:cell_methods = "time: mean" ; dvidtd:time_rep = "instantaneous" ; // global attributes: :title = "sea ice model output for CICE" ; :contents = "Diagnostic and Prognostic Variables" ; :source = "Los Alamos Sea Ice Model, undefined_version_name" ; :comment = "All years have exactly 365 days" ; :comment2 = "File written on model date 20011005" ; :comment3 = "seconds elapsed into model date: 7200" ; :conventions = "CF-1.0" ; :history = "This dataset was created on 2021-01-05 at 00:23:45.9" ; :io_flavor = "io_netcdf" ; data: time = 277.083333333333 ; }The file in question is attached, with the namelist and diagnostic output of the run that created it.
I thought for a moment that it was because of the
coordinatesattributes that is recorded as "TLON TLAT time", i.e. in a diffferent order than the variables' shape, which are all(time, nj, ni), but CICE standalone also produces files like these and Panoply does not complain....iceh_inst.2001-10-05-07200.tar.gz
The text was updated successfully, but these errors were encountered: