diff --git a/interpolator/include/interpolator.inc b/interpolator/include/interpolator.inc index 754bae35cf..c774593561 100644 --- a/interpolator/include/interpolator.inc +++ b/interpolator/include/interpolator.inc @@ -145,7 +145,7 @@ 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 @@ -153,6 +153,15 @@ 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 @@ -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" @@ -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) @@ -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 & @@ -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 + 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 @@ -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 !-------------------------------------------------------------------- @@ -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), & + 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 !--------------------------------------------------------------------- @@ -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 diff --git a/interpolator/interpolator.F90 b/interpolator/interpolator.F90 index f598e2e56e..a00cf6b7c0 100644 --- a/interpolator/interpolator.F90 +++ b/interpolator/interpolator.F90 @@ -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, & @@ -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 @@ -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) -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 diff --git a/test_fms/interpolator/test_interpolator2.F90 b/test_fms/interpolator/test_interpolator2.F90 index 7a923a0a34..0fa62c3668 100644 --- a/test_fms/interpolator/test_interpolator2.F90 +++ b/test_fms/interpolator/test_interpolator2.F90 @@ -1,4 +1,4 @@ - !*********************************************************************** +!*********************************************************************** !* GNU Lesser General Public License !* !* This file is part of the GFDL Flexible Modeling System (FMS). @@ -31,8 +31,12 @@ program test_interpolator2 register_field, register_variable_attribute, & register_axis, & write_data, close_file, open_file - use mpp_mod, only: mpp_error, FATAL - use time_manager_mod, only: time_type, set_date, set_calendar_type, time_manager_init + use mpp_mod, only: mpp_error, FATAL, WARNING + use time_manager_mod, only: time_type, set_calendar_type, time_manager_init, & + get_date_no_leap, get_date_julian, set_date, set_time, & + set_date_no_leap, set_date_julian, operator(/), & + operator(+), operator(-), time_type_to_real, increment_time, & + leap_year, days_in_month, print_date, print_time use fms_mod, only: fms_init use constants_mod, only: PI use platform_mod, only: r4_kind, r8_kind @@ -41,15 +45,15 @@ program test_interpolator2 implicit none - character(100), parameter :: ncfile='immadeup.o3.climatology.nc' !< fake climatology file. - integer, parameter :: calendar_type=2 !< JULIAN calendar + character(100), parameter :: ncfile='immadeup.o3.climatology.nc' !< fake climatology file integer, parameter :: lkind=TEST_INTP_KIND_ - real(r8_kind), parameter :: tol=1.e-5_r8_kind !< the interpolation methods are not perfect. - !! Will not get perfectly agreeing answers + !> the interpolation methods are not perfect.Will not get perfectly agreeing answers + real(r8_kind) :: tol + integer :: calendar_type !> climatology related variables and arrays (made up data) - integer :: nlonlat !< number of latitude and longitudinal center coordinates - integer :: nlonlatb !< number of latitude and longitudinal boundary coordinates + integer :: nlonlat !< number of latitude and longitudinal center coordinates in file + integer :: nlonlatb !< number of latitude and longitudinal boundary coordinates in file integer :: ntime !< number of time slices integer :: npfull !< number of p levels integer :: nphalf !< number of half p levels @@ -57,7 +61,7 @@ program test_interpolator2 real(TEST_INTP_KIND_), allocatable :: lon(:) !< climatology coordinates real(TEST_INTP_KIND_), allocatable :: latb(:) !< climatology coordinates real(TEST_INTP_KIND_), allocatable :: lonb(:) !< climatology coordinates - real(TEST_INTP_KIND_), allocatable :: clim_time (:) !< climatology time + real(r8_kind), allocatable :: clim_time (:) !< climatology time real(TEST_INTP_KIND_), allocatable :: pfull(:) !< climatology p level real(TEST_INTP_KIND_), allocatable :: phalf(:) !< climatology p half level real(TEST_INTP_KIND_), allocatable :: ozone(:,:,:,:) !< climatology ozone data @@ -69,50 +73,64 @@ program test_interpolator2 real(TEST_INTP_KIND_), allocatable :: latb_mod(:,:) !< model coordinates real(TEST_INTP_KIND_), allocatable :: lonb_mod(:,:) !< model coordinates - type(interpolate_type) :: o3 !< recyclable interpolate_type - - call fms_init - call time_manager_init - call set_calendar_type(calendar_type) + !> array holding model times + type(time_type), allocatable :: model_time_julian(:), model_time_noleap(:) - !> set data - call set_write_data(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=5, npfull_in=3) + type(interpolate_type) :: o3 !< recyclable interpolate_type - !> test interpolator_init - write(*,*) '===== test_interpolator_init =====' - call test_interpolator_init(o3) + !> whether the file input is yearly, daily data + logical :: yearly, daily + logical :: noleap - !> test interpolator 2D-4D - write(*,*) '===== test_intepolator =======' - call test_interpolator(o3) + logical :: test_file_daily_julian=.true., test_file_daily_noleap=.false. + logical :: test_file_yearly_noleap=.false., test_file_yearly_julian=.false. + logical :: test_file_no_time=.false. + integer :: nml_unit_var=99 + character(*), parameter :: nml_file='test_interpolator.nml' + NAMELIST / test_interpolator_nml / test_file_daily_noleap, test_file_daily_julian, & + test_file_yearly_noleap, test_file_yearly_julian, test_file_no_time - !> test interpolate_type_eq - !! This test has been commented out and will be included - !! in the testing suite once fileobj cp is added into - !! test_interpolate_type_eq - !write(*,*) '===== test_interpolate_type_eq =====' - !call test_interpolate_type_eq() + if(lkind==r4_kind) tol=1.e-4_r8_kind + if(lkind==r8_kind) tol=1.e-6_r8_kind - !> test query_interpolator - write(*,*) '===== test_query_interpolator =====' - call test_query_interpolator() + open(unit=nml_unit_var, file=nml_file) + read(unit=nml_unit_var, nml=test_interpolator_nml) + close(nml_unit_var) - !> test interpolator end - write(*,*) '===== test_interpolator_end =====' - call test_interpolator_end(o3) + call fms_init + call time_manager_init + call write_header - !> deallocate all arrays used to write the .nc file and used for model coordinates - call deallocate_arrays() + !> set data + call set_parameters_wrapper + call set_and_write_data + + if(.not.test_file_no_time) then + !> test interpolator when model calendar is JULIAN + calendar_type=2 + call set_calendar_type(calendar_type) + call run_test_set + !--------------------------------------------------------- + !> test interpolator when model calendar is NOLEAP + calendar_type=4 + call set_calendar_type(calendar_type) + call run_test_set + end if !> test interpolator_no_time_axis - !! Write out new set of data that will have a time axis, but will have "0" time points - !! because that's how interpolator_init is set up. - call set_write_data(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=0, npfull_in=3) - call test_interpolator_init(o3) - write(*,*) '===== test_intepolator_no_time_axis =======' - call test_interpolator_no_time_axis(o3) + if(test_file_no_time) then + write(*,*) 'test_intepolator_no_time_axis' + calendar_type=2 !< still need to set model calendar + call set_calendar_type(calendar_type) + call test_interpolator_init(o3) + call test_interpolator_no_time_axis(o3) + call test_interpolator_end(o3) + end if contains + +#include "test_interpolator_write_climatology.inc" + !===============================================! subroutine test_interpolator_init(clim_type) @@ -129,7 +147,7 @@ subroutine test_interpolator_init(clim_type) end subroutine test_interpolator_init !===============================================! - subroutine test_interpolator(clim_type) + subroutine test_interpolator(clim_type, model_time) !> call the variants of interpolator (4D-2d) that interpolates data at a given time-point !! The tests here do not test the "no_axis" interpolator routines @@ -138,55 +156,59 @@ subroutine test_interpolator(clim_type) implicit none type(interpolate_type), intent(inout) :: clim_type - real(TEST_INTP_KIND_), dimension(nlonlat,nlonlat,npfull,1) :: interp_data !< last column, there is only one field - real(TEST_INTP_KIND_), dimension(nlonlat,nlonlat,nphalf) :: phalf - type(time_type) :: model_time + type(time_type), dimension(ntime), intent(in) :: model_time + type(time_type) :: test_time + real(TEST_INTP_KIND_), dimension(nlonlat_mod,nlonlat_mod,npfull,1) :: interp_data ! test interpolator_4D_r4/8 - call interpolator(clim_type, model_time, phalf, interp_data, 'ozone') + call interpolator(clim_type, test_time, phalf_in, interp_data, 'ozone') do i=1, npfull - do j=1, nlonlat - do k=1, nlonlat - call check_answers(interp_data(k,j,i,1), ozone(k,j,i,itime), tol, 'test interpolator_4D') + do j=1, nlonlat_mod + do k=1, nlonlat_mod + call check_answers(interp_data(k,j,i,1), answer, tol, 'test interpolator_4D') end do end do end do !> test interpolator_3_r4/8 - call interpolator(clim_type, model_time, phalf, interp_data(:,:,:,1), 'ozone') + call interpolator(clim_type, test_time, phalf_in, interp_data(:,:,:,1), 'ozone') do i=1, npfull - do j=1, nlonlat - do k=1, nlonlat - call check_answers(interp_data(k,j,i,1), ozone(k,j,i,itime), tol, 'test interpolator_3D') + do j=1, nlonlat_mod + do k=1, nlonlat_mod + call check_answers(interp_data(k,j,i,1), answer, tol, 'test interpolator_3D') end do end do end do !> test interpolator_2D_r4/8 - call interpolator(clim_type, model_time, interp_data(:,:,1,1), 'ozone') + call interpolator(clim_type, test_time, interp_data(:,:,1,1), 'ozone') do j=1, nlonlat_mod do k=1, nlonlat_mod - call check_answers(interp_data(k,j,1,1), ozone(k,j,1,itime), tol, 'test interpolator_2D') + call check_answers(interp_data(k,j,1,1), answer, tol, 'test interpolator_2D') end do end do !> Test obtain_interpolator_time_slices - call obtain_interpolator_time_slices(clim_type,model_time) - call interpolator(clim_type, model_time, interp_data(:,:,1,1), 'ozone') + call obtain_interpolator_time_slices(clim_type,test_time) + call interpolator(clim_type, test_time, interp_data(:,:,1,1), 'ozone') call unset_interpolator_time_flag(clim_type) do j=1, nlonlat_mod do k=1, nlonlat_mod - call check_answers(interp_data(k,j,1,1), ozone(k,j,1,itime), tol, 'test interpolator_2D') + call check_answers(interp_data(k,j,1,1), answer, tol, 'test interpolator_2D') end do end do @@ -214,18 +236,17 @@ subroutine test_interpolator_no_time_axis(clim_type) type(interpolate_type) :: clim_type - real(TEST_INTP_KIND_), dimension(nlonlat,nlonlat,nphalf-1,1) :: interp_data !< last column, there is only one field - real(TEST_INTP_KIND_), dimension(nlonlat,nlonlat,nphalf) :: phalf + real(TEST_INTP_KIND_), dimension(nlonlat,nlonlat,npfull,1) :: interp_data !< last column, there is only one field + real(TEST_INTP_KIND_), dimension(nlonlat,nlonlat,nphalf) :: phalf_in integer :: i, j, k - phalf(:,:,1)=0.0000_lkind - phalf(:,:,2)=0.0002_lkind - phalf(:,:,3)=0.0004_lkind - phalf(:,:,4)=0.0005_lkind + do i=1, nphalf + phalf_in(:,:,i)=phalf(i) + end do !> test interpolator_4D_no_time_axis_r4/8 - call interpolator(clim_type, phalf, interp_data, 'ozone') - do i=1, nphalf-1 + call interpolator(clim_type, phalf_in, interp_data, 'ozone') + do i=1, npfull do j=1, nlonlat do k=1, nlonlat call check_answers(interp_data(k,j,i,1), ozone(k,j,i,1), tol, 'test interpolator_4D_no_time_axis') @@ -234,8 +255,8 @@ subroutine test_interpolator_no_time_axis(clim_type) end do !> test interpolator_3D_no_time_axis_r4/8 - call interpolator(clim_type, phalf, interp_data(:,:,:,1), 'ozone') - do i=1, nphalf-1 + call interpolator(clim_type, phalf_in, interp_data(:,:,:,1), 'ozone') + do i=1, npfull do j=1, nlonlat do k=1, nlonlat call check_answers(interp_data(k,j,i,1), ozone(k,j,i,1), tol, 'test interpolator_3D_no_time_axis') @@ -256,14 +277,15 @@ end subroutine test_interpolator_no_time_axis subroutine test_interpolate_type_eq !> This subroutine tests interpolaote_type_eq (assignment = operator) - !! The success of "=" is insured by checking to see if interpolation with o3_copy succeds. + !! The success of "=" is insured by checking to see if interpolation with o3_copy succeeds. implicit none type(interpolate_type) :: o3_copy o3_copy = o3 - call test_interpolator(o3_copy) + if(calendar_type==2) call test_interpolator(o3_copy, model_time_julian) + if(calendar_type==4) call test_interpolator(o3_copy, model_time_noleap) end subroutine test_interpolate_type_eq !===============================================! @@ -290,6 +312,37 @@ subroutine test_query_interpolator end subroutine test_query_interpolator !===============================================! + subroutine run_test_set + + implicit none + integer :: i + + if(calendar_type==2) write(*,*) "** MODEL CALENDAR JULIAN ** MODEL CALENDAR JULIAN ** MODEL CALENDAR JULIAN **" + if(calendar_type==4) write(*,*) "** MODEL CALENDAR NOLEAP ** MODEL CALENDAR NOLEAP ** MODEL CALENDAR NOLEAP **" + + write(*,*) '1. interpolator_init' + call test_interpolator_init(o3) + + !> test interpolator 2D-4D + write(*,*) '2. interpolator4D-2D' + if(calendar_type==2) call test_interpolator(o3, model_time_julian) + if(calendar_type==4) call test_interpolator(o3, model_time_noleap) + + !> test interpolate_type_eq + !! This test has been commented out and will be included in the testing suite once fileobj cp is added into + write(*,*) '3. skipping interpolate_type_eq until bug in interpolator is fixed' + !call test_interpolate_type_eq() + + !> test query_interpolator + write(*,*) '4. query_interpolator' + call test_query_interpolator() + + !> test interpolator end + write(*,*) '5. interpolator_end' + call test_interpolator_end(o3) + + end subroutine run_test_set + !===============================================! subroutine check_answers(results, answers, tol, whoami) implicit none @@ -298,12 +351,52 @@ subroutine check_answers(results, answers, tol, whoami) character(*) :: whoami if (real(abs(results-answers),r8_kind).gt.tol) then + !if (results.ne.answers) then write(*,*) ' EXPECTED ', answers, ' but computed ', results call mpp_error(FATAL, trim(whoami)) end if end subroutine check_answers !===============================================! -#include "test_interpolator_write_climatology.inc" + subroutine write_header + + + if(test_file_daily_noleap) & + write(*,"(////10x,a,i0////)") & + " ** DAILY FILE CAL NOLEAP ** DAILY FILE CAL NOLEAP ** DAILY FILE CAL NOLEAP ** ", lkind + if(test_file_daily_julian) & + write(*,"(////10x,a,i0/////)") & + ' ** DAILY FILE CAL JULIAN ** DAILY FILE CAL JULIAN ** DAILY FILE CAL JULIAN ** ', lkind + if(test_file_yearly_noleap) & + write(*,"(////10x,a,i0/////)") & + ' ** YEARLY FILE CAL NOLEAP ** YEARLY FILE CAL NOLEAP ** YEARLY FILE CAL NOLEAP ** ', lkind + if(test_file_yearly_julian) & + write(*,"(////10x,a,i0/////)") & + ' ** YEARLY FILE CAL JULIAN ** YEARLY FILE CAL JULIAN ** YEARLY FILE CAL JULIAN ** ', lkind + if(test_file_no_time) & + write(*,"(////10x,a/////)") " ** NO TIME AXIS ** NO TIME AXIS ** NO TIME AXIS **" + + end subroutine write_header + !===============================================! + subroutine set_parameters_wrapper + + if(test_file_daily_noleap) then + call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=240, npfull_in=3, & + daily_in=.true., yearly_in=.false., noleap_in=.true.) + else if(test_file_daily_julian) then + call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=240, npfull_in=3, & + daily_in=.true., yearly_in=.false., noleap_in=.false.) + else if(test_file_yearly_noleap) then + call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=240, npfull_in=3, & + daily_in=.false., yearly_in=.true., noleap_in=.true.) + else if(test_file_yearly_julian) then + call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=240, npfull_in=3, & + daily_in=.false., yearly_in=.true., noleap_in=.false.) + else if(test_file_no_time) then + call set_parameters(nlonlat_in=10, nlonlat_mod_in=10, ntime_in=0, npfull_in=3, & + daily_in=.true., yearly_in=.false., noleap_in=.false.) + end if + end subroutine set_parameters_wrapper + !===============================================! end program test_interpolator2 diff --git a/test_fms/interpolator/test_interpolator2.sh b/test_fms/interpolator/test_interpolator2.sh index b1a6110027..7248b92c22 100755 --- a/test_fms/interpolator/test_interpolator2.sh +++ b/test_fms/interpolator/test_interpolator2.sh @@ -53,13 +53,87 @@ cat <<_EOF > input.nml _EOF # Run test -test_expect_success "test interpolator" 'mpirun -n 2 ./test_interpolator' +test_expect_success "test interpolator" 'mpirun -n 1 ./test_interpolator' -#Run the second set of interpolator tests + +#Run the daily interpolator tests when the file calendar is in units of days and calendar type is NOLEAP +cat < test_interpolator.nml +&test_interpolator_nml +test_file_daily_noleap=.false. +test_file_daily_julian=.true. +test_file_yearly_noleap=.false. +test_file_yearly_julian=.false. +test_file_no_time=.false. +/ +EOF +mkdir -p INPUT +test_expect_success "test_interpolator2 file data daily julian r4 unit tests" 'mpirun -n 1 ./test_interpolator2_r4' +test_expect_success "test_interpolator2 file data daily julian r8 unit tests" 'mpirun -n 1 ./test_interpolator2_r8' +rm -rf INPUT *.nc test_interpolator.nml + + +#Run the daily interpolator tests when the file calendar is in units of days and calendar type is JULIAN +cat < test_interpolator.nml +&test_interpolator_nml +test_file_daily_noleap=.true. +test_file_daily_julian=.false. +test_file_yearly_noleap=.false. +test_file_yearly_julian=.false. +test_file_no_time=.false. +/ +EOF +mkdir -p INPUT +test_expect_success "test_interpolator2 file data daily noleap r4 unit tests" 'mpirun -n 1 ./test_interpolator2_r4' +test_expect_success "test_interpolator2 file data daily noleap r8 unit tests" 'mpirun -n 1 ./test_interpolator2_r8' +rm -rf INPUT *.nc test_interpolator.nml + + +#Run the yearly interpolator tests when the file calendar is in units of years and calendar type is NOLEAP +cat < test_interpolator.nml +&test_interpolator_nml +test_file_daily_noleap=.false. +test_file_daily_julian=.false. +test_file_yearly_noleap=.true. +test_file_yearly_julian=.false. +test_file_no_time=.false. +/ +EOF +mkdir -p INPUT +test_expect_success "test_interpolator2 file data yearly noleap r4 unit tests" 'mpirun -n 1 ./test_interpolator2_r4' +test_expect_success "test_interpolator2 file data yearly noleap r8 unit tests" 'mpirun -n 1 ./test_interpolator2_r8' +rm -rf INPUT *.nc test_interpolator.nml + + +#Run the yearly interpolator tests when the file calendar is in units of years and calendar type is JULIAN +cat < test_interpolator.nml +&test_interpolator_nml +test_file_daily_noleap=.false. +test_file_daily_julian=.false. +test_file_yearly_noleap=.false. +test_file_yearly_julian=.true. +test_file_no_time=.false. +/ +EOF +mkdir -p INPUT +test_expect_success "test_interpolator2 file data yearly julian r4 unit tests" 'mpirun -n 1 ./test_interpolator2_r4' +test_expect_success "test_interpolator2 file data yearly julian r8 unit tests" 'mpirun -n 1 ./test_interpolator2_r8' +rm -rf INPUT *.nc test_interpolator.nml + + +#Run no_time_axis +cat < test_interpolator.nml +&test_interpolator_nml +test_file_daily_noleap=.false. +test_file_daily_julian=.false. +test_file_yearly_noleap=.false. +test_file_yearly_julian=.false. +test_file_no_time=.true. +/ +EOF mkdir -p INPUT -test_expect_success "test_interpolator2 r4 unit tests" 'mpirun -n 1 ./test_interpolator2_r4' -test_expect_success "test_interpolator2 r8 unit tests" 'mpirun -n 1 ./test_interpolator2_r8' +test_expect_success "test_interpolator2 file data no time axis r4 unit tests" 'mpirun -n 1 ./test_interpolator2_r4' +test_expect_success "test_interpolator2 file data no time axis r8 unit tests" 'mpirun -n 1 ./test_interpolator2_r8' +rm -rf INPUT *.nc test_interpolator.nml -rm -rf INPUT *.nc # remove any leftover io files to save space test_done diff --git a/test_fms/interpolator/test_interpolator_write_climatology.inc b/test_fms/interpolator/test_interpolator_write_climatology.inc index 3d6ce90758..c05c5deed5 100644 --- a/test_fms/interpolator/test_interpolator_write_climatology.inc +++ b/test_fms/interpolator/test_interpolator_write_climatology.inc @@ -24,9 +24,6 @@ subroutine write_climatology_file implicit none type(FmsNetcdfFile_t) :: fileobj - integer :: i,j,k,l, itime - integer :: stringlen - character(100) :: mystring ! write netcdf file if(open_file(fileobj, 'INPUT/'//trim(ncfile), 'overwrite')) then @@ -44,9 +41,10 @@ subroutine write_climatology_file !if(ntime /=0 ) then call register_field(fileobj, 'time', 'double', dimensions=(/'time'/)) call register_string_attribute(fileobj, 'time', 'axis', 'T') - call register_string_attribute(fileobj, 'time', 'calendar', 'noleap') - call register_string_attribute(fileobj, 'time', 'climatology', '1979-01-01 00:00:00, 1998-01-01 00:00:00') - call register_string_attribute(fileobj, 'time', 'units', 'days since 1849-01-01 00:00:00') + if(noleap) call register_string_attribute(fileobj, 'time', 'calendar', 'noleap') + if(.not.noleap) call register_string_attribute(fileobj, 'time', 'calendar', 'julian') + if(daily) call register_string_attribute(fileobj, 'time', 'units', 'days since 1849-01-01 01:01:01') + if(yearly) call register_string_attribute(fileobj, 'time', 'units', 'years since 1849-01-01 01:01:01') !end if !lon axes @@ -119,52 +117,31 @@ subroutine register_string_attribute(fileobj, variable, attribute, att_value) end subroutine register_string_attribute !===============================================! -subroutine allocate_arrays() - - implicit none - integer :: ntime2 - - ntime2=ntime - if(ntime==0) ntime2=1 - - allocate( lat(nlonlat), & - lon(nlonlat), & - latb(nlonlatb), & - lonb(nlonlatb), & - clim_time(ntime2),& - pfull(npfull), & - phalf(nphalf), & - ozone(nlonlat, nlonlat, npfull, ntime2) ) - - allocate( lat_mod(nlonlat_mod, nlonlat_mod), & - lon_mod(nlonlat_mod, nlonlat_mod), & - latb_mod(nlonlatb_mod, nlonlatb_mod), & - lonb_mod(nlonlatb_mod, nlonlatb_mod) ) - -end subroutine allocate_arrays -!===============================================! subroutine deallocate_arrays() implicit none - deallocate( lat, & - lon, & - latb, & - lonb, & - clim_time, & + deallocate( lat, & + lon, & + latb, & + lonb, & + clim_time, & + model_time_noleap, & + model_time_julian, & pfull, & phalf, & - lat_mod, & - lon_mod, & - latb_mod, & - lonb_mod, & + lat_mod, & + lon_mod, & + latb_mod, & + lonb_mod, & ozone ) end subroutine deallocate_arrays !===============================================! -subroutine set_write_data(nlonlat_in, nlonlat_mod_in, ntime_in, npfull_in) +subroutine set_parameters(nlonlat_in, nlonlat_mod_in, ntime_in, npfull_in, daily_in, yearly_in, noleap_in) implicit none integer, intent(in) :: nlonlat_in, nlonlat_mod_in, ntime_in, npfull_in + logical, intent(in) :: daily_in, yearly_in, noleap_in nlonlat = nlonlat_in !< number of latitude and longitudinal center coordinates nlonlatb = nlonlat_in+1 !< number of latitude and longitudinal boundary coordinates @@ -173,25 +150,117 @@ subroutine set_write_data(nlonlat_in, nlonlat_mod_in, ntime_in, npfull_in) ntime = ntime_in !< number of time slices npfull = npfull_in !< number of p levels nphalf = npfull_in+1 !< number of half p levels - call allocate_arrays() + + daily = daily_in + yearly = yearly_in + noleap = noleap_in + +end subroutine set_parameters +!===============================================! +subroutine set_and_write_data + + implicit none + if( ntime /= 0) call set_clim_time() call set_latlon_b() + call set_latlon_b_mod() call set_pfullhalf() call set_ozone() - call write_climatology_file() + call write_climatology_file - call set_latlon_b_mod() - -end subroutine set_write_data +end subroutine set_and_write_data !===============================================! subroutine set_clim_time() implicit none integer :: i - do i=1, ntime - clim_time(i)=2*i - end do + type(time_type) :: base_time + integer :: l, ii, yr, mo, dy, hr, mn, sc + integer :: ntime2 + + if(allocated(clim_time)) deallocate(clim_time) + if(allocated(model_time_julian)) deallocate(model_time_julian) + if(allocated(model_time_noleap)) deallocate(model_time_noleap) + allocate(clim_time(ntime)) + allocate(model_time_julian(ntime), model_time_noleap(ntime)) + + !write(*,*) ' -- SETTING TIME -- ' + + hr = 1 ; mn = 1 ; sc = 1 + + if(test_file_daily_noleap) then + ! base_time must match that in the time attribute + base_time=set_date_no_leap(1849,1,1,hr,mn,sc) + yr = 1849 ; mo = 1 ; dy = 15 + do i=1, ntime + yr = yr + 1 + mo = mo + 1 ; if( mo > 12 ) mo=1 + dy = 15 + model_time_julian(i)=set_date_julian(yr, mo, dy, hr, mn, sc) + model_time_noleap(i)=set_date_no_leap(yr, mo, dy, hr, mn, sc) + clim_time(i)=time_type_to_real(model_time_noleap(i)-base_time)/86400._r8_kind + !call print_time(model_time_noleap(i)) + end do + else if(test_file_daily_julian) then + ! base_time must match that in the time attribute + base_time=set_date_julian(1849,1,1,hr,mn,sc) + yr = 1849 ; mo = 1 ; dy = 15 + do i=1, ntime + yr = yr + 1 + mo = mo + 1 ; if( mo > 12 ) mo=1 + dy = 15 + model_time_julian(i)=set_date_julian(yr, mo, dy, hr, mn, sc) + model_time_noleap(i)=set_date_no_leap(yr, mo, dy, hr, mn, sc) + clim_time(i)=time_type_to_real(model_time_julian(i)-base_time)/86400._r8_kind + !call print_time(model_time_julian(i)) + end do + else if(test_file_yearly_noleap) then + ! base_time must match that in the time attribute + base_time=set_date_no_leap(1849,1,1,hr,mn,sc) + yr = 1849 ; mo = 1 ; dy = 15 + do i=1, ntime + yr = yr + 1 + mo = mo + 1 ; if( mo > 12 ) mo=1 + dy = 15 + model_time_julian(i)=set_date_julian(yr, mo, dy, hr, mn, sc) + model_time_noleap(i)=set_date_no_leap(yr, mo, dy, hr, mn, sc) + clim_time(i)=real(yr-1849,r8_kind) + do ii=1, mo-1 + select case(ii) + case(1,3,5,7,8,10,12) ; clim_time(i) = clim_time(i) + 31.0_r8_kind/365._r8_kind + case(2) ; clim_time(i) = clim_time(i) + 28._r8_kind/365._r8_kind + case(4,6,9,11) ; clim_time(i) = clim_time(i) + 30.0_r8_kind/365._r8_kind + end select + end do + clim_time(i) = clim_time(i) + real(dy-1,r8_kind)/365._r8_kind + !call print_time(model_time_noleap(i)) + end do + else if(test_file_yearly_julian) then + ! base_time must match that in the time attribute + base_time=set_date_julian(1849,1,1,hr,mn,sc) + yr = 1849 ; mo = 1 ; dy = 15 + call set_calendar_type(2) + do i=1, ntime + yr = yr + 1 + mo = mo + 1 ; if( mo > 12 ) mo=1 + dy = 15 + model_time_julian(i)=set_date_julian(yr,mo,dy,hr,mn,sc) + model_time_noleap(i)=set_date_no_leap(yr,mo,dy,hr,mn,sc) + clim_time(i)=real(yr-1849,lkind) + l=0 ; if(leap_year(model_time_julian(i))) l=l+1 + do ii=1, mo-1 + select case(ii) + case(1,3,5,7,8,10,12) ; clim_time(i) = clim_time(i) + 31.0_r8_kind/real(365+l,r8_kind) + case(2) ; clim_time(i) = clim_time(i) + real(28+l,r8_kind)/real(365+l,r8_kind) + case(4,6,9,11) ; clim_time(i) = clim_time(i) + 30.0_r8_kind/real(365+l,r8_kind) + end select + end do + clim_time(i) = clim_time(i) + real(dy-1,r8_kind)/real(365+l,r8_kind) + !call print_time(model_time_julian(i)) + end do + call set_calendar_type(0) + end if end subroutine set_clim_time !===============================================! @@ -200,6 +269,14 @@ subroutine set_latlon_b() implicit none integer :: i + if(allocated(lat)) deallocate(lat) + if(allocated(lon)) deallocate(lon) + if(allocated(latb)) deallocate(latb) + if(allocated(lonb)) deallocate(lonb) + + allocate(lat(nlonlat), lon(nlonlat)) + allocate(latb(nlonlat+1), lonb(nlonlat+1)) + do i=1, nlonlat lat(i)= real(2*i-1,TEST_INTP_KIND_) lon(i)= real(2*i-1,TEST_INTP_KIND_) @@ -222,20 +299,26 @@ subroutine set_latlon_b_mod implicit none integer :: i, j - !> the model grid and the climatology grid are the same - do i=1, nlonlat - lon_mod(:,i) = real(2*i-1,TEST_INTP_KIND_) - lat_mod(i,:) = real(2*i-1,TEST_INTP_KIND_) - end do + if(allocated(lat_mod)) deallocate(lat_mod) + if(allocated(lon_mod)) deallocate(lon_mod) + if(allocated(latb_mod)) deallocate(latb_mod) + if(allocated(lonb_mod)) deallocate(lonb_mod) - lonb_mod(1:nlonlatb,1)=1.0_lkind - latb_mod(1,1:nlonlatb)=1.0_lkind - do i=1, nlonlat-1 - lonb_mod(1:nlonlatb,i+1)=2*i - latb_mod(i+1,1:nlonlatb)=2*i + allocate(lat_mod(nlonlat_mod,nlonlat_mod), lon_mod(nlonlat_mod,nlonlat_mod)) + allocate(latb_mod(nlonlatb_mod,nlonlatb_mod), lonb_mod(nlonlatb_mod,nlonlatb_mod)) + + !> nlonlat = nlonlat_mod + ! the model coordinates are the same as the file coordinates + do i=1, nlonlat_mod + lat_mod(i,:) = lat(i) + lon_mod(:,i) = lon(i) + end do + !> nlonlatb_mod = nlonlatb + !! the model coordinates are the same as the file coordinates + do i=1, nlonlatb_mod + latb_mod(i,:) = latb(i) + lonb_mod(:,i) = lonb(i) end do - lonb_mod(1:nlonlatb,nlonlatb)=real(2*nlonlat-1,TEST_INTP_KIND_) - latb_mod(nlonlatb,1:nlonlatb)=real(2*nlonlat-1,TEST_INTP_KIND_) !> convert from degrees to radians lon_mod = lon_mod*real(PI,TEST_INTP_KIND_)/180.0_lkind @@ -251,6 +334,10 @@ subroutine set_pfullhalf() implicit none integer :: i + if(allocated(pfull)) deallocate(pfull) + if(allocated(phalf)) deallocate(phalf) + allocate(pfull(npfull), phalf(nphalf)) + do i=1, npfull pfull(i) = 0.0001_lkind * real(i-1,TEST_INTP_KIND_) end do @@ -271,12 +358,13 @@ subroutine set_ozone() ntime2=ntime if(ntime == 0 ) ntime2=1 - + if(allocated(ozone)) deallocate(ozone) + allocate(ozone(nlonlat, nlonlat, npfull, ntime2)) do i=1, ntime2 do j=1, npfull do k=1, nlonlat do l=1, nlonlat - ozone(l,k,j,i)=real(i,TEST_INTP_KIND_) + ozone(l,k,j,i)= real(i,TEST_INTP_KIND_) end do end do end do