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

Add support for annual data in interpolator_mod #1389

Merged
merged 7 commits into from
Dec 14, 2023
Merged

Add support for annual data in interpolator_mod #1389

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
738981a
add yearly support to interpolator
Oct 11, 2023
158ed5c
unit test cleanup
Oct 19, 2023
b61e1af
cleanup interpolator.inc
Oct 19, 2023
a3b03a5
fix mistake
Oct 19, 2023
19a9788
remove write statements for debugging
Oct 19, 2023
e3c9e5c
remove period
Dec 13, 2023
ffadb8c
Merge remote-tracking branch 'origin/main' into annual_interp
Dec 13, 2023
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
149 changes: 95 additions & 54 deletions interpolator/include/interpolator.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,14 +145,23 @@ integer :: model_calendar
integer :: yr, mo, dy, hr, mn, sc
integer :: n
type(time_type) :: Julian_time, Noleap_time
real(FMS_INTP_KIND_), allocatable :: time_in(:)
real(r8_kind), allocatable :: time_in(:)
real(FMS_INTP_KIND_), allocatable, save :: agrid_mod(:,:,:)
integer :: nx, ny
type(FmsNetcdfFile_t) :: fileobj

integer, parameter :: lkind=FMS_INTP_KIND_
real(FMS_INTP_KIND_), parameter :: lPI=real(PI,FMS_INTP_KIND_)

!> variables used to set time
logical :: yearly !< flags to indicate if time data is in units of months or years
integer :: num_years !< number of years
integer :: base_year, base_month, base_day, base_hour, base_minute, base_second
integer :: nn !< counter
logical :: noleap_file_calendar !< is the file calendar noleap or julian
real(r8_kind) :: num_days, frac_year !< variables for yearly=.true.
type(time_type) :: n_time !< temporary time

clim_type%separate_time_vary_calc = .false.

num_fields = 0
Expand Down Expand Up @@ -332,6 +341,7 @@ if(dimension_exists(fileobj, "time")) then
filehr = 0
filemin = 0
filesec = 0
yearly = .false.
select case(units(:3))
case('day')
fileunits = units(12:) !Assuming "days since YYYY-MM-DD HH:MM:SS"
Expand Down Expand Up @@ -363,8 +373,24 @@ if(dimension_exists(fileobj, "time")) then
read(fileunits(12:13), *) filehr
read(fileunits(15:16), *) filemin
read(fileunits(18:19), *) filesec
case( 'yea')
fileunits = units(13:) !Assuming "years since YYYY-MM-DD HH:MM:SS"
if ( len_trim(fileunits) < 19 ) then
write(error_mesg, '(A49,A,A51,A)' ) &
'Interpolator_init : Incorrect time units in file ', &
trim(file_name), '. Expecting years since YYYY-MM-DD HH:MM:SS, found', &
trim(units)
call mpp_error(FATAL,error_mesg)
endif
read(fileunits(1:4) , *) fileyr
read(fileunits(6:7) , *) filemon
read(fileunits(9:10) , *) fileday
read(fileunits(12:13), *) filehr
read(fileunits(15:16), *) filemin
read(fileunits(18:19), *) filesec
yearly = .true.
case default
call mpp_error(FATAL,'Interpolator_init : Time units not recognised in file '//file_name)
call mpp_error(FATAL,'Interpolator_init : Time units not recognized in file '//file_name)
end select

clim_type%climatological_year = (fileyr == 0)
Expand All @@ -375,32 +401,35 @@ if(dimension_exists(fileobj, "time")) then
! if file date has a non-zero year in the base time, determine that
! base_time based on the netcdf info.
!----------------------------------------------------------------------
if ( (model_calendar == JULIAN .and. &
& trim(adjustl(lowercase(file_calendar))) == 'julian') .or. &
& (model_calendar == NOLEAP .and. &
& trim(adjustl(lowercase(file_calendar))) == 'noleap') ) then
call mpp_error (NOTE, 'interpolator_mod: Model and file&
& calendars are the same for file ' // &
& trim(file_name) // '; no calendar conversion &
&needed')
base_time = set_date (fileyr, filemon, fileday, filehr, &
filemin,filesec)
else if ( (model_calendar == JULIAN .and. &
& trim(adjustl(lowercase(file_calendar))) == 'noleap')) then
call mpp_error (NOTE, 'interpolator_mod: Using julian &
noleap_file_calendar=.false.
if ( (model_calendar == JULIAN .and. trim(adjustl(lowercase(file_calendar))) == 'julian')) then
call mpp_error (NOTE, 'interpolator_mod: Model and file&
& calendars are the same for file ' // &
& trim(file_name) // '; no calendar conversion &
&needed')
base_time = set_date (fileyr, filemon, fileday, filehr,filemin,filesec)
noleap_file_calendar=.false.
else if((model_calendar == NOLEAP .and. trim(adjustl(lowercase(file_calendar))) == 'noleap') ) then
call mpp_error (NOTE, 'interpolator_mod: Model and file&
& calendars are the same for file ' // &
& trim(file_name) // '; no calendar conversion &
&needed')
base_time = set_date (fileyr, filemon, fileday, filehr,filemin,filesec)
noleap_file_calendar=.true.
else if ( (model_calendar == JULIAN .and. trim(adjustl(lowercase(file_calendar))) == 'noleap')) then
call mpp_error (NOTE, 'interpolator_mod: Using julian &
&model calendar and noleap file calendar&
& for file ' // trim(file_name) // &
&'; calendar conversion needed')
base_time = set_date_no_leap (fileyr, filemon, fileday, &
& filehr, filemin, filesec)
else if ( (model_calendar == NOLEAP .and. &
& trim(adjustl(lowercase(file_calendar))) == 'julian')) then
base_time = set_date_no_leap (fileyr, filemon, fileday,filehr, filemin, filesec)
noleap_file_calendar=.true.
else if ( (model_calendar == NOLEAP .and. trim(adjustl(lowercase(file_calendar))) == 'julian')) then
call mpp_error (NOTE, 'interpolator_mod: Using noleap &
&model calendar and julian file calendar&
& for file ' // trim(file_name) // &
&'; calendar conversion needed')
base_time = set_date_julian (fileyr, filemon, fileday, &
& filehr, filemin, filesec)
base_time = set_date_julian (fileyr, filemon, fileday,filehr, filemin, filesec)
noleap_file_calendar=.false.
else
call mpp_error (FATAL , 'interpolator_mod: Model and file&
& calendars ( ' // trim(file_calendar) // ' ) differ &
Expand All @@ -425,17 +454,38 @@ if(dimension_exists(fileobj, "time")) then
if (ntime > 0) then
allocate(time_in(ntime), clim_type%time_slice(ntime))
allocate(clim_type%clim_times(12,(ntime+11)/12))
time_in = 0.0_lkind
time_in = 0.0_r8_kind
clim_type%time_slice = set_time(0,0) + base_time
clim_type%clim_times = set_time(0,0) + base_time
call fms2_io_read_data(fileobj, "time", time_in)
ntime_in = ntime
!> convert the number of years passed to days if yearly=.true.
if(yearly) then
do n = 1, ntime
if(.not.noleap_file_calendar) then
! Julian calendar
num_years = int(time_in(n))
frac_year = time_in(n) - real(num_years, r8_kind)
call get_date_julian(base_time, base_year, base_month, base_day, base_hour, base_minute, base_second)
num_days = 0.0_r8_kind
do nn=1, num_years
if( mod(base_year+nn-1,4)==0) num_days = num_days + 366._r8_kind
rem1776 marked this conversation as resolved.
Show resolved Hide resolved
if( mod(base_year+nn-1,4)/=0) num_days = num_days + 365._r8_kind
end do
if( mod(base_year+num_years,4)==0) num_days = num_days + 366._r8_kind*frac_year
if( mod(base_year+num_years,4)/=0) num_days = num_days + 365._r8_kind*frac_year
else
num_days = time_in(n)*365._r8_kind
end if
time_in(n)=num_days
end do
end if
! determine whether the data is a continuous set of monthly values or
! a series of annual cycles spread throughout the period of data
non_monthly = .false.
do n = 1, ntime-1
! Assume that the times in the data file correspond to days only.
if (time_in(n+1) > (time_in(n) + 32._lkind)) then
if (time_in(n+1) > (time_in(n) + 32._r8_kind)) then
non_monthly = .true.
exit
endif
Expand All @@ -458,9 +508,8 @@ if(dimension_exists(fileobj, "time")) then
!! time_interp_list with the optional argument modtime=YEAR, so that
!! the time that is needed in time_slice is the displacement into the
!! year, not the displacement from a base_time.
clim_type%time_slice(n) = &
set_time( INT( (time_in(n)-real(INT(time_in(n)),FMS_INTP_KIND_)) &
*real(SECONDS_PER_DAY,FMS_INTP_KIND_)), INT(time_in(n)))
clim_type%time_slice(n) = set_time( INT( (time_in(n)-real(INT(time_in(n)),r8_kind))*SECONDS_PER_DAY), &
INT(time_in(n)))
else

!--------------------------------------------------------------------
Expand All @@ -472,59 +521,51 @@ if(dimension_exists(fileobj, "time")) then
! include the base_time; values will be generated relative to the
! "real" time.
!--------------------------------------------------------------------
if ( (model_calendar == JULIAN .and. &
& trim(adjustl(lowercase(file_calendar))) == 'julian') .or. &
& (model_calendar == NOLEAP .and. &
& trim(adjustl(lowercase(file_calendar))) == 'noleap') ) then
if ( (model_calendar == JULIAN .and. trim(adjustl(lowercase(file_calendar))) == 'julian') .or. &
& (model_calendar == NOLEAP .and. trim(adjustl(lowercase(file_calendar))) == 'noleap') ) then

!---------------------------------------------------------------------
! no calendar conversion needed.
!---------------------------------------------------------------------
clim_type%time_slice(n) = &
set_time( INT( (time_in(n)-real(INT(time_in(n)),FMS_INTP_KIND_)) &
*real(SECONDS_PER_DAY,FMS_INTP_KIND_)), INT(time_in(n))) &
+ base_time
set_time( INT( (time_in(n)-real(INT(time_in(n)),r8_kind))*SECONDS_PER_DAY), &
INT(time_in(n))) + base_time

!---------------------------------------------------------------------
! convert file times from noleap to julian.
!---------------------------------------------------------------------
else if ( (model_calendar == JULIAN .and. &
& trim(adjustl(lowercase(file_calendar))) == 'noleap')) then
Noleap_time = set_time (0, INT(time_in(n))) + base_time
call get_date_no_leap (Noleap_time, yr, mo, dy, hr, &
mn, sc)
clim_type%time_slice(n) = set_date_julian (yr, mo, dy, &
hr, mn, sc)
else if ( (model_calendar == JULIAN .and. trim(adjustl(lowercase(file_calendar))) == 'noleap')) then
Noleap_time = set_time( INT((time_in(n)-real(INT(time_in(n)),r8_kind))*SECONDS_PER_DAY), &
Copy link
Contributor Author

Choose a reason for hiding this comment

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

Changed to also take into account fractional number of days (converted to seconds) when setting time for greater accuracy. The same is done for the case where the model calendar and the file calendar are the same (see above).

INT(time_in(n))) + base_time
!Noleap_time = set_time (0, INT(time_in(n))) + base_time
call get_date_no_leap (Noleap_time, yr, mo, dy, hr, mn, sc)
clim_type%time_slice(n) = set_date_julian (yr, mo, dy, hr, mn, sc)
if (n == 1) then
call print_date (clim_type%time_slice(1), &
str= 'for file ' // trim(file_name) // ', the &
&first time slice is mapped to :')
str= 'for file ' // trim(file_name) // ', the first time slice is mapped to :')
endif
if (n == ntime) then
call print_date (clim_type%time_slice(ntime), &
str= 'for file ' // trim(file_name) // ', the &
&last time slice is mapped to:')
str= 'for file ' // trim(file_name) // ', the last time slice is mapped to:')
endif


!---------------------------------------------------------------------
! convert file times from julian to noleap.
!---------------------------------------------------------------------
else if ( (model_calendar == NOLEAP .and. &
& trim(adjustl(lowercase(file_calendar))) == 'julian')) then
Julian_time = set_time (0, INT(time_in(n))) + base_time
else if ( (model_calendar == NOLEAP .and. trim(adjustl(lowercase(file_calendar))) == 'julian')) then
Julian_time = set_time( INT( (time_in(n)-real(INT(time_in(n)),r8_kind))*SECONDS_PER_DAY), &
INT(time_in(n))) + base_time
!Julian_time = set_time (0, INT(time_in(n))) + base_time
call get_date_julian (Julian_time, yr, mo, dy, hr, mn, sc)
clim_type%time_slice(n) = set_date_no_leap (yr, mo, dy, &
hr, mn, sc)
clim_type%time_slice(n) = set_date_no_leap (yr, mo, dy,hr, mn, sc)
if (n == 1) then
call print_date (clim_type%time_slice(1), &
str= 'for file ' // trim(file_name) // ', the &
&first time slice is mapped to :')
str= 'for file ' // trim(file_name) // ', the first time slice is mapped to :')
endif
if (n == ntime) then
call print_date (clim_type%time_slice(ntime), &
str= 'for file ' // trim(file_name) // ', the &
&last time slice is mapped to:')
str= 'for file ' // trim(file_name) // ', the last time slice is mapped to:')
endif

!---------------------------------------------------------------------
Expand All @@ -540,7 +581,7 @@ if(dimension_exists(fileobj, "time")) then
else
allocate(time_in(1), clim_type%time_slice(1))
allocate(clim_type%clim_times(1,1))
time_in = 0.0_lkind
time_in = 0.0_r8_kind
clim_type%time_slice = set_time(0,0) + base_time
clim_type%clim_times(1,1) = set_time(0,0) + base_time
endif
Expand Down
17 changes: 8 additions & 9 deletions interpolator/interpolator.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,9 +63,12 @@ module interpolator_mod
use time_manager_mod, only : time_type, &
set_time, &
set_date, &
get_date, &
time_type_to_real, &
days_in_year, &
get_calendar_type, &
leap_year, &
JULIAN, NOLEAP, &
get_date, &
get_date_julian, set_date_no_leap, &
set_date_julian, get_date_no_leap, &
print_date, &
Expand Down Expand Up @@ -435,8 +438,8 @@ subroutine interpolate_type_eq (Out, In)

Out%interph = In%interph
if (allocated(In%time_slice)) Out%time_slice = In%time_slice
Out%file_name = In%file_name
Out%time_flag = In%time_flag
Out%file_name = In%file_name
Out%time_flag = In%time_flag
Out%level_type = In%level_type
Out%is = In%is
Out%ie = In%ie
Expand Down Expand Up @@ -708,18 +711,14 @@ subroutine interpolator_end(clim_type)
deallocate(clim_type%r8_type%nmon_pyear)
end if
endif
if (allocated(clim_type%indexm)) deallocate(clim_type%indexm)
Copy link
Contributor Author

Choose a reason for hiding this comment

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

This set of deallocation was specified twice in interpolator_end. Removed it to make FMS more perfect.

if (allocated(clim_type%indexp)) deallocate(clim_type%indexp)
if (allocated(clim_type%climatology)) deallocate(clim_type%climatology)
if (allocated(clim_type%clim_times)) deallocate(clim_type%clim_times)

clim_type%r4_type%is_allocated=.false.
clim_type%r8_type%is_allocated=.false.

!! RSH mod
if( .not. (clim_type%TIME_FLAG .eq. LINEAR .and. &
if( .not.(clim_type%TIME_FLAG .eq. LINEAR .and. read_all_on_init) &
.and. (clim_type%TIME_FLAG.ne.NOTIME) ) then
! read_all_on_init)) .or. clim_type%TIME_FLAG .eq. BILINEAR ) then
read_all_on_init) ) then
call close_file(clim_type%fileobj)
endif

Expand Down
Loading
Loading