diff --git a/data_override/data_override.F90 b/data_override/data_override.F90 index da81222702..aeb9054157 100644 --- a/data_override/data_override.F90 +++ b/data_override/data_override.F90 @@ -1,3 +1,4 @@ + !*********************************************************************** !* GNU Lesser General Public License !* @@ -45,8 +46,10 @@ module data_override_mod assignment(=) use time_interp_external2_mod, only:time_interp_external_init, & time_interp_external, & + time_interp_external_bridge, & init_external_field, & get_external_field_size, & + get_time_axis, & set_override_region, & reset_src_data_region, & NO_REGION, INSIDE_REGION, OUTSIDE_REGION, & @@ -56,7 +59,7 @@ module data_override_mod use mpp_domains_mod, only : domain2d, mpp_get_compute_domain, NULL_DOMAIN2D,operator(.NE.),operator(.EQ.) use mpp_domains_mod, only : mpp_get_global_domain, mpp_get_data_domain use mpp_domains_mod, only : domainUG, mpp_pass_SG_to_UG, mpp_get_UG_SG_domain, NULL_DOMAINUG -use time_manager_mod, only: time_type +use time_manager_mod, only: time_type, OPERATOR(>), OPERATOR(<) use fms2_io_mod, only : FmsNetcdfFile_t, open_file, close_file, & read_data, fms2_io_init, variable_exists, & get_mosaic_tile_file @@ -75,10 +78,16 @@ module data_override_mod character(len=128) :: fieldname_code !< fieldname used in user's code (model) character(len=128) :: fieldname_file !< fieldname used in the netcdf data file character(len=512) :: file_name !< name of netCDF data file + character(len=512) :: prev_file_name !< name of netCDF data file for previous segment + character(len=512) :: next_file_name !< name of netCDF data file for next segment character(len=128) :: interpol_method !< interpolation method (default "bilinear") real :: factor !< For unit conversion, default=1, see OVERVIEW above real :: lon_start, lon_end, lat_start, lat_end integer :: region_type + logical :: multifile = .false. + type(time_type), dimension(:), pointer :: time_records => NULL() + type(time_type), dimension(:), pointer :: time_prev_records => NULL() + type(time_type), dimension(:), pointer :: time_next_records => NULL() end type data_type !> Private type for holding various data fields for performing data overrides @@ -87,6 +96,8 @@ module data_override_mod character(len=3) :: gridname character(len=128) :: fieldname integer :: t_index !< index for time interp + integer :: pt_index !< previous index for time interp + integer :: nt_index !< next index for time interp type(horiz_interp_type), allocatable :: horz_interp(:) !< index for horizontal spatial interp integer :: dims(4) !< dimensions(x,y,z,t) of the field in filename integer :: comp_domain(4) !< istart,iend,jstart,jend for compute domain @@ -232,6 +243,9 @@ subroutine data_override_init(Atm_domain_in, Ocean_domain_in, Ice_domain_in, Lan default_table%file_name = 'none' default_table%factor = 1. default_table%interpol_method = 'bilinear' + default_table%multifile = .false. + default_table%prev_file_name = 'none' + default_table%next_file_name = 'none' #ifdef use_yaml call read_table_yaml(data_table) @@ -357,6 +371,7 @@ subroutine read_table(data_table) integer :: iunit integer :: io_status + integer :: index_1col, index_2col character(len=256) :: record type(data_type) :: data_entry @@ -401,6 +416,20 @@ subroutine read_table(data_table) read(record,*,err=99) data_entry%gridname, data_entry%fieldname_code, data_entry%fieldname_file, & data_entry%file_name, data_entry%interpol_method, data_entry%factor, region, & & region_type + + if (index(data_entry%file_name, ":") .ne. 0) then + data_entry%multifile = .true. + index_1col = index(data_entry%file_name, ":") + index_2col = index(data_entry%file_name(index_1col+1:), ":") + data_entry%prev_file_name = data_entry%file_name(1:index_1col-1) + data_entry%next_file_name = data_entry%file_name(index_1col+index_2col+1:) + ! once previous/next files are filled in, overwrite current + data_entry%file_name = data_entry%file_name(index_1col+1:index_1col+index_2col-1) + else + data_entry%multifile = .false. + data_entry%prev_file_name = "" + data_entry%next_file_name = "" + endif if (data_entry%interpol_method == 'default') then data_entry%interpol_method = default_table%interpol_method endif @@ -443,6 +472,20 @@ subroutine read_table(data_table) ntable_lima = ntable_lima + 1 read(record,*,err=99) data_entry%gridname, data_entry%fieldname_code, data_entry%fieldname_file, & data_entry%file_name, ongrid, data_entry%factor + if (index(data_entry%file_name, ":") .ne. 0) then + data_entry%multifile = .true. + index_1col = index(data_entry%file_name, ":") + index_2col = index(data_entry%file_name(index_1col+1:), ":") + if (index_2col .eq. 0) call mpp_error(FATAL, "data_override: must provide either 1 or 3 forcing files") + data_entry%prev_file_name = data_entry%file_name(1:index_1col-1) + data_entry%next_file_name = data_entry%file_name(index_1col+index_2col+1:) + ! once previous/next files are filled in, overwrite current + data_entry%file_name = data_entry%file_name(index_1col+1:index_1col+index_2col-1) + else + data_entry%multifile = .false. + data_entry%prev_file_name = "" + data_entry%next_file_name = "" + endif if(ongrid) then data_entry%interpol_method = 'none' else @@ -457,6 +500,20 @@ subroutine read_table(data_table) ntable_new=ntable_new+1 read(record,*,err=99) data_entry%gridname, data_entry%fieldname_code, data_entry%fieldname_file, & data_entry%file_name, data_entry%interpol_method, data_entry%factor + if (index(data_entry%file_name, ":") .ne. 0) then + index_1col = index(data_entry%file_name, ":") + index_2col = index(data_entry%file_name(index_1col+1:), ":") + data_entry%multifile = .true. + if (index_2col .eq. 0) call mpp_error(FATAL, "data_override: must provide either 1 or 3 forcing files") + data_entry%prev_file_name = data_entry%file_name(1:index_1col-1) + data_entry%next_file_name = data_entry%file_name(index_1col+index_2col+1:) + ! once previous/next files are filled in, overwrite current + data_entry%file_name = data_entry%file_name(index_1col+1:index_1col+index_2col-1) + else + data_entry%multifile = .false. + data_entry%prev_file_name = "" + data_entry%next_file_name = "" + endif if (data_entry%interpol_method == 'default') then data_entry%interpol_method = default_table%interpol_method endif @@ -691,11 +748,15 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde character(len=512) :: filename !< file containing source data character(len=512) :: filename2 !< file containing source data + character(len=512) :: prevfilename, prevfilename2 !< file containing source data for previous segment + character(len=512) :: nextfilename, nextfilename2 !< file containing source data for next segment character(len=128) :: fieldname !< fieldname used in the data file integer :: i,j - integer :: dims(4) + integer :: dims(4), prev_dims(4), next_dims(4) integer :: index1 !< field index in data_table integer :: id_time !< index for time interp in override array + integer :: id_time_prev=-1 !< index for time interp bridge in override array + integer :: id_time_next=-1 !< index for time interp bridge in override array integer :: axis_sizes(4) character(len=32) :: axis_names(4) real, dimension(:,:), pointer :: lon_local =>NULL() !< of output (target) grid cells @@ -705,6 +766,7 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde logical :: data_file_is_2D = .false. !< data in netCDF file is 2D logical :: ongrid, use_comp_domain type(domain2D) :: domain + type(time_type) :: first_record, last_record integer :: curr_position !< position of the field currently processed in override_array real :: factor integer, dimension(4) :: comp_domain = 0 !< istart,iend,jstart,jend for compute domain @@ -717,6 +779,7 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde real :: lat_min, lat_max integer :: is_src, ie_src, js_src, je_src logical :: exists + logical :: multifile type(FmsNetcdfFile_t) :: fileobj integer :: startingi !< Starting x index for the compute domain relative to the input buffer integer :: endingi !< Ending x index for the compute domain relative to the input buffer @@ -750,6 +813,7 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde fieldname = data_table(index1)%fieldname_file ! fieldname in netCDF data file factor = data_table(index1)%factor + multifile = data_table(index1)%multifile if(fieldname == "") then data = factor @@ -758,6 +822,8 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde else filename = data_table(index1)%file_name if (filename == "") call mpp_error(FATAL,'data_override: filename not given in data_table') + if (multifile) prevfilename = data_table(index1)%prev_file_name + if (multifile) nextfilename = data_table(index1)%next_file_name endif ongrid = (data_table(index1)%interpol_method == 'none') @@ -832,22 +898,82 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde call get_mosaic_tile_file(filename,filename2,.false.,domain) filename = filename2 endif + if (multifile) then + inquire(file=trim(prevfilename),EXIST=exists) + if (.not. exists) then + call get_mosaic_tile_file(prevfilename,prevfilename2,.false.,domain) + prevfilename = prevfilename2 + endif + inquire(file=trim(nextfilename),EXIST=exists) + if (.not. exists) then + call get_mosaic_tile_file(nextfilename,nextfilename2,.false.,domain) + nextfilename = nextfilename2 + endif + endif !--- we always only pass data on compute domain id_time = init_external_field(filename,fieldname,domain=domain,verbose=.false., & use_comp_domain=use_comp_domain, nwindows=nwindows, ongrid=ongrid) + + if (multifile) then + id_time_prev = -1 + if (prevfilename /= 'none') then + id_time_prev = init_external_field(prevfilename,fieldname,domain=domain, & + verbose=.false.,use_comp_domain=use_comp_domain, & + nwindows = nwindows, ongrid=ongrid) + prev_dims = get_external_field_size(id_time_prev) + allocate(data_table(index1)%time_prev_records(prev_dims(4))) + call get_time_axis(id_time_prev,data_table(index1)%time_prev_records) + endif + id_time_next = -1 + if (nextfilename /= 'none') then + id_time_next = init_external_field(nextfilename,fieldname,domain=domain, & + verbose=.false.,use_comp_domain=use_comp_domain, & + nwindows = nwindows, ongrid=ongrid) + next_dims = get_external_field_size(id_time_next) + allocate(data_table(index1)%time_next_records(next_dims(4))) + call get_time_axis(id_time_next,data_table(index1)%time_next_records) + endif + endif + dims = get_external_field_size(id_time) override_array(curr_position)%dims = dims if(id_time<0) call mpp_error(FATAL,'data_override:field not found in init_external_field 1') override_array(curr_position)%t_index = id_time + override_array(curr_position)%pt_index = id_time_prev + override_array(curr_position)%nt_index = id_time_next else !ongrid=false id_time = init_external_field(filename,fieldname,domain=domain, axis_names=axis_names,& axis_sizes=axis_sizes, verbose=.false.,override=.true.,use_comp_domain=use_comp_domain, & nwindows = nwindows) + + if (multifile) then + id_time_prev = -1 + if (prevfilename /= 'none') then + id_time_prev = init_external_field(prevfilename,fieldname,domain=domain, axis_names=axis_names,& + axis_sizes=axis_sizes, verbose=.false.,override=.true.,use_comp_domain=use_comp_domain, & + nwindows = nwindows) + prev_dims = get_external_field_size(id_time_prev) + allocate(data_table(index1)%time_prev_records(prev_dims(4))) + call get_time_axis(id_time_prev,data_table(index1)%time_prev_records) + endif + id_time_next = -1 + if (nextfilename /= 'none') then + id_time_next = init_external_field(nextfilename,fieldname,domain=domain, axis_names=axis_names,& + axis_sizes=axis_sizes, verbose=.false.,override=.true.,use_comp_domain=use_comp_domain, & + nwindows = nwindows) + next_dims = get_external_field_size(id_time_next) + allocate(data_table(index1)%time_next_records(next_dims(4))) + call get_time_axis(id_time_next,data_table(index1)%time_next_records) + endif + endif + dims = get_external_field_size(id_time) override_array(curr_position)%dims = dims if(id_time<0) call mpp_error(FATAL,'data_override:field not found in init_external_field 2') override_array(curr_position)%t_index = id_time + override_array(curr_position)%pt_index = id_time_prev + override_array(curr_position)%nt_index = id_time_next ! get lon and lat of the input (source) grid, assuming that axis%data contains ! lat and lon of the input grid (in degrees) @@ -915,6 +1041,10 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde override_array(curr_position)%js_src = js_src override_array(curr_position)%je_src = je_src call reset_src_data_region(id_time, is_src, ie_src, js_src, je_src) + if (multifile) then + call reset_src_data_region(id_time_prev, is_src, ie_src, js_src, je_src) + call reset_src_data_region(id_time_next, is_src, ie_src, js_src, je_src) + endif ! Find the index of lon_start, lon_end, lat_start and lat_end in the input grid (nearest points) if( data_table(index1)%region_type .NE. NO_REGION ) then @@ -941,6 +1071,10 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde jstart = jstart - js_src + 1 jend = jend - js_src + 1 call set_override_region(id_time, data_table(index1)%region_type, istart, iend, jstart, jend) + if (multifile) then + call set_override_region(id_time_prev, data_table(index1)%region_type, istart, iend, jstart, jend) + call set_override_region(id_time_next, data_table(index1)%region_type, istart, iend, jstart, jend) + endif deallocate(lon_tmp, lat_tmp) endif @@ -961,6 +1095,8 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde endif !9 Get id_time previously stored in override_array id_time = override_array(curr_position)%t_index + id_time_prev = override_array(curr_position)%pt_index + id_time_next = override_array(curr_position)%nt_index endif !$OMP END CRITICAL @@ -1033,6 +1169,14 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde if(dims(3) .NE. 1 .and. (size(data,3) .NE. dims(3))) & call mpp_error(FATAL, "data_override: dims(3) .NE. 1 and size(data,3) .NE. dims(3)") + + dims = get_external_field_size(id_time) + if (.not. associated(data_table(index1)%time_records)) allocate(data_table(index1)%time_records(dims(4))) + call get_time_axis(id_time,data_table(index1)%time_records) + + first_record = data_table(index1)%time_records(1) + last_record = data_table(index1)%time_records(dims(4)) + if(ongrid) then if (.not. use_comp_domain) then !< Determine the size of the halox and the part of `data` that is in the compute domain @@ -1047,27 +1191,143 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde !10 do time interp to get data in compute_domain if(data_file_is_2D) then if (use_comp_domain) then - call time_interp_external(id_time,time,data(:,:,1),verbose=.false., & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + + if ((timelast_record)) then + if (multifile) then ! bridging between files + if (timelast_record) then + ! next file must be init and time must be between last record of current file and + ! first record of next file + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + ! bridge with next file + call time_interp_external_bridge(id_time,id_time_next,time,data(:,:,1),verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + else ! first_record <= time <= last_record, contrary to first if block + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else ! multifile = false + call mpp_error(FATAL, 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else ! standard behavior + call time_interp_external(id_time,time,data(:,:,1),verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif + else !> If this in an ongrid case and you are not in the compute domain, send in `data` to be the correct !! size - call time_interp_external(id_time,time,data(startingi:endingi,startingj:endingj,1),verbose=.false., & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + + if ((timelast_record)) then + if (multifile) then ! bridging between files + if (timelast_record) then + ! next file must be init and time must be between last record of current file and + ! first record of next file + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + ! bridge with next file + call time_interp_external_bridge(id_time,id_time_next,time,& + data(startingi:endingi,startingj:endingj,1),verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + else ! first_record <= time <= last_record, contrary to first if block + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else ! multifile = false + call mpp_error(FATAL, 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else ! standard behavior + call time_interp_external(id_time,time,data(startingi:endingi,startingj:endingj,1),verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif ! end bridge + end if + data(:,:,1) = data(:,:,1)*factor do i = 2, size(data,3) data(:,:,i) = data(:,:,1) end do else if (use_comp_domain) then - call time_interp_external(id_time,time,data,verbose=.false., & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + + if ((timelast_record)) then + if (multifile) then ! bridging between files, see previous blocks for comments + if (timelast_record) then + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + call time_interp_external_bridge(id_time,id_time_next,time,data,verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + else + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else + call mpp_error(FATAL, 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else + call time_interp_external(id_time,time,data,verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif + else !> If this in an ongrid case and you are not in the compute domain, send in `data` to be the correct !! size - call time_interp_external(id_time,time,data(startingi:endingi,startingj:endingj,:),verbose=.false., & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + if ((timelast_record)) then + if (multifile) then ! bridging between files, see previous blocks for comments + if (timelast_record) then + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + call time_interp_external_bridge(id_time,id_time_next,time,& + data(startingi:endingi,startingj:endingj,:),verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + else + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else + call mpp_error(FATAL, 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else + call time_interp_external(id_time,time,data(startingi:endingi,startingj:endingj,:),verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif + end if data = data*factor endif @@ -1075,9 +1335,37 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde ! do time interp to get global data if(data_file_is_2D) then if( data_table(index1)%region_type == NO_REGION ) then - call time_interp_external(id_time,time,data(:,:,1),verbose=.false., & - horz_interp=override_array(curr_position)%horz_interp(window_id), & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + + if ((timelast_record)) then + if (multifile) then ! bridging between files, see previous blocks for comments + if (timelast_record) then + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + call time_interp_external_bridge(id_time,id_time_next,time,data(:,:,1),verbose=.false., & + horz_interp=override_array(curr_position)%horz_interp(window_id), & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + else + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else + call mpp_error(FATAL, & + 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else + call time_interp_external(id_time,time,data(:,:,1),verbose=.false., & + horz_interp=override_array(curr_position)%horz_interp(window_id), & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif + data(:,:,1) = data(:,:,1)*factor do i = 2, size(data,3) data(:,:,i) = data(:,:,1) @@ -1085,10 +1373,40 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde else allocate(mask_out(size(data,1), size(data,2),1)) mask_out = .false. - call time_interp_external(id_time,time,data(:,:,1),verbose=.false., & - horz_interp=override_array(curr_position)%horz_interp(window_id), & - mask_out =mask_out(:,:,1), & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + + if ((timelast_record)) then + if (multifile) then ! bridging between files, see previous blocks for comments + if (timelast_record) then + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + call time_interp_external_bridge(id_time,id_time_next,time,data(:,:,1),verbose=.false., & + horz_interp=override_array(curr_position)%horz_interp(window_id), & + mask_out =mask_out(:,:,1), & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + else + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else + call mpp_error(FATAL,& + 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else + call time_interp_external(id_time,time,data(:,:,1),verbose=.false., & + horz_interp=override_array(curr_position)%horz_interp(window_id), & + mask_out =mask_out(:,:,1), & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif + where(mask_out(:,:,1)) data(:,:,1) = data(:,:,1)*factor end where @@ -1101,17 +1419,75 @@ subroutine data_override_3d(gridname,fieldname_code,data,time,override,data_inde endif else if( data_table(index1)%region_type == NO_REGION ) then - call time_interp_external(id_time,time,data,verbose=.false., & - horz_interp=override_array(curr_position)%horz_interp(window_id), & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + + if ((timelast_record)) then + if (multifile) then ! bridging between files, see previous blocks for comments + if (timelast_record) then + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + call time_interp_external_bridge(id_time,id_time_next,time,data,verbose=.false., & + horz_interp=override_array(curr_position)%horz_interp(window_id), & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + else + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else + call mpp_error(FATAL,& + 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else + call time_interp_external(id_time,time,data,verbose=.false., & + horz_interp=override_array(curr_position)%horz_interp(window_id), & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif + data = data*factor else allocate(mask_out(size(data,1), size(data,2), size(data,3)) ) mask_out = .false. - call time_interp_external(id_time,time,data,verbose=.false., & - horz_interp=override_array(curr_position)%horz_interp(window_id), & - mask_out =mask_out, & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + + if ((timelast_record)) then + if (multifile) then ! bridging between files, see previous blocks for comments + if (timelast_record) then + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + call time_interp_external_bridge(id_time,id_time_next,time,data,verbose=.false., & + horz_interp=override_array(curr_position)%horz_interp(window_id), & + mask_out =mask_out, & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + else + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else + call mpp_error(FATAL,& + 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else + call time_interp_external(id_time,time,data,verbose=.false., & + horz_interp=override_array(curr_position)%horz_interp(window_id), & + mask_out =mask_out, & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif + where(mask_out) data = data*factor end where @@ -1135,13 +1511,20 @@ subroutine data_override_0d(gridname,fieldname_code,data,time,override,data_inde real, intent(out) :: data !< output data array returned by this call integer, intent(in), optional :: data_index + type(time_type) :: first_record, last_record character(len=512) :: filename !< file containing source data + character(len=512) :: prevfilename !< file containing source data for bridge + character(len=512) :: nextfilename !< file containing source data for bridge character(len=128) :: fieldname !< fieldname used in the data file integer :: index1 !< field index in data_table + integer :: dims(4), prev_dims(4), next_dims(4) integer :: id_time !< index for time interp in override array + integer :: id_time_prev=-1 !index for time interp in override array + integer :: id_time_next=-1 !index for time interp in override array integer :: curr_position !< position of the field currently processed in override_array integer :: i real :: factor + logical :: multifile if(.not.module_is_initialized) & call mpp_error(FATAL,'Error: need to call data_override_init first') @@ -1167,6 +1550,7 @@ subroutine data_override_0d(gridname,fieldname_code,data,time,override,data_inde fieldname = data_table(index1)%fieldname_file ! fieldname in netCDF data file factor = data_table(index1)%factor + multifile = data_table(index1)%multifile if(fieldname == "") then data = factor @@ -1175,6 +1559,8 @@ subroutine data_override_0d(gridname,fieldname_code,data,time,override,data_inde else filename = data_table(index1)%file_name if (filename == "") call mpp_error(FATAL,'data_override: filename not given in data_table') + if (multifile) prevfilename = data_table(index1)%prev_file_name + if (multifile) nextfilename = data_table(index1)%next_file_name endif !3 Check if fieldname has been previously processed @@ -1203,8 +1589,54 @@ subroutine data_override_0d(gridname,fieldname_code,data,time,override,data_inde id_time = override_array(curr_position)%t_index endif !if curr_position < 0 + + if (multifile) then + id_time_prev = -1 + if (prevfilename /= 'none') then + id_time_prev = init_external_field(prevfilename,fieldname,verbose=.false.) + prev_dims = get_external_field_size(id_time_prev) + allocate(data_table(index1)%time_prev_records(prev_dims(4))) + call get_time_axis(id_time_prev,data_table(index1)%time_prev_records) + endif + id_time_next = -1 + if (nextfilename /= 'none') then + id_time_next = init_external_field(nextfilename,fieldname,verbose=.false.) + next_dims = get_external_field_size(id_time_next) + allocate(data_table(index1)%time_next_records(next_dims(4))) + call get_time_axis(id_time_next,data_table(index1)%time_next_records) + endif + endif + + !10 do time interp to get data in compute_domain - call time_interp_external(id_time, time, data, verbose=.false.) + + first_record = data_table(index1)%time_records(1) + last_record = data_table(index1)%time_records(dims(4)) + + if ((timelast_record)) then + if (multifile) then ! bridging between files + if (timelast_record) then + if (id_time_next<0) call mpp_error(FATAL,'data_override:next file needed with multifile') + if (time>data_table(index1)%time_next_records(1)) call mpp_error(FATAL, & + 'data_override: time_interp_external_bridge should only be called to bridge with next file') + call time_interp_external_bridge(id_time, id_time_next,time,data,verbose=.false.) + else + call mpp_error(FATAL, 'data_override: this should never happen') + endif + else + call mpp_error(FATAL, 'data_override: current time outside bounds, maybe add previous/next files to data_table') + endif + else + call time_interp_external(id_time,time,data,verbose=.false.) + endif + + data = data*factor !$OMP END SINGLE diff --git a/time_interp/time_interp_external2.F90 b/time_interp/time_interp_external2.F90 index fbe2f9e6f1..db869b0e4e 100644 --- a/time_interp/time_interp_external2.F90 +++ b/time_interp/time_interp_external2.F90 @@ -47,7 +47,7 @@ module time_interp_external2_mod use mpp_mod, only : mpp_error,FATAL,WARNING,mpp_pe, stdout, stdlog, NOTE use mpp_mod, only : input_nml_file, mpp_npes, mpp_root_pe, mpp_broadcast, mpp_get_current_pelist use time_manager_mod, only : time_type, get_date, set_date, operator ( >= ) , operator ( + ) , days_in_month, & - operator( - ), operator ( / ) , days_in_year, increment_time, & + operator( - ), operator ( / ) , operator ( // ), days_in_year, increment_time, & set_time, get_time, operator( > ), get_calendar_type, NO_CALENDAR use get_cal_time_mod, only : get_cal_time use mpp_domains_mod, only : domain2d, mpp_get_compute_domain, mpp_get_data_domain, & @@ -85,6 +85,7 @@ module time_interp_external2_mod time_interp_external_exit, get_external_field_size, get_time_axis, get_external_field_missing public set_override_region, reset_src_data_region public get_external_fileobj + public time_interp_external_bridge private find_buf_index,& set_time_modulo @@ -147,6 +148,12 @@ module time_interp_external2_mod module procedure time_interp_external_3d end interface + interface time_interp_external_bridge + module procedure time_interp_external_bridge_0d + module procedure time_interp_external_bridge_2d + module procedure time_interp_external_bridge_3d + end interface + !> @addtogroup time_interp_external2_mod !> @{ @@ -1009,6 +1016,259 @@ subroutine time_interp_external_0d(index, time, data, verbose) end subroutine time_interp_external_0d + + subroutine time_interp_external_bridge_2d(index1, index2, time, data_in, interp, verbose,horz_interp, mask_out, & + is_in, ie_in, js_in, je_in, window_id) + + integer, intent(in) :: index1, index2 + type(time_type), intent(in) :: time + real, dimension(:,:), intent(inout) :: data_in + integer, intent(in), optional :: interp + logical, intent(in), optional :: verbose + type(horiz_interp_type),intent(in), optional :: horz_interp + logical, dimension(:,:), intent(out), optional :: mask_out ! set to true where output data is valid + integer, intent(in), optional :: is_in, ie_in, js_in, je_in + integer, intent(in), optional :: window_id + + real , dimension(size(data_in,1), size(data_in,2), 1) :: data_out + logical, dimension(size(data_in,1), size(data_in,2), 1) :: mask3d + + data_out(:,:,1) = data_in(:,:) ! fill initial values for the portions of array that are not touched by 3d routine + call time_interp_external_bridge_3d(index1, index2, time, data_out, interp, verbose, horz_interp, mask3d, & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + data_in(:,:) = data_out(:,:,1) + if (PRESENT(mask_out)) mask_out(:,:) = mask3d(:,:,1) + + return + end subroutine time_interp_external_bridge_2d + + subroutine time_interp_external_bridge_3d(index1, index2, time, data, interp,verbose,horz_interp, mask_out, is_in, ie_in, js_in, je_in, window_id) + + integer, intent(in) :: index1, index2 + type(time_type), intent(in) :: time + real, dimension(:,:,:), intent(inout) :: data + integer, intent(in), optional :: interp + logical, intent(in), optional :: verbose + type(horiz_interp_type), intent(in), optional :: horz_interp + logical, dimension(:,:,:), intent(out), optional :: mask_out ! set to true where output data is valid + integer, intent(in), optional :: is_in, ie_in, js_in, je_in + integer, intent(in), optional :: window_id + + type(time_type) :: time1, time2 + integer :: dims1(4), dims2(4) + integer :: first_rec, last_rec + integer :: nx, ny, nz, interp_method, t1, t2 + integer :: i1, i2, isc, iec, jsc, jec, mod_time + integer :: isc1, iec1, jsc1, jec1, isc2, iec2, jsc2, jec2 + integer :: yy, mm, dd, hh, min, ss + character(len=256) :: err_msg, filename + + integer :: isw, iew, jsw, jew, nxw, nyw + ! these are boundaries of the updated portion of the "data" argument + ! they are calculated using sizes of the "data" and isc,iec,jsc,jsc + ! fileds from respective input field, to center the updated portion + ! in the output array + + real :: w1,w2 + logical :: verb + character(len=16) :: message1, message2 + + nx = size(data,1) + ny = size(data,2) + nz = size(data,3) + + interp_method = LINEAR_TIME_INTERP + if (PRESENT(interp)) interp_method = interp + verb=.false. + if (PRESENT(verbose)) verb=verbose + if (debug_this_module) verb = .true. + + if (index1 < 1) & + call mpp_error(FATAL,'invalid index in call to time_interp_external_bridge -- field was not initialized or failed to initialize') + if (index2 < 1) & + call mpp_error(FATAL,'invalid index in call to time_interp_external_bridge -- field was not initialized or failed to initialize') + + isc1=field(index1)%isc;iec1=field(index1)%iec + jsc1=field(index1)%jsc;jec1=field(index1)%jec + + isc2=field(index2)%isc;iec2=field(index2)%iec + jsc2=field(index2)%jsc;jec2=field(index2)%jec + + if ((isc1 /= isc2) .or. (iec1 /= iec2) .or. (jsc1 /= jsc2) .or. (jec1 /= jec2)) then + call mpp_error(FATAL, 'time_interp_external_bridge: dimensions must be the same in both index1 and index2 ') + endif + + isc=isc1 ; iec=iec1 ; jsc=jsc1 ; jec=jec1 + + if (trim(field(index2)%name) /= trim(field(index1)%name)) then + call mpp_error(FATAL, 'time_interp_external_bridge: cannot bridge between different fields.' // & + 'field1='//trim(field(index1)%name)//' field2='//trim(field(index2)%name)) + endif + + if ((field(index1)%numwindows == 1) .and. (field(index2)%numwindows == 1)) then + nxw = iec-isc+1 + nyw = jec-jsc+1 + else + if( .not. present(is_in) .or. .not. present(ie_in) .or. .not. present(js_in) .or. .not. present(je_in) ) then + call mpp_error(FATAL, 'time_interp_external_bridge: is_in, ie_in, js_in and je_in must be present ' // & + 'when numwindows > 1, field='//trim(field(index1)%name)) + endif + nxw = ie_in - is_in + 1 + nyw = je_in - js_in + 1 + isc = isc + is_in - 1 + iec = isc + ie_in - is_in + jsc = jsc + js_in - 1 + jec = jsc + je_in - js_in + endif + + isw = (nx-nxw)/2+1; iew = isw+nxw-1 + jsw = (ny-nyw)/2+1; jew = jsw+nyw-1 + + if (nx < nxw .or. ny < nyw .or. nz < field(index1)%siz(3)) then + write(message1,'(i6,2i5)') nx,ny,nz + call mpp_error(FATAL,'field '//trim(field(index1)%name)//' Array size mismatch in time_interp_external_bridge.'// & + ' Array "data" is too small. shape(data)='//message1) + endif + if (nx < nxw .or. ny < nyw .or. nz < field(index2)%siz(3)) then + write(message1,'(i6,2i5)') nx,ny,nz + call mpp_error(FATAL,'field '//trim(field(index2)%name)//' Array size mismatch in time_interp_external_bridge.'// & + ' Array "data" is too small. shape(data)='//message1) + endif + + if(PRESENT(mask_out)) then + if (size(mask_out,1) /= nx .or. size(mask_out,2) /= ny .or. size(mask_out,3) /= nz) then + write(message1,'(i6,2i5)') nx,ny,nz + write(message2,'(i6,2i5)') size(mask_out,1),size(mask_out,2),size(mask_out,3) + call mpp_error(FATAL,'field '//trim(field(index1)%name)//' array size mismatch in time_interp_external_bridge.'// & + ' Shape of array "mask_out" does not match that of array "data".'// & + ' shape(data)='//message1//' shape(mask_out)='//message2) + endif + endif + + if ((field(index1)%have_modulo_times) .or. (field(index2)%have_modulo_times)) then + call mpp_error(FATAL, 'time_interp_external_bridge: field '//trim(field(index1)%name)//' array cannot have modulo time') + endif + if ((field(index1)%modulo_time) .or. (field(index2)%modulo_time)) then + call mpp_error(FATAL, 'time_interp_external_bridge: field '//trim(field(index1)%name)//' array cannot have modulo time') + endif + + dims1 = get_external_field_size(index1) + dims2 = get_external_field_size(index2) + + t1 = dims1(4) + t2 = 1 + + time1 = field(index1)%time(t1) + time2 = field(index2)%time(t2) + w2= (time - time1) // (time2 - time1) + w1 = 1.0-w2 + + if (verb) then + call get_date(time,yy,mm,dd,hh,min,ss) + write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') & + 'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss + write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2 + endif + + call load_record(field(index1),t1,horz_interp, is_in, ie_in ,js_in, je_in, window_id) + call load_record(field(index2),t2,horz_interp, is_in, ie_in ,js_in, je_in, window_id) + i1 = find_buf_index(t1,field(index1)%ibuf) + i2 = find_buf_index(t2,field(index2)%ibuf) + if(i1<0.or.i2<0) & + call mpp_error(FATAL,'time_interp_external_bridge : records were not loaded correctly in memory') + + if (verb) then + write(outunit,*) 'ibuf= ',field(index2)%ibuf + write(outunit,*) 'i1,i2= ',i1, i2 + endif + + if( (field(index1)%region_type == NO_REGION) .and. (field(index2)%region_type == NO_REGION) ) then + where(field(index1)%mask(isc:iec,jsc:jec,:,i1).and.field(index2)%mask(isc:iec,jsc:jec,:,i2)) + data(isw:iew,jsw:jew,:) = field(index1)%data(isc:iec,jsc:jec,:,i1)*w1 + & + field(index2)%data(isc:iec,jsc:jec,:,i2)*w2 + elsewhere + data(isw:iew,jsw:jew,:) = field(index1)%missing + end where + else + where(field(index1)%mask(isc:iec,jsc:jec,:,i1).and.field(index2)%mask(isc:iec,jsc:jec,:,i2)) + data(isw:iew,jsw:jew,:) = field(index1)%data(isc:iec,jsc:jec,:,i1)*w1 + & + field(index2)%data(isc:iec,jsc:jec,:,i2)*w2 + end where + endif + + if(PRESENT(mask_out)) & + mask_out(isw:iew,jsw:jew,:) = & + field(index1)%mask(isc:iec,jsc:jec,:,i1).and.& + field(index2)%mask(isc:iec,jsc:jec,:,i2) + + end subroutine time_interp_external_bridge_3d + + subroutine time_interp_external_bridge_0d(index1, index2, time, data, verbose) + + integer, intent(in) :: index1, index2 + type(time_type), intent(in) :: time + real, intent(inout) :: data + logical, intent(in), optional :: verbose + + integer :: t1, t2 + integer :: i1, i2, mod_time + integer :: yy, mm, dd, hh, min, ss + character(len=256) :: err_msg, filename + type(time_type) :: time1, time2 + integer :: dims1(4), dims2(4) + + real :: w1,w2 + logical :: verb + + verb=.false. + if (PRESENT(verbose)) verb=verbose + if (debug_this_module) verb = .true. + + if (index1 < 0) & + call mpp_error(FATAL,'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize') + if (index2 < 0) & + call mpp_error(FATAL,'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize') + + if ((field(index1)%have_modulo_times) .or. (field(index2)%have_modulo_times)) then + call mpp_error(FATAL, 'time_interp_external_bridge: field '//trim(field(index1)%name)//' array cannot have modulo time') + endif + if ((field(index1)%modulo_time) .or. (field(index2)%modulo_time)) then + call mpp_error(FATAL, 'time_interp_external_bridge: field '//trim(field(index1)%name)//' array cannot have modulo time') + endif + + dims1 = get_external_field_size(index1) + dims2 = get_external_field_size(index2) + + t1 = dims1(4) + t2 = 1 + + time1 = field(index1)%time(t1) + time2 = field(index2)%time(t2) + w2= (time - time1) // (time2 - time1) + w1 = 1.0-w2 + + if (verb) then + call get_date(time,yy,mm,dd,hh,min,ss) + write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') & + 'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss + write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2 + endif + call load_record_0d(field(index2),t1) + call load_record_0d(field(index2),t2) + i1 = find_buf_index(t1,field(index2)%ibuf) + i2 = find_buf_index(t2,field(index2)%ibuf) + + if(i1<0.or.i2<0) & + call mpp_error(FATAL,'time_interp_external : records were not loaded correctly in memory') + data = field(index2)%data(1,1,1,i1)*w1 + field(index2)%data(1,1,1,i2)*w2 + if (verb) then + write(outunit,*) 'ibuf= ',field(index2)%ibuf + write(outunit,*) 'i1,i2= ',i1, i2 + endif + + end subroutine time_interp_external_bridge_0d + + subroutine set_time_modulo(Time) type(time_type), intent(inout), dimension(:) :: Time