diff --git a/data_override/README.MD b/data_override/README.MD index 8790d94bbf..f9e19464aa 100644 --- a/data_override/README.MD +++ b/data_override/README.MD @@ -24,6 +24,13 @@ If it is desired to interpolate the data to a region of the model grid. The foll - **lat_start:** The starting longitude in the same units as the grid data in the file - **lon_end:** The ending longitude in the same units as the grid data in the file +If it is desired to use multiple(3) input netcdf files instead of 1. The following **optional** keys are available. +- **is_multi_file:** Set to `True` is using the multi-file feature +- **prev_file_name:** The name of the first file in the set +- **next_file_name:** The name of the third file in the set + +Note that **file_name** must be the second file in the set. **prev_file_name** and/or **next_file_name** are required if **is_multi_file** is set to `True` + #### 2. How to use it? In order to use the yaml data format, [libyaml](https://github.com/yaml/libyaml) needs to be installed and linked with FMS. Additionally, FMS must be compiled with -Duse_yaml macro. If using autotools, you can add `--with-yaml`, which will add the macro for you and check that libyaml is linked correctly. ``` diff --git a/data_override/include/data_override.inc b/data_override/include/data_override.inc index 305e58ae6c..a4b5feec12 100644 --- a/data_override/include/data_override.inc +++ b/data_override/include/data_override.inc @@ -30,6 +30,7 @@ use horiz_interp_mod, only : horiz_interp_init, horiz_interp_new, horiz_interp_t assignment(=) use time_interp_external2_mod, only: time_interp_external_init, & time_interp_external, & + time_interp_external_bridge, get_time_axis, & init_external_field, & get_external_field_size, & set_override_region, & @@ -41,7 +42,7 @@ use axis_utils2_mod, only : nearest_index, axis_edges 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 @@ -64,6 +65,12 @@ type data_type real(FMS_DATA_OVERRIDE_KIND_) :: factor !< For unit conversion, default=1, see OVERVIEW above real(FMS_DATA_OVERRIDE_KIND_) :: lon_start, lon_end, lat_start, lat_end integer :: region_type + logical :: multifile = .false. + 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 + 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 @@ -72,6 +79,8 @@ type override_type 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 @@ -198,6 +207,9 @@ end if default_table%file_name = 'none' default_table%factor = 1._lkind default_table%interpol_method = 'bilinear' + default_table%multifile = .false. + default_table%prev_file_name = '' + default_table%next_file_name = '' #ifdef use_yaml if (use_data_table_yaml) then @@ -355,6 +367,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 @@ -399,6 +412,22 @@ 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:), ":") + if (index_2col .eq. 0) call mpp_error(FATAL, "data_override: when bridging over forcing files, " & + // "central forcing files must be preceded AND followed by the column (:) separator") + 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 @@ -441,6 +470,21 @@ 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: when bridging over forcing files, " & + // "central forcing files must be preceded AND followed by the column (:) separator") + 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 @@ -455,6 +499,21 @@ 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: when bridging over forcing files, " & + // "central forcing files must be preceded AND followed by the column (:) separator") + 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 @@ -517,6 +576,23 @@ subroutine read_table_yaml(data_table) call get_value_from_key(file_id, entry_id(i), "fieldname_file", data_table(i)%fieldname_file, & & is_optional=.true.) + data_table(i)%multifile = .false. + call get_value_from_key(file_id, entry_id(i), "is_multi_file", data_table(i)%multifile, & + & is_optional=.true.) + + if (data_table(i)%multifile) then + data_table(i)%prev_file_name = "" + data_table(i)%next_file_name = "" + call get_value_from_key(file_id, entry_id(i), "prev_file_name", data_table(i)%prev_file_name, & + & is_optional=.true.) + call get_value_from_key(file_id, entry_id(i), "next_file_name", data_table(i)%next_file_name, & + & is_optional=.true.) + if (trim(data_table(i)%prev_file_name) .eq. "" .and. trim(data_table(i)%next_file_name) .eq. "") & + call mpp_error(FATAL, "The prev_file_name and next_file_name must be present if is_multi_file. "//& + "Check your data_table.yaml entry for field:"//trim(data_table(i)%gridname)//":"//& + trim(data_table(i)%fieldname_code)) + endif + data_table(i)%file_name = "" call get_value_from_key(file_id, entry_id(i), "file_name", data_table(i)%file_name, & & is_optional=.true.) @@ -627,13 +703,23 @@ subroutine DATA_OVERRIDE_0D_(gridname,fieldname_code,data_out,time,override,data real(FMS_DATA_OVERRIDE_KIND_), intent(out) :: data_out !< output data array returned by this call integer, intent(in), optional :: data_index + type(time_type) :: first_record !< first record of "current" file + type(time_type) :: last_record !< last record of "current" file character(len=512) :: filename !< file containing source data + character(len=512) :: prevfilename !< file containing previous source data, when using multiple files + character(len=512) :: nextfilename !< file containing next source data, when using multiple files character(len=128) :: fieldname !< fieldname used in the data file integer :: index1 !< field index in data_table + integer :: dims(4) + integer :: prev_dims(4) !< dimensions of previous source data, when using multiple files + integer :: next_dims(4) !< dimensions of next source data, when using multiple files integer :: id_time !< index for time interp in override array + integer :: id_time_prev=-1 !< time index for previous file, when using multiple files + integer :: id_time_next=-1 !< time index for next file, when using multiple files integer :: curr_position !< position of the field currently processed in override_array integer :: i real(FMS_DATA_OVERRIDE_KIND_) :: factor + logical :: multifile !< use multiple consecutive files for override if(.not.module_is_initialized) & call mpp_error(FATAL,'Error: need to call data_override_init first') @@ -659,6 +745,7 @@ subroutine DATA_OVERRIDE_0D_(gridname,fieldname_code,data_out,time,override,data 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_out = factor @@ -667,6 +754,8 @@ subroutine DATA_OVERRIDE_0D_(gridname,fieldname_code,data_out,time,override,data 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 @@ -695,8 +784,71 @@ subroutine DATA_OVERRIDE_0D_(gridname,fieldname_code,data_out,time,override,data id_time = override_array(curr_position)%t_index endif !if curr_position < 0 + + ! if using consecutive files for data_override, get time axis for previous and next files + ! and check spatial dims for consistency + if_multi1: if (multifile) then + id_time_prev = -1 + if_prev1: if (trim(prevfilename) /= '') then + id_time_prev = init_external_field(prevfilename,fieldname,verbose=.false.) + dims = get_external_field_size(id_time) + prev_dims = get_external_field_size(id_time_prev) + ! check consistency of spatial dims + if ((prev_dims(1) .ne. dims(1)) .or. (prev_dims(2) .ne. dims(2)) .or. & + (prev_dims(3) .ne. dims(3))) then + call mpp_error(FATAL, 'data_override: dimensions mismatch between consecutive forcing files') + endif + allocate(data_table(index1)%time_prev_records(prev_dims(4))) + call get_time_axis(id_time_prev,data_table(index1)%time_prev_records) + endif if_prev1 + id_time_next = -1 + if_next1: if (trim(nextfilename) /= '') then + id_time_next = init_external_field(nextfilename,fieldname,verbose=.false.) + dims = get_external_field_size(id_time) + next_dims = get_external_field_size(id_time_next) + ! check consistency of spatial dims + if ((next_dims(1) .ne. dims(1)) .or. (next_dims(2) .ne. dims(2)) .or. & + (next_dims(3) .ne. dims(3))) then + call mpp_error(FATAL, 'data_override: dimensions mismatch between consecutive forcing files') + endif + allocate(data_table(index1)%time_next_records(next_dims(4))) + call get_time_axis(id_time_next,data_table(index1)%time_next_records) + endif if_next1 + endif if_multi1 + + !10 do time interp to get data in compute_domain - call time_interp_external(id_time, time, data_out, verbose=.false.) + + first_record = data_table(index1)%time_records(1) + last_record = data_table(index1)%time_records(dims(4)) + + ! if using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi2: if (multifile) then + if_time2: 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_out,verbose=.false.) + else ! first_record < time < last_record, do not use bridge + call time_interp_external(id_time,time,data_out,verbose=.false.) + endif if_time2 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,data_out,verbose=.false.) + endif if_multi2 + + data_out = data_out*factor !$OMP END SINGLE @@ -749,11 +901,19 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d character(len=512) :: filename !< file containing source data character(len=512) :: filename2 !< file containing source data + character(len=512) :: prevfilename !< file containing source data for previous file + character(len=512) :: prevfilename2 !< file containing source data for previous file + character(len=512) :: nextfilename !< file containing source data for next file + character(len=512) :: nextfilename2 !< file containing source data for next file character(len=128) :: fieldname !< fieldname used in the data file integer :: i,j integer :: dims(4) + integer :: prev_dims(4) !< dimensions of previous source data, when using multiple files + integer :: next_dims(4) !< dimensions of next source data, when using multiple files integer :: index1 !< field index in data_table integer :: id_time !< index for time interp in override array + integer :: id_time_prev=-1 !< time index for previous file, when using multiple files + integer :: id_time_next=-1 !< time index for next file, when using multiple files integer :: axis_sizes(4) character(len=32) :: axis_names(4) real(FMS_DATA_OVERRIDE_KIND_), dimension(:,:), pointer :: lon_local =>NULL() !< of output (target) grid cells @@ -763,6 +923,8 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d 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 !< first record of "current" file + type(time_type) :: last_record !< last record of "current" file integer :: curr_position !< position of the field currently processed in override_array real(FMS_DATA_OVERRIDE_KIND_) :: factor integer, dimension(4) :: comp_domain = 0 !< istart,iend,jstart,jend for compute domain @@ -775,6 +937,7 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d real(FMS_DATA_OVERRIDE_KIND_) :: lat_min, lat_max integer :: is_src, ie_src, js_src, je_src logical :: exists + logical :: multifile !< use multiple consecutive files for override 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 @@ -808,6 +971,7 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d fieldname = data_table(index1)%fieldname_file ! fieldname in netCDF data file factor = data_table(index1)%factor + multifile = data_table(index1)%multifile if(fieldname == "") then return_data = factor @@ -816,6 +980,8 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d 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') @@ -893,21 +1059,103 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d filename = filename2 endif + ! if using consecutive files for data_override, get file names + if_multi3: if (multifile) then + if_prev3: if (trim(prevfilename) /= '') then + inquire(file=trim(prevfilename),EXIST=exists) + if (.not. exists) then + call get_mosaic_tile_file(prevfilename,prevfilename2,.false.,domain) + prevfilename = prevfilename2 + endif + endif if_prev3 + if_next3: if (trim(nextfilename) /= '') then + inquire(file=trim(nextfilename),EXIST=exists) + if (.not. exists) then + call get_mosaic_tile_file(nextfilename,nextfilename2,.false.,domain) + nextfilename = nextfilename2 + endif + endif if_next3 + endif if_multi3 + !--- 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 using consecutive files for data_override, get time axis for previous and next files + ! and check spatial dims for consistency + if_multi4: if (multifile) then + id_time_prev = -1 + if_prev4:if (trim(prevfilename) /= '') then + id_time_prev = init_external_field(prevfilename,fieldname,domain=domain, & + verbose=.false.,use_comp_domain=use_comp_domain, & + nwindows = nwindows, ongrid=ongrid) + dims = get_external_field_size(id_time) + prev_dims = get_external_field_size(id_time_prev) + ! check consistency of spatial dims + if ((prev_dims(1) .ne. dims(1)) .or. (prev_dims(2) .ne. dims(2)) .or. & + (prev_dims(3) .ne. dims(3))) then + call mpp_error(FATAL, 'data_override: dimensions mismatch between consecutive forcing files') + endif + allocate(data_table(index1)%time_prev_records(prev_dims(4))) + call get_time_axis(id_time_prev,data_table(index1)%time_prev_records) + endif if_prev4 + id_time_next = -1 + if_next4: if (trim(nextfilename) /= '') then + id_time_next = init_external_field(nextfilename,fieldname,domain=domain, & + verbose=.false.,use_comp_domain=use_comp_domain, & + nwindows = nwindows, ongrid=ongrid) + dims = get_external_field_size(id_time) + next_dims = get_external_field_size(id_time_next) + ! check consistency of spatial dims + if ((next_dims(1) .ne. dims(1)) .or. (next_dims(2) .ne. dims(2)) .or. & + (next_dims(3) .ne. dims(3))) then + call mpp_error(FATAL, 'data_override: dimensions mismatch between consecutive forcing files') + endif + allocate(data_table(index1)%time_next_records(next_dims(4))) + call get_time_axis(id_time_next,data_table(index1)%time_next_records) + endif if_next4 + endif if_multi4 + 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 using consecutive files for data_override, get time axis for previous and next files + ! and check spatial dims for consistency + if_multi5: if (multifile) then + id_time_prev = -1 + if_prev5: if (trim(prevfilename) /= '') 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 if_prev5 + id_time_next = -1 + if_next5: if (trim(nextfilename) /= '') 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 if_next5 + endif if_multi5 + 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) @@ -975,6 +1223,14 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d 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 + if (trim(prevfilename) /= '') then + call reset_src_data_region(id_time_prev, is_src, ie_src, js_src, je_src) + endif + if (trim(nextfilename) /= '') then + call reset_src_data_region(id_time_next, is_src, ie_src, js_src, je_src) + endif + 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 @@ -1001,10 +1257,18 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d 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 + if (trim(prevfilename) /= '') then + call set_override_region(id_time_prev, data_table(index1)%region_type, istart, iend, jstart, jend) + endif + if (trim(nextfilename) /= '') then + call set_override_region(id_time_next, data_table(index1)%region_type, istart, iend, jstart, jend) + endif + endif deallocate(lon_tmp, lat_tmp) endif - endif + endif !ongrid else !curr_position >0 dims = override_array(curr_position)%dims comp_domain = override_array(curr_position)%comp_domain @@ -1021,6 +1285,8 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d 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 @@ -1093,6 +1359,14 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d if(dims(3) .NE. 1 .and. (size(return_data,3) .NE. dims(3))) & call mpp_error(FATAL, "data_override: dims(3) .NE. 1 and size(return_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 @@ -1107,13 +1381,85 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d !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,return_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 using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi6: if (multifile) then + if_time6: 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,return_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, do not use bridge + call time_interp_external(id_time,time,return_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 if_time6 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,return_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 if_multi6 + 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,return_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 using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi7: if (multifile) then + if_time7: 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,& + return_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, do not use bridge + call time_interp_external(id_time,time,return_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 if_time7 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,return_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 if_multi7 + end if return_data(:,:,1) = return_data(:,:,1)*factor do i = 2, size(return_data,3) @@ -1121,13 +1467,73 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d end do else if (use_comp_domain) then - call time_interp_external(id_time,time,return_data,verbose=.false., & - is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + + ! if using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi8: if (multifile) then + if_time8: 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,return_data,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, do not use bridge + call time_interp_external(id_time,time,return_data,verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif if_time8 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,return_data,verbose=.false., & + is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id) + endif if_multi8 + 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,return_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 using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi9: if (multifile) then + if_time9: 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,& + return_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 ! first_record <= time <= last_record, do not use bridge + call time_interp_external(id_time,time,return_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 if_time9 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,return_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 if_multi9 + end if return_data = return_data*factor endif @@ -1135,9 +1541,41 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d ! 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,return_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 using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi10: if (multifile) then + if_time10: 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,return_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 ! first_record <= time <= last_record, do not use bridge + call time_interp_external(id_time,time,return_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 if_time10 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,return_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 if_multi10 + return_data(:,:,1) = return_data(:,:,1)*factor do i = 2, size(return_data,3) return_data(:,:,i) = return_data(:,:,1) @@ -1145,10 +1583,45 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d else allocate(mask_out(size(return_data,1), size(return_data,2),1)) mask_out = .false. - call time_interp_external(id_time,time,return_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 using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi11: if (multifile) then + if_time11: 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,return_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 ! first_record <= time <= last_record, do not use bridge + call time_interp_external(id_time,time,return_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 if_time11 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,return_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 if_multi11 + where(mask_out(:,:,1)) return_data(:,:,1) = return_data(:,:,1)*factor end where @@ -1161,17 +1634,84 @@ subroutine DATA_OVERRIDE_3D_(gridname,fieldname_code,return_data,time,override,d endif else if( data_table(index1)%region_type == NO_REGION ) then - call time_interp_external(id_time,time,return_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) - return_data = return_data*factor + + ! if using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi12: if (multifile) then + if_time12: 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,return_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 ! first_record <= time <= last_record, do not use bridge + call time_interp_external(id_time,time,return_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 if_time12 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,return_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 if_multi12 + + return_data = return_data*factor else allocate(mask_out(size(return_data,1), size(return_data,2), size(return_data,3)) ) mask_out = .false. - call time_interp_external(id_time,time,return_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 using consecutive files, allow to perform time interpolation between the last record of previous + ! file and first record of current file OR between the last record of current file and first record of + ! next file hence "bridging" over files. + if_multi13: if (multifile) then + if_time13: 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,return_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 ! first_record <= time <= last_record, do not use bridge + call time_interp_external(id_time,time,return_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 if_time13 + else ! standard behavior + if ((timelast_record)) then + call mpp_error(WARNING, & + 'data_override: current time outside bounds, use [previous]:current:[next] files in data_table') + endif + call time_interp_external(id_time,time,return_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 if_multi13 + where(mask_out) return_data = return_data*factor end where diff --git a/time_interp/Makefile.am b/time_interp/Makefile.am index a4955bc20a..9097254cbc 100644 --- a/time_interp/Makefile.am +++ b/time_interp/Makefile.am @@ -37,6 +37,9 @@ libtime_interp_la_SOURCES = \ include/time_interp_r4.fh \ include/time_interp_r8.fh \ include/time_interp.inc \ + include/time_interp_external2_bridge_r4.fh \ + include/time_interp_external2_bridge_r8.fh \ + include/time_interp_external2_bridge.inc \ include/time_interp_external2_r4.fh \ include/time_interp_external2_r8.fh \ include/time_interp_external2.inc diff --git a/time_interp/include/time_interp_external2_bridge.inc b/time_interp/include/time_interp_external2_bridge.inc new file mode 100644 index 0000000000..e550b42296 --- /dev/null +++ b/time_interp/include/time_interp_external2_bridge.inc @@ -0,0 +1,305 @@ +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the GFDL Flexible Modeling System (FMS). +!* +!* FMS is free software: you can redistribute it and/or modify it under +!* the terms of the GNU Lesser General Public License as published by +!* the Free Software Foundation, either version 3 of the License, or (at +!* your option) any later version. +!* +!* FMS is distributed in the hope that it will be useful, but WITHOUT +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +!* for more details. +!* +!* You should have received a copy of the GNU Lesser General Public +!* License along with FMS. If not, see . +!*********************************************************************** +!> @ingroup time_interp +!> @addtogroup time_interp_external2_mod +!> @{ + + !> @brief 2D interpolation for @ref time_interp_external_bridge + 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 !< index of first external field + integer, intent(in) :: index2 !< index of second external field + type(time_type), intent(in) :: time !< target time for data + real(FMS_TI_KIND_), dimension(:,:), intent(inout) :: data_in !< global or local data array + integer, intent(in), optional :: interp !< hardcoded to linear + logical, intent(in), optional :: verbose !< flag for debugging + type(horiz_interp_type),intent(in), optional :: horz_interp !< horizontal interpolation type + 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 !< horizontal indices for load_record + integer, intent(in), optional :: window_id !< harcoded to 1 in load_record + + real(FMS_TI_KIND_), dimension(size(data_in,1), size(data_in,2), 1) :: data_out !< 3d version of data_in + logical, dimension(size(data_in,1), size(data_in,2), 1) :: mask3d !< 3d version of mask_out + + 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(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_ + + !> 3D interpolation for @ref time_interp_external + !! Provide data from external file interpolated to current model time. + !! Data may be local to current processor or global, depending on + !! "init_external_field" flags. + subroutine TIME_INTERP_EXTERNAL_BRIDGE_3D_(index1, index2, time, time_data, interp,verbose,horz_interp, mask_out, & + & is_in, ie_in, js_in, je_in, window_id) + + integer, intent(in) :: index1 !< index of first external field + integer, intent(in) :: index2 !< index of second external field + type(time_type), intent(in) :: time !< target time for data + real(FMS_TI_KIND_), dimension(:,:,:), intent(inout) :: time_data !< global or local data array + integer, intent(in), optional :: interp !< hardcoded to linear + logical, intent(in), optional :: verbose !< flag for debugging + type(horiz_interp_type), intent(in), optional :: horz_interp !< horizontal interpolation type + logical, dimension(:,:,:), intent(out), optional :: mask_out !< set to true where output data is valid + integer, intent(in), optional :: is_in, ie_in !< x horizontal indices for load_record + integer, intent(in), optional :: js_in, je_in !< y horizontal indices for load_record + integer, intent(in), optional :: window_id !< harcoded to 1 in load_record + + type(time_type) :: time1 !< time type associated with index1 of external field + type(time_type) :: time2 !< time type associated with index2 of external field + integer :: dims1(4) !< dimensions XYZT of index1 + integer :: dims2(4) !< dimensions XYZT of index2 + integer :: nx, ny, nz !< size in X,Y,Z of array + integer :: interp_method !< hardcoded to linear in time + integer :: t1, t2 !< temporary to store time index + integer :: i1, i2 !< temporary to store time index + integer :: isc, iec, jsc, jec !< start/end arrays in X,Y + integer :: isc1, iec1, jsc1, jec1 !< start/end arrays in X,Y for field1 + integer :: isc2, iec2, jsc2, jec2 !< start/end arrays in X,Y for field2 + integer :: yy, mm, dd, hh, minu, ss !< year, month, day, hour, minute, second for date + + 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(FMS_TI_KIND_) :: w1 !< interp weight for index1 + real(FMS_TI_KIND_) :: w2 !< interp weight for index2 + logical :: verb !< verbose + character(len=16) :: message1 !< temp string + character(len=16) :: message2 !< temp string + integer, parameter :: kindl = FMS_TI_KIND_ + + nx = size(time_data,1) + ny = size(time_data,2) + nz = size(time_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 .or. index1 > num_fields) & + call mpp_error(FATAL,'invalid index1 in call to time_interp_external_bridge:' // & + '-- field was not initialized or failed to initialize') + if (index2 < 1 .or. index2 > num_fields) & + call mpp_error(FATAL,'invalid index2 in call to time_interp_external_bridge:' // & + ' -- field was not initialized or failed to initialize') + + isc1=loaded_fields(index1)%isc;iec1=loaded_fields(index1)%iec + jsc1=loaded_fields(index1)%jsc;jec1=loaded_fields(index1)%jec + isc2=loaded_fields(index2)%isc;iec2=loaded_fields(index2)%iec + jsc2=loaded_fields(index2)%jsc;jec2=loaded_fields(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(loaded_fields(index2)%name) /= trim(loaded_fields(index1)%name)) then + call mpp_error(FATAL, 'time_interp_external_bridge: cannot bridge between different fields.' // & + 'field1='//trim(loaded_fields(index1)%name)//' field2='//trim(loaded_fields(index2)%name)) + endif + + if ((loaded_fields(index1)%numwindows == 1) .and. (loaded_fields(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: is_in, ie_in, js_in and je_in must be present ' // & + 'when numwindows > 1, field='//trim(loaded_fields(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 < loaded_fields(index1)%siz(3)) then + write(message1,'(i6,2i5)') nx,ny,nz + call mpp_error(FATAL,'field '//trim(loaded_fields(index1)%name)// & + ' Array size mismatch in time_interp_external_bridge.'// & + ' Array "data" is too small. shape(time_data)='//message1) + endif + if (nx < nxw .or. ny < nyw .or. nz < loaded_fields(index2)%siz(3)) then + write(message1,'(i6,2i5)') nx,ny,nz + call mpp_error(FATAL,'field '//trim(loaded_fields(index2)%name)// & + ' Array size mismatch in time_interp_external_bridge.'// & + ' Array "data" is too small. shape(time_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(loaded_fields(index1)%name)// & + ' array size mismatch in time_interp_external_bridge.'// & + ' Shape of array "mask_out" does not match that of array "data".'// & + ' shape(time_data)='//message1//' shape(mask_out)='//message2) + endif + endif + + if ((loaded_fields(index1)%have_modulo_times) .or. (loaded_fields(index2)%have_modulo_times)) then + call mpp_error(FATAL, 'time_interp_external_bridge: field '//trim(loaded_fields(index1)%name)// & + ' array cannot have modulo time') + endif + if ((loaded_fields(index1)%modulo_time) .or. (loaded_fields(index2)%modulo_time)) then + call mpp_error(FATAL, 'time_interp_external_bridge: field '//trim(loaded_fields(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 = loaded_fields(index1)%time(t1) + time2 = loaded_fields(index2)%time(t2) + w2= (time - time1) // (time2 - time1) + w1 = 1.0_kindl-w2 + + if (verb) then + call get_date(time,yy,mm,dd,hh,minu,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,':',minu,':',ss + write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2 + endif + + call load_record(loaded_fields(index1),t1,horz_interp, is_in, ie_in ,js_in, je_in, window_id) + call load_record(loaded_fields(index2),t2,horz_interp, is_in, ie_in ,js_in, je_in, window_id) + i1 = find_buf_index(t1,loaded_fields(index1)%ibuf) + i2 = find_buf_index(t2,loaded_fields(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= ',loaded_fields(index2)%ibuf + write(outunit,*) 'i1,i2= ',i1, i2 + endif + + if((loaded_fields(index1)%region_type == NO_REGION) .and. (loaded_fields(index2)%region_type == NO_REGION) ) then + where(loaded_fields(index1)%mask(isc:iec,jsc:jec,:,i1).and.loaded_fields(index2)%mask(isc:iec,jsc:jec,:,i2)) + time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index1)%domain_data(isc:iec,jsc:jec,:,i1), kindl)*w1 + & + real(loaded_fields(index2)%domain_data(isc:iec,jsc:jec,:,i2), kindl)*w2 + elsewhere + time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index1)%missing, kindl) + end where + else + where(loaded_fields(index1)%mask(isc:iec,jsc:jec,:,i1).and.loaded_fields(index2)%mask(isc:iec,jsc:jec,:,i2)) + time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index1)%domain_data(isc:iec,jsc:jec,:,i1), kindl)*w1 + & + real(loaded_fields(index2)%domain_data(isc:iec,jsc:jec,:,i2), kindl)*w2 + end where + endif + + if(PRESENT(mask_out)) & + mask_out(isw:iew,jsw:jew,:) = & + loaded_fields(index1)%mask(isc:iec,jsc:jec,:,i1).and.& + loaded_fields(index2)%mask(isc:iec,jsc:jec,:,i2) + + end subroutine TIME_INTERP_EXTERNAL_BRIDGE_3D_ + + subroutine TIME_INTERP_EXTERNAL_BRIDGE_0D_(index1, index2, time, time_data, verbose) + + integer, intent(in) :: index1 !< index of first external field + integer, intent(in) :: index2 !< index of second external field + type(time_type), intent(in) :: time !< target time for data + real(FMS_TI_KIND_), intent(inout) :: time_data !< global or local data array + logical, intent(in), optional :: verbose !< flag for debugging + + integer :: t1, t2 !< temporary to store time index + integer :: i1, i2 !< temporary to store time index + integer :: yy, mm, dd, hh, minu, ss !< year, month, day, hour, minute, second for date + type(time_type) :: time1 !< time type associated with index1 of external field + type(time_type) :: time2 !< time type associated with index2 of external field + integer :: dims1(4) !< dimensions XYZT of index1 + integer :: dims2(4) !< dimensions XYZT of index2 + + real(FMS_TI_KIND_) :: w1 !< interp weight for index1 + real(FMS_TI_KIND_) :: w2 !< interp weight for index2 + logical :: verb !< verbose + integer, parameter :: kindl = FMS_TI_KIND_ + + verb=.false. + if (PRESENT(verbose)) verb=verbose + if (debug_this_module) verb = .true. + + if (index1 < 0 .or. index1 > num_fields) & + call mpp_error(FATAL,'invalid index1 in call to time_interp_ext' // & + ' -- field was not initialized or failed to initialize') + if (index2 < 0 .or. index2 > num_fields) & + call mpp_error(FATAL,'invalid index2 in call to time_interp_ext' // & + ' -- field was not initialized or failed to initialize') + + if ((loaded_fields(index1)%have_modulo_times) .or. (loaded_fields(index2)%have_modulo_times)) then + call mpp_error(FATAL, 'time_interp_external_bridge: field '//trim(loaded_fields(index1)%name)// & + ' array cannot have modulo time') + endif + if ((loaded_fields(index1)%modulo_time) .or. (loaded_fields(index2)%modulo_time)) then + call mpp_error(FATAL, 'time_interp_external_bridge: field '//trim(loaded_fields(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 = loaded_fields(index1)%time(t1) + time2 = loaded_fields(index2)%time(t2) + w2= (time - time1) // (time2 - time1) + w1 = 1.0_kindl-w2 + + if (verb) then + call get_date(time,yy,mm,dd,hh,minu,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,':',minu,':',ss + write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2 + endif + call load_record_0d(loaded_fields(index2),t1) + call load_record_0d(loaded_fields(index2),t2) + i1 = find_buf_index(t1,loaded_fields(index2)%ibuf) + i2 = find_buf_index(t2,loaded_fields(index2)%ibuf) + + if(i1<0.or.i2<0) & + call mpp_error(FATAL,'time_interp_external : records were not loaded correctly in memory') + time_data = real(loaded_fields(index2)%domain_data(1,1,1,i1), FMS_TI_KIND_)*w1 & + + real(loaded_fields(index2)%domain_data(1,1,1,i2), FMS_TI_KIND_)*w2 + if (verb) then + write(outunit,*) 'ibuf= ',loaded_fields(index2)%ibuf + write(outunit,*) 'i1,i2= ',i1, i2 + endif + + end subroutine TIME_INTERP_EXTERNAL_BRIDGE_0D_ + + +! ============================================================================ diff --git a/time_interp/include/time_interp_external2_bridge_r4.fh b/time_interp/include/time_interp_external2_bridge_r4.fh new file mode 100644 index 0000000000..1e8fafe4d1 --- /dev/null +++ b/time_interp/include/time_interp_external2_bridge_r4.fh @@ -0,0 +1,31 @@ +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the GFDL Flexible Modeling System (FMS). +!* +!* FMS is free software: you can redistribute it and/or modify it under +!* the terms of the GNU Lesser General Public License as published by +!* the Free Software Foundation, either version 3 of the License, or (at +!* your option) any later version. +!* +!* FMS is distributed in the hope that it will be useful, but WITHOUT +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +!* for more details. +!* +!* You should have received a copy of the GNU Lesser General Public +!* License along with FMS. If not, see . +!*********************************************************************** +#undef FMS_TI_KIND_ +#define FMS_TI_KIND_ r4_kind + +#undef TIME_INTERP_EXTERNAL_BRIDGE_0D_ +#define TIME_INTERP_EXTERNAL_BRIDGE_0D_ time_interp_external_bridge_0d_r4 + +#undef TIME_INTERP_EXTERNAL_BRIDGE_2D_ +#define TIME_INTERP_EXTERNAL_BRIDGE_2D_ time_interp_external_bridge_2d_r4 + +#undef TIME_INTERP_EXTERNAL_BRIDGE_3D_ +#define TIME_INTERP_EXTERNAL_BRIDGE_3D_ time_interp_external_bridge_3d_r4 + +#include "time_interp_external2_bridge.inc" diff --git a/time_interp/include/time_interp_external2_bridge_r8.fh b/time_interp/include/time_interp_external2_bridge_r8.fh new file mode 100644 index 0000000000..ef6d2fe209 --- /dev/null +++ b/time_interp/include/time_interp_external2_bridge_r8.fh @@ -0,0 +1,31 @@ +!*********************************************************************** +!* GNU Lesser General Public License +!* +!* This file is part of the GFDL Flexible Modeling System (FMS). +!* +!* FMS is free software: you can redistribute it and/or modify it under +!* the terms of the GNU Lesser General Public License as published by +!* the Free Software Foundation, either version 3 of the License, or (at +!* your option) any later version. +!* +!* FMS is distributed in the hope that it will be useful, but WITHOUT +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +!* for more details. +!* +!* You should have received a copy of the GNU Lesser General Public +!* License along with FMS. If not, see . +!*********************************************************************** +#undef FMS_TI_KIND_ +#define FMS_TI_KIND_ r8_kind + +#undef TIME_INTERP_EXTERNAL_BRIDGE_0D_ +#define TIME_INTERP_EXTERNAL_BRIDGE_0D_ time_interp_external_bridge_0d_r8 + +#undef TIME_INTERP_EXTERNAL_BRIDGE_2D_ +#define TIME_INTERP_EXTERNAL_BRIDGE_2D_ time_interp_external_bridge_2d_r8 + +#undef TIME_INTERP_EXTERNAL_BRIDGE_3D_ +#define TIME_INTERP_EXTERNAL_BRIDGE_3D_ time_interp_external_bridge_3d_r8 + +#include "time_interp_external2_bridge.inc" diff --git a/time_interp/time_interp_external2.F90 b/time_interp/time_interp_external2.F90 index 67e5127188..02fad81f4b 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 @@ -150,6 +151,15 @@ module time_interp_external2_mod module procedure time_interp_external_3d_r8 end interface + interface time_interp_external_bridge + module procedure time_interp_external_bridge_0d_r4 + module procedure time_interp_external_bridge_2d_r4 + module procedure time_interp_external_bridge_3d_r4 + module procedure time_interp_external_bridge_0d_r8 + module procedure time_interp_external_bridge_2d_r8 + module procedure time_interp_external_bridge_3d_r8 + end interface + !> @addtogroup time_interp_external2_mod !> @{ @@ -1129,6 +1139,9 @@ end subroutine time_interp_external_exit #include "time_interp_external2_r4.fh" #include "time_interp_external2_r8.fh" +#include "time_interp_external2_bridge_r4.fh" +#include "time_interp_external2_bridge_r8.fh" + end module time_interp_external2_mod !> @} ! close documentation grouping