Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix for time units attribute when IAU used with netcdf format #2

Closed
wants to merge 4 commits into from
Closed

fix for time units attribute when IAU used with netcdf format #2

Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6bb6800
fix for time units attribute when IAU used with netcdf format
jswhit2 Oct 17, 2019
b05be21
modify time units in module_fcst_grid_comp instead of module_write_ne…
jswhit2 Nov 1, 2019
623dc0c
remove unused vars
jswhit2 Nov 1, 2019
0e695bd
Merge remote-tracking branch 'upstream/develop' into feature/nc_time_…
jswhit2 Nov 4, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
73 changes: 70 additions & 3 deletions io/module_write_netcdf.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ module module_write_netcdf

use esmf
use netcdf
use module_fv3_io_def,only : ideflate, nbits, &
use module_fv3_io_def,only : ideflate, nbits, iau_offset, &
output_grid,dx,dy,lon1,lat1,lon2,lat2

implicit none
Expand Down Expand Up @@ -60,8 +60,6 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
real(8) :: varr8val
character(len=ESMF_MAXSTR) :: varcval

character(128) :: time_units

integer :: ncerr
integer :: ncid
integer :: oldMode
Expand Down Expand Up @@ -520,6 +518,10 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc)
real(ESMF_KIND_R4), allocatable :: valueListR4(:)
real(ESMF_KIND_R8), allocatable :: valueListR8(:)
character(len=ESMF_MAXSTR), allocatable :: valueListC(:)
character(nf90_max_name) :: time_units
integer idate(6)
real(8) rinc(5)
integer idat(8),jdat(8)
!
call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", &
attnestflag=ESMF_ATTNEST_OFF, name=dim_name, &
Expand Down Expand Up @@ -555,8 +557,73 @@ subroutine add_dim(ncid, dim_name, dimid, grid, rc)
end if

call get_grid_attr(grid, dim_name, ncid, dim_varid, rc)

! if iau_offset set, update time units attribute.
if ( trim(dim_name) == "time" .and. iau_offset > 0) then
! get time units just written to file
ncerr = nf90_get_att(ncid, dim_varid, 'units', time_units); NC_ERR_STOP(ncerr)
! idate: year,month,day,hour,minute,second
idate = get_idate_from_time_units(time_units)
! idat: year,month,day,time zone,hour,minute,second,millsecond
idat = 0
idat(1:3) = idate(1:3)
idat(5:7) = idate(4:6)
! rinc: days,hours,minutes,seconds,milliseconds
rinc = 0; rinc(2) = iau_offset
call w3movdat(rinc,idat,jdat)
! update idate using iau_offset
idate(1:3)=jdat(1:3)
idate(4:6)=jdat(5:7)
time_units = get_time_units_from_idate(idate)
! rewrite time units attribute.
ncerr = nf90_put_att(ncid, dim_varid, 'units', trim(time_units)); NC_ERR_STOP(ncerr)
endif
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can this manipulation of 'units' attribute be done outside of this subroutine? At the place where data and attributes are added to ESMF_FieldBundle. This routine (write_netcdf) should be as generic as possible. It shouldn't "know" anything about iau_offset.

Copy link
Contributor

@jswhit jswhit Nov 1, 2019 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ideally yes, but I'm afraid that model may use this ESMF units attribute internally. Maybe @junwang-noaa can confirm whether it's used or not?

Copy link
Collaborator

@junwang-noaa junwang-noaa Nov 1, 2019 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I looked at the code again. Write grid component actually uses IO_BASETIME ( ESMF_Time) and IO_CURRTIMEDIFF (ESMF_TimeInterval) in its internal state to specify output initial date and the forecast hours related to initial output date. These two variable can be changed from forecast state due to various reasons, and in current code write grid component already adjusts the two time variables if iau_offset>0. Instead getting the initial output time from write grid, I think you can pass these two time variables into write_netcdf subroutine.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The "time:units" attribute is set in module_fcst_grid_comp.F90. I see iau_offset is available in that module. Can "time:units" attribute be created with correct value in that module and passed to write_netcdf routine via FieldBundle. Instead of writing incorrect value, reading it back, parsing individual time components, adjusting based on iau_offset, constructing new "hours since ..." string and finally overwriting the one already written.
Passing IO_BASETIME and IO_CURRTIMEDIFF instead of iau_offset and doing the same or similar manipulation is equally bad.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"Time:units" in fcst_grid_comp is to control the forecast integration and it is correct to provide time information from forecast point of view (starting integration from 2019103121 for 2019110100 cycle, but we decided only to take forecast results for 2019110100 history output). In my point of view "Time:units" on forecast grid comp should not be changed, because it is write grid comp decided what to output with the data. In the case of running with IAU, it's 2019110100 is the corresponding f00. That is why we use the two time variables IO_BASETIME and IO_CURRTIMEDIFF owned by write grid comp to have consistent time stamp for all the tasks done on write grid comp(output files, inline_post etc).

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The "time:units" is a string attribute added to the export state for the purpose of passing that information (and many other) to the routine that writes that state to a (netcdf) file (convention="NetCDF"). I do not see how it can "control the forecast". Where exactly is "time:units" attribute used, except in the write_netcdf routine? If it's used then that's also not the best way of controlling the forecast. We use various esmf objects (clocks, time, time intervals and alarms) for that, not string attributes.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I should say "time:units" on forecast grid component shows the time stamp on forecast integration, no matter what downstream routine(s) will use it, we should not change it. Write grid comp should be allowed to make decision on what to output.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I implemented Dusan's suggestion and it seems to work.


end subroutine add_dim

function get_idate_from_time_units(time_units) result(idate)
! return integer array with year,month,day,hour,minute,second
! parsed from time units attribute.
integer idate(6)
character(len=nf90_max_name), intent(in) :: time_units
integer ipos1,ipos2
ipos1 = scan(time_units,"since",back=.true.)+1
ipos2 = scan(time_units,"-",back=.false.)-1
read(time_units(ipos1:ipos2),*) idate(1)
ipos1 = ipos2+2; ipos2=ipos1+1
read(time_units(ipos1:ipos2),*) idate(2)
ipos1 = ipos2+2; ipos2=ipos1+1
read(time_units(ipos1:ipos2),*) idate(3)
ipos1 = scan(time_units,":")-2
ipos2 = ipos1+1
read(time_units(ipos1:ipos2),*) idate(4)
ipos1 = ipos2+2
ipos2 = ipos1+1
read(time_units(ipos1:ipos2),*) idate(5)
ipos1 = ipos2+2
ipos2 = ipos1+1
read(time_units(ipos1:ipos2),*) idate(6)
end function get_idate_from_time_units

function get_time_units_from_idate(idate, time_measure) result(time_units)
! create time units attribute of form 'hours since YYYY-MM-DD HH:MM:SS'
! from integer array with year,month,day,hour,minute,second
! optional argument 'time_measure' can be used to change 'hours' to
! 'days', 'minutes', 'seconds' etc.
character(len=*), intent(in), optional :: time_measure
integer, intent(in) :: idate(6)
character(len=12) :: timechar
character(len=nf90_max_name) :: time_units
if (present(time_measure)) then
timechar = trim(time_measure)
else
timechar = 'hours'
endif
write(time_units,101) idate
101 format(' since ',i4.4,'-',i2.2,'-',i2.2,' ',&
i2.2,':',i2.2,':',i2.2)
time_units = trim(adjustl(timechar))//time_units
end function get_time_units_from_idate
!
!----------------------------------------------------------------------------------------
subroutine nccheck(status)
Expand Down