[CLM-only] Restart files do NOT have _FillValue/missing_value attributes on fields #88
Comments
ekluzek
added
the
priority: low
|
Erik Kluzek < erik > - 2011-08-24 14:12:36 -0600 mkarbinit sets t_soisno to spval and then defines the non-snow layers to something for everything besides lake. Hydrology2Mod then has the following... But, for t_soisno(c,j) it should set it to spval rather than zero. |
|
Tim Hoar < thoar > - 2011-08-24 14:28:23 -0600 There may be more to it ... I'm looking at values for H2OSOI_LIQ and sometimes 0, // H2OSOI_LIQ(1,1602) |
|
Erik Kluzek < erik > - 2011-08-24 14:47:41 -0600 (In reply to comment #2)
Actually, that's a good point. T_SOISNO physically can't be zero, so it should be set to spval, but the others can be zero so I thought they should be left alone. But, really the zero is being used to mark them as missing values NOT to a literal value of zero. The other values are also initialized to spval, and then set to something in mkarbinit. So here (in Hydrology2) they should all be set to spval and _FillValue should equal spval on the file. I also found several places where 1.e36 was hardcoded in rather than using spval... biogeochem/DryDepVelocity.F90: crs = 1.e36_r8 |
|
Erik Kluzek < erik > - 2011-08-25 13:34:38 -0600 Tim also notes some values where snow isn't where T_SOISNO is set to a negative value. Presumably, this is because it's being set to zero, and then there is an operation over it that should be checking SNLSNO that isn't... Keith is going to try to find out where this is happening.
|
|
Keith Oleson < oleson > - 2011-08-25 16:44:35 -0600 I was able to replicate YongFei's case bfb except that I did not get negative values for T_SOISNO in the restart files. YongFei now reports that her negative values are probably due to the assimilation procedure. Regardless, I think that the variables mentioned by Erik below could probably be set to spval instead of zero. I did a quick test where I set them to spval in Hydrology2 and got bfb answers. This is very limited testing of course. Zeros can still appear in the restart file however if the initial file used for the simulation has zeros and the snow layers don't change state during the simulation. |
|
Erik Kluzek < erik > - 2011-09-21 12:12:20 -0600 I tried setting fill_value=nan for several fields on the restart file that are initialized to nan. This worked except it blew up in DEBUG mode. So I took it back out. So either fields that will be on the restart file need to be initialized to spval, or I need to convert them from nan to spval when writing and from spval to nan when reading. Consistently initializing fields to nan and only changing them to spval for I/O would be a good thing. |
|
Erik Kluzek < erik > - 2011-09-27 12:04:32 -0600 Fields that Tim would like us to especially pay attention to... Erik, Dave; As soon as these variables have a _FillValue (?and a missing_value?) attribute, I'm happy to check the range of values to see if there are (perhaps) unanticipated uses of other fill values. I'll be working on Matlab-based diagnostic scripts, so I will be doing this in parallel to Erik's effort. To all on the 'cc' list; the CLM restart files are easier to ingest into DART if they have a _FillValue specified. If you think you want particular variables to be part of the DART state vector - speak up now. [Andy, Bill - a' priori knowledge of the valid_range might be enough. I don't think we're going to get CLM to add _FillValue and (perhaps) valid_range attributes to the 5kajillion carbon species.] Here is my list of 11 CLM variables that could use netcdf attributes of '_FillValue'. Tim Hoar, Associate Scientist |
|
Erik Kluzek < erik > - 2011-09-27 12:05:28 -0600 This was marked as fixed and shouldn't have been... |
|
Erik Kluzek < erik > - 2011-09-27 14:59:16 -0600 Of the fields in the list only: T_SOISNO, H2OSOI_LIQ/ICE are set to spval everything else is set to nan. |
|
Tim Hoar < thoar > - 2011-09-28 10:55:00 -0600 (In reply to comment #9)
OK then ... I need to be able to replace nan with my DART missing code, ... do what I need to do ... |
|
Erik Kluzek < erik > - 2012-01-23 10:01:40 -0700 I'm starting to work on this in clm4_0_38. I checked in code to handle it and it works on bluefire, but fails with intel. So I'm going to back it out until I can get a solution that works on all platforms. I am leaving in an option to ncdio_pio cnvrtnan2fill that will convert any spval read to nan, and convert any nan on write to spval so that FillValue functions correctly. Here's the code changes to get back to the previous version. Without the changes below, the convrtnan2fill option is never used. [yong:lnd/clm/src] erik% svn diff
Index: main/accumulMod.F90
===================================================================
--- main/accumulMod.F90 (revision 33981)
+++ main/accumulMod.F90 (working copy)
@@ -27,7 +27,6 @@
use abortutils , only: endrun
use clm_varctl , only: iulog
use nanMod , only: bigint
- use shr_infnan_mod, only: shr_infnan_isnan
!
! !PUBLIC TYPES:
implicit none
@@ -425,7 +424,7 @@
! !IROUTINE: update_accum_field_sl
!
! !INTERFACE:
- subroutine update_accum_field_sl (name, field, nstep, filternan)
+ subroutine update_accum_field_sl (name, field, nstep)
!
! !DESCRIPTION:
! Accumulate single level field over specified time interval.
@@ -436,7 +435,6 @@
character(len=*), intent(in) :: name !field name
real(r8), pointer, dimension(:) :: field !field values for current time step
integer , intent(in) :: nstep !time step index
- logical , intent(in), optional :: filternan ! filter out any NaN's in the input field
!
! !REVISION HISTORY:
! Created by Sam Levis
@@ -449,8 +447,6 @@
integer :: i,k,nf !indices
integer :: accper !temporary accumulation period
integer :: beg,end !subgrid beginning,ending indices
- integer :: filter(sizeof(field)) !filter of
- integer :: nfilter ! number of points to filter
!------------------------------------------------------------------------
! find field index. return if "name" is not on list
@@ -475,23 +471,6 @@
call endrun
endif
- ! Get filter
- if ( filternan ) then
- nfilter = 0
- do i = beg, end
- if ( .not. shr_infnan_isnan( field(i) ) )then
- nfilter = nfilter + 1
- filter(nfilter) = i
- end if
- end do
- else
- nfilter = 0
- do i = beg, end
- nfilter = nfilter + 1
- filter(nfilter) = i
- end do
- end if
-
! accumulate field
if (accum(nf)%acctype /= 'timeavg' .AND. &
@@ -515,11 +494,11 @@
!normalize at end of accumulation period
if ((mod(nstep,accum(nf)%period) == 1) .and. (nstep /= 0)) then
- accum(nf)%val(filter(:nfilter),1) = 0._r8
+ accum(nf)%val(beg:end,1) = 0._r8
end if
- accum(nf)%val(filter(:nfilter),1) = accum(nf)%val(filter(:nfilter),1) + field(filter(:nfilter))
+ accum(nf)%val(beg:end,1) = accum(nf)%val(beg:end,1) + field(beg:end)
if (mod(nstep,accum(nf)%period) == 0) then
- accum(nf)%val(filter(:nfilter),1) = accum(nf)%val(filter(:nfilter),1) / accum(nf)%period
+ accum(nf)%val(beg:end,1) = accum(nf)%val(beg:end,1) / accum(nf)%period
endif
else if (accum(nf)%acctype == 'runmean') then
@@ -527,19 +506,18 @@
!running mean - reset accumulation period until greater than nstep
accper = min (nstep,accum(nf)%period)
- accum(nf)%val(filter(:nfilter),1) = ((accper-1)*accum(nf)%val(filter(:nfilter),1) + field(filter(:nfilter))) / accper
+ accum(nf)%val(beg:end,1) = ((accper-1)*accum(nf)%val(beg:end,1) + field(beg:end)) / accper
else if (accum(nf)%acctype == 'runaccum') then
!running accumulation field reset at trigger -99999
- do i = 1, nfilter
- k = filter(i)
+ do k = beg,end
if (nint(field(k)) == -99999) then
accum(nf)%val(k,1) = 0._r8
end if
end do
- accum(nf)%val(filter(:nfilter),1) = min(max(accum(nf)%val(filter(:nfilter),1) + field(filter(:nfilter)), 0._r8), 99999._r8)
+ accum(nf)%val(beg:end,1) = min(max(accum(nf)%val(beg:end,1) + field(beg:end), 0._r8), 99999._r8)
end if
Index: main/accFldsMod.F90
===================================================================
--- main/accFldsMod.F90 (revision 33981)
+++ main/accFldsMod.F90 (working copy)
@@ -578,9 +578,9 @@
do p = begp,endp
rbufslp(p) = fsun(p)
end do
- call update_accum_field ('FSUN24', rbufslp, nstep, filternan=.true.)
- call extract_accum_field ('FSUN24', fsun24, nstep )
- call update_accum_field ('FSUN240', rbufslp, nstep, filternan=.true.)
+ call update_accum_field ('FSUN24', rbufslp, nstep)
+ call extract_accum_field ('FSUN24', fsun24, nstep)
+ call update_accum_field ('FSUN240', rbufslp, nstep)
call extract_accum_field ('FSUN240', fsun240, nstep)
! Accumulate and extract elai_p (heald 04/06)
Index: biogeophys/BiogeophysRestMod.F90
===================================================================
--- biogeophys/BiogeophysRestMod.F90 (revision 33981)
+++ biogeophys/BiogeophysRestMod.F90 (working copy)
@@ -59,6 +59,7 @@
use initSurfAlbMod , only : do_initsurfalb
use clm_time_manager, only : is_first_step
use SNICARMod , only : snw_rds_min
+ use shr_infnan_mod , only : shr_infnan_isnan
use mkarbinitMod , only : perturbIC
use clm_time_manager, only : is_restart
!
@@ -345,12 +346,12 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='frac_sno', xtype=ncd_double, &
- dim1name='column', fill_value=spval, &
+ dim1name='column',&
long_name='fraction of ground covered by snow (0 to 1)',units='unitless')
else if (flag == 'read' .or. flag == 'write') then
call ncd_io(varname='frac_sno', data=cptr%cps%frac_sno, &
dim1name=namec, &
- ncid=ncid, flag=flag, readvar=readvar, cnvrtnan2fill=.true.)
+ ncid=ncid, flag=flag, readvar=readvar)
if (flag == 'read' .and. .not. readvar) then
if (is_restart()) call endrun()
end if
@@ -361,14 +362,14 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='DZSNO', xtype=ncd_double, &
dim1name='column', dim2name='levsno', switchdim=.true., &
- long_name='snow layer thickness', units='m', fill_value=spval)
+ long_name='snow layer thickness', units='m')
else if (flag == 'read' .or. flag == 'write') then
allocate(temp2d(begc:endc,-nlevsno+1:0))
if (flag == 'write') then
temp2d(begc:endc,-nlevsno+1:0) = cptr%cps%dz(begc:endc,-nlevsno+1:0)
end if
call ncd_io(varname='DZSNO', data=temp2d, &
- dim1name=namec, switchdim=.true., cnvrtnan2fill=.true., &
+ dim1name=namec, switchdim=.true., &
lowerb2=-nlevsno+1, upperb2=0, ncid=ncid, flag=flag, readvar=readvar)
if (flag == 'read' .and. .not. readvar) then
if (is_restart()) call endrun()
@@ -384,14 +385,14 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='ZSNO', xtype=ncd_double, &
dim1name='column', dim2name='levsno', switchdim=.true., &
- long_name='snow layer depth', units='m', fill_value=spval)
+ long_name='snow layer depth', units='m')
else if (flag == 'read' .or. flag == 'write') then
allocate(temp2d(begc:endc,-nlevsno+1:0))
if (flag == 'write') then
temp2d(begc:endc,-nlevsno+1:0) = cptr%cps%z(begc:endc,-nlevsno+1:0)
end if
call ncd_io(varname='ZSNO', data=temp2d, &
- dim1name=namec, switchdim=.true., cnvrtnan2fill=.true., &
+ dim1name=namec, switchdim=.true., &
lowerb2=-nlevsno+1, upperb2=0, ncid=ncid, flag=flag, readvar=readvar)
if (flag == 'read' .and. .not. readvar) then
if (is_restart()) call endrun()
@@ -407,14 +408,14 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='ZISNO', xtype=ncd_double, &
dim1name='column', dim2name='levsno', switchdim=.true., &
- long_name='snow interface depth', units='m', fill_value=spval)
+ long_name='snow interface depth', units='m')
else if (flag == 'read' .or. flag == 'write') then
allocate(temp2d(begc:endc,-nlevsno:-1))
if (flag == 'write') then
temp2d(begc:endc,-nlevsno:-1) = cptr%cps%zi(begc:endc,-nlevsno:-1)
end if
call ncd_io(varname='ZISNO', data=temp2d, &
- dim1name=namec, switchdim=.true., cnvrtnan2fill=.true., &
+ dim1name=namec, switchdim=.true., &
lowerb2=-nlevsno, upperb2=-1, ncid=ncid, flag=flag, readvar=readvar)
if (flag == 'read' .and. .not. readvar) then
if (is_restart()) call endrun()
@@ -928,13 +929,13 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='H2OSOI_LIQ', xtype=ncd_double, &
dim1name='column', dim2name='levtot', switchdim=.true., &
- long_name='liquid water', units='kg/m2', fill_value=spval)
+ long_name='liquid water', units='kg/m2')
else if (flag == 'read' .or. flag == 'write') then
call ncd_io(varname='H2OSOI_LIQ', data=cptr%cws%h2osoi_liq, &
- dim1name=namec, switchdim=.true., ncid=ncid, flag=flag, readvar=readvar)
- if (flag=='read' .and. .not. readvar) then
- if (is_restart()) call endrun()
- end if
+ dim1name=namec, switchdim=.true., ncid=ncid, flag=flag, readvar=readvar)
+ if (flag=='read' .and. .not. readvar) then
+ if (is_restart()) call endrun()
+ end if
end if
! column water state variable - h2osoi_ice
@@ -942,8 +943,7 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='H2OSOI_ICE', xtype=ncd_double, &
dim1name='column', dim2name='levtot', switchdim=.true., &
- long_name='ice lens', units='kg/m2', fill_value=spval)
-
+ long_name='ice lens', units='kg/m2')
else if (flag == 'read' .or. flag == 'write') then
call ncd_io(varname='H2OSOI_ICE', data=cptr%cws%h2osoi_ice, &
dim1name=namec, switchdim=.true., ncid=ncid, flag=flag, readvar=readvar)
@@ -956,7 +956,7 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='T_GRND', xtype=ncd_double, &
- dim1name='column', fill_value=spval, &
+ dim1name='column', &
long_name='ground temperature', units='K')
else if (flag == 'read' .or. flag == 'write') then
call ncd_io(varname='T_GRND', data=cptr%ces%t_grnd, &
@@ -1242,7 +1242,7 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='T_SOISNO', xtype=ncd_double, &
dim1name='column', dim2name='levtot', switchdim=.true., &
- long_name='soil-snow temperature', units='K', fill_value=spval)
+ long_name='soil-snow temperature', units='K')
else if (flag == 'read' .or. flag == 'write') then
call ncd_io(varname='T_SOISNO', data=cptr%ces%t_soisno, &
dim1name=namec, switchdim=.true., ncid=ncid, flag=flag, readvar=readvar)
@@ -1256,7 +1256,7 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='T_LAKE', xtype=ncd_double, &
dim1name='column', dim2name='levlak', switchdim=.true., &
- long_name='lake temperature', units='K', fill_value=spval)
+ long_name='lake temperature', units='K')
else if (flag == 'read' .or. flag == 'write') then
call ncd_io(varname='T_LAKE', data=cptr%ces%t_lake, &
dim1name=namec, switchdim=.true., ncid=ncid, flag=flag, readvar=readvar)
@@ -1372,15 +1372,21 @@
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname='fsun', xtype=ncd_double, &
- dim1name='pft', fill_value=spval, &
+ dim1name='pft', &
long_name='sunlit fraction of canopy', units='')
else if (flag == 'read' .or. flag == 'write') then
call ncd_io(varname='fsun', data=pptr%pps%fsun, &
dim1name=namep, &
- ncid=ncid, flag=flag, readvar=readvar, cnvrtnan2fill=.true. )
+ ncid=ncid, flag=flag, readvar=readvar)
if (flag=='read' )then
if ( .not. readvar) then
if (is_restart()) call endrun()
+ else
+ do p = begp, endp
+ if ( shr_infnan_isnan( pptr%pps%fsun(p) ) )then
+ pptr%pps%fsun(p) = spval
+ end if
+ end do
end if
end if
end if
@@ -1505,7 +1511,7 @@
varname = 'T_REF2M'
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname=varname, xtype=ncd_double, &
- dim1name='pft', fill_value=spval, &
+ dim1name='pft', &
long_name='2m height surface air temperature', units='K')
else if (flag == 'read' .or. flag == 'write') then
call ncd_io(varname=varname, data=pptr%pes%t_ref2m, & |
|
Erik Kluzek < erik > - 2012-01-23 10:16:03 -0700 Here's the error on mirage... (seq_mct_drv) : Model initialization complete forrtl: severe (174): SIGSEGV, segmentation fault occurred This is the update of TREVAV, in assessing accum(nf)%acctype... |
billsacks
added
the
enhancement
|
I was about to open a new issue with title Here's a summary of the issue that I was about to post: Suggestions? Ideas? Comments? |
|
@slevisconsulting Thank you for your continued careful work on this, which is leading to the discovery of potential issues. I can't tell from your comment: are the 0 values actually causing problems, or are they possibly in unused points (inactive columns or non-existent soil / snow layers)? Before deciding on a solution to the problem, I think it would help to better understand the cause. Is this happening in snow layers or soil layers? Over what landunit/column type(s)? And what aspect of the interpolation – horizontal or vertical – is creating these 0 values? It looks like t_soisno is initialized to spval everywhere in this code: CTSM/src/biogeophys/TemperatureType.F90 Lines 685 to 688 in ff1ae2c The only place where I see an explicit zeroing of t_soisno is here: CTSM/src/biogeophys/SnowHydrologyMod.F90 Line 2892 in ff1ae2c So is it just that these zeroed snow layers get copied into the output file? If so, have you checked whether the 0 values you're seeing are only in non-existent snow layers (based on snlsno in the given column)? I believe that it is expected that these 0 values would make their way into the interpolated finidat file, but they should only appear in non-existent snow layers. If you're seeing something different, we should talk more. If this does turn out to be a problem, then putting variable-specific logic in init_interp should be a last resort solution: Mariana and I worked hard to remove this variable-specific logic from init_interp a few years ago because it was hard to maintain and led to bugs. One solution might be to find where t_soisno is set to 0 in the code (maybe just the one above line) and set it to spval there instead, if that would solve the problem. |
|
The near-0K TSOI affects the calculated FSH, so I think we have an actual problem:
The near-0K TSOI values appear down to the deepest soil layer before we get into the bedrock. Line 2892 that you referenced is in the loop Replacing all T_SOISNO=0 with 283 or NaN or 1e36 directly in the finidat file fixed the problem, so I think that init_interp receives 0K and treats it as a valid value. To avoid variable-specific logic, how about unit-specific logic that applies to degrees K? |
|
@billsacks' comment during this morning's mtg suggests to me that the T_SOISNO=0 values end up in the WRF-CTSM(CLMSP) finidat file after the manipulation with the .ncl script (here: /glade/work/slevis/git_wrf/ctsm_init) by which we introduce WRF initial conditions to finidat. Looking at the .ncl script tells me that we do not explicitly and intentionally initialize T_SOISNO=0 (except for snow layers), so I suspect that T_SOISNO=0 enters via WRF's TSLB variable. For now it is clear that I may proceed with WRF-CTSM(CLMSP) simulations once I ensure that T_SOISNO=0 is replaced with NaNs in my finidat. Whether init_interp's handling of 0K should be changed has not been resolved. |
|
I checked 7 of our out-of-the-box initial conditions files with this python code (depends on numpy and Netcdf4, both available on cheyenne if you load def print_stats(filename):
dat = Dataset(filename)
tsoi = dat.variables['T_SOISNO'][:]
levsno = dat.dimensions['levsno'].size
print('levsno = {}'.format(levsno))
print('number of points with tsoi == 0: {}'.format(np.sum(tsoi == 0)))
print('number of snow points: {}'.format(np.sum(tsoi[:,:levsno] == 0)))
print('number of soil points: {}'.format(np.sum(tsoi[:,levsno:] == 0)))
print('min T in soil: {}'.format(np.min(tsoi[:,levsno:])))All files had a lot of 0 values over snow layers (as expected), but none of them had any 0 values in any soil layers. It would be worth double-checking this for the original source file you started with before merging in the WRF values. As long as you see the same thing, then I think the culprit is the WRF variable, and maybe the solution is to fix this in the ncl script or as a post-processing step on this merged file, rather than changing any CTSM code to handle this scenario. If you do add any CTSM code, what I'd suggest is adding an error check in the restart routine in TemperatureType that ensures that there are no super-cold temperature values after reading the initial conditions file, and aborting if there are. The minimum values on our out-of-the-box files are about 210 K, but probably a check for any values <= 0 would catch any common issues (with 0 values or unit problems). To be clear: I don't feel this is necessary, but I'd be happy with this error checking being added if you feel it would help catch future errors early. |
|
Thanks @billsacks . |
### Description of changes Improve support for user mods in *.streams.xml file. If the file is found in the run directory it will be used for that data components streams. I have added an xsd file as well to use with xmllint to assure the syntax of the file. Due to namelist change in cmeps for multiple ice sheets all namelists changed and so new baselines, may12, were generated. ### Specific notes Contributors other than yourself, if any: CMEPS Issues Fixed (include github issue #): ESCOMP#79 Are there dependencies on other component PRs - [ ] CIME (list) - [ ] CMEPS (list) Are changes expected to change answers? - [X] bit for bit - [ ] different at roundoff level - [ ] more substantial Any User Interface Changes (namelist or namelist defaults changes)? - [X] Yes Indirectly due to changes in glc (multiple ice sheet support) - [ ] No Testing performed: - [X] (required) aux_cdeps - machines and compilers: cheyenne intel - details (e.g. failed tests): see issue ESCOMP#89 - [ ] (optional) CESM prealpha test - machines and compilers - details (e.g. failed tests): Hashes used for testing: - [X] CIME - repository to check out: https://github.com/ESCOMP/CESM.git - branch: master - hash: 353730e - [X] CMEPS - repository to check out: https://github.com/ESCOMP/CESM.git - branch: master - hash: cf0d47ba - [ ] CESM - repository to check out: https://github.com/ESCOMP/CESM.git - branch: nuopc_dev - hash:
### Description of changes These should have been part of ESCOMP#88 and were used in the testing of that PR but missed in the commit.
Erik Kluzek < erik > - 2011-08-24 14:02:01 -0600
Bugzilla Id: 1401
Bugzilla CC: dlawren, oleson, rfisher, sacks, thoar,
Tim Hoar noticed that the restart files don't have _FillValue/missing_value attributes on fields as they should. The history files do, but restart files do NOT. Also we noticed that T_SOISNO has some values set to 0.0 while others are set to 1.e+36. Both of which are really missing values, but the code shouldn't set them to two different values.
The text was updated successfully, but these errors were encountered: