FMS  2025.04
Flexible Modeling System
time_interp_external2_bridge.inc
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 !> @ingroup time_interp
19 !> @addtogroup time_interp_external2_mod
20 !> @{
21 
22  !> @brief 1D interpolation for @ref time_interp_external_bridge
23  subroutine time_interp_external_bridge_1d_(index1, index2, time, data_in, interp, verbose,horz_interp, mask_out, &
24  is_in, ie_in, window_id)
25 
26  integer, intent(in) :: index1 !< index of first external field
27  integer, intent(in) :: index2 !< index of second external field
28  type(time_type), intent(in) :: time !< target time for data
29  real(FMS_TI_KIND_), dimension(:), intent(inout) :: data_in !< global or local data array
30  integer, intent(in), optional :: interp !< hardcoded to linear
31  logical, intent(in), optional :: verbose !< flag for debugging
32  type(horiz_interp_type),intent(in), optional :: horz_interp !< horizontal interpolation type
33  logical, dimension(:), intent(out), optional :: mask_out !< set to true where output data is valid
34  integer, intent(in), optional :: is_in, ie_in !< horizontal indices for load_record
35  integer, intent(in), optional :: window_id !< harcoded to 1 in load_record
36 
37  real(FMS_TI_KIND_), dimension(size(data_in,1), 1, 1) :: data_out !< 3d version of data_in
38  logical, dimension(size(data_in,1), 1, 1) :: mask3d !< 3d version of mask_out
39 
40  data_out(:,1,1) = data_in(:) ! fill initial values for the portions of array that are not touched by 3d routine
41  call time_interp_external_bridge(index1, index2, time, data_out, interp, verbose, horz_interp, mask3d, &
42  is_in=is_in,ie_in=ie_in,js_in=1,je_in=1,window_id=window_id)
43  data_in(:) = data_out(:,1,1)
44  if (PRESENT(mask_out)) mask_out(:) = mask3d(:,1,1)
45 
46  end subroutine time_interp_external_bridge_1d_
47 
48 
49  !> @brief 2D interpolation for @ref time_interp_external_bridge
50  subroutine time_interp_external_bridge_2d_(index1, index2, time, data_in, interp, verbose,horz_interp, mask_out, &
51  is_in, ie_in, js_in, je_in, window_id)
52 
53  integer, intent(in) :: index1 !< index of first external field
54  integer, intent(in) :: index2 !< index of second external field
55  type(time_type), intent(in) :: time !< target time for data
56  real(FMS_TI_KIND_), dimension(:,:), intent(inout) :: data_in !< global or local data array
57  integer, intent(in), optional :: interp !< hardcoded to linear
58  logical, intent(in), optional :: verbose !< flag for debugging
59  type(horiz_interp_type),intent(in), optional :: horz_interp !< horizontal interpolation type
60  logical, dimension(:,:), intent(out), optional :: mask_out !< set to true where output data is valid
61  integer, intent(in), optional :: is_in, ie_in, js_in, je_in !< horizontal indices for load_record
62  integer, intent(in), optional :: window_id !< harcoded to 1 in load_record
63 
64  real(FMS_TI_KIND_), dimension(size(data_in,1), size(data_in,2), 1) :: data_out !< 3d version of data_in
65  logical, dimension(size(data_in,1), size(data_in,2), 1) :: mask3d !< 3d version of mask_out
66 
67  data_out(:,:,1) = data_in(:,:) ! fill initial values for the portions of array that are not touched by 3d routine
68  call time_interp_external_bridge(index1, index2, time, data_out, interp, verbose, horz_interp, mask3d, &
69  is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id)
70  data_in(:,:) = data_out(:,:,1)
71  if (PRESENT(mask_out)) mask_out(:,:) = mask3d(:,:,1)
72 
73  return
74  end subroutine time_interp_external_bridge_2d_
75 
76  !> 3D interpolation for @ref time_interp_external
77  !! Provide data from external file interpolated to current model time.
78  !! Data may be local to current processor or global, depending on
79  !! "init_external_field" flags.
80  subroutine time_interp_external_bridge_3d_(index1, index2, time, time_data, interp,verbose,horz_interp, mask_out, &
81  & is_in, ie_in, js_in, je_in, window_id)
82 
83  integer, intent(in) :: index1 !< index of first external field
84  integer, intent(in) :: index2 !< index of second external field
85  type(time_type), intent(in) :: time !< target time for data
86  real(FMS_TI_KIND_), dimension(:,:,:), intent(inout) :: time_data !< global or local data array
87  integer, intent(in), optional :: interp !< hardcoded to linear
88  logical, intent(in), optional :: verbose !< flag for debugging
89  type(horiz_interp_type), intent(in), optional :: horz_interp !< horizontal interpolation type
90  logical, dimension(:,:,:), intent(out), optional :: mask_out !< set to true where output data is valid
91  integer, intent(in), optional :: is_in, ie_in !< x horizontal indices for load_record
92  integer, intent(in), optional :: js_in, je_in !< y horizontal indices for load_record
93  integer, intent(in), optional :: window_id !< harcoded to 1 in load_record
94 
95  type(time_type) :: time1 !< time type associated with index1 of external field
96  type(time_type) :: time2 !< time type associated with index2 of external field
97  integer :: dims1(4) !< dimensions XYZT of index1
98  integer :: dims2(4) !< dimensions XYZT of index2
99  integer :: nx, ny, nz !< size in X,Y,Z of array
100  integer :: interp_method !< hardcoded to linear in time
101  integer :: t1, t2 !< temporary to store time index
102  integer :: i1, i2 !< temporary to store time index
103  integer :: isc, iec, jsc, jec !< start/end arrays in X,Y
104  integer :: isc1, iec1, jsc1, jec1 !< start/end arrays in X,Y for field1
105  integer :: isc2, iec2, jsc2, jec2 !< start/end arrays in X,Y for field2
106  integer :: yy, mm, dd, hh, minu, ss !< year, month, day, hour, minute, second for date
107 
108  integer :: isw, iew, jsw, jew, nxw, nyw !< these are boundaries of the updated portion of the "data" argument
109  !! they are calculated using sizes of the "data" and isc,iec,jsc,jsc
110  !! fileds from respective input field, to center the updated portion
111  !! in the output array
112 
113  real(FMS_TI_KIND_) :: w1 !< interp weight for index1
114  real(FMS_TI_KIND_) :: w2 !< interp weight for index2
115  logical :: verb !< verbose
116  character(len=16) :: message1 !< temp string
117  character(len=16) :: message2 !< temp string
118  integer, parameter :: kindl = fms_ti_kind_
119 
120  nx = size(time_data,1)
121  ny = size(time_data,2)
122  nz = size(time_data,3)
123 
124  interp_method = linear_time_interp
125  if (PRESENT(interp)) interp_method = interp
126  verb=.false.
127  if (PRESENT(verbose)) verb=verbose
128  if (debug_this_module) verb = .true.
129 
130  if (index1 < 1 .or. index1 > num_fields) &
131  call mpp_error(fatal,'invalid index1 in call to time_interp_external_bridge:' // &
132  '-- field was not initialized or failed to initialize')
133  if (index2 < 1 .or. index2 > num_fields) &
134  call mpp_error(fatal,'invalid index2 in call to time_interp_external_bridge:' // &
135  ' -- field was not initialized or failed to initialize')
136 
137  isc1=loaded_fields(index1)%isc;iec1=loaded_fields(index1)%iec
138  jsc1=loaded_fields(index1)%jsc;jec1=loaded_fields(index1)%jec
139  isc2=loaded_fields(index2)%isc;iec2=loaded_fields(index2)%iec
140  jsc2=loaded_fields(index2)%jsc;jec2=loaded_fields(index2)%jec
141 
142  if ((isc1 /= isc2) .or. (iec1 /= iec2) .or. (jsc1 /= jsc2) .or. (jec1 /= jec2)) then
143  call mpp_error(fatal, 'time_interp_external_bridge: dimensions must be the same in both index1 and index2 ')
144  endif
145 
146  isc=isc1 ; iec=iec1 ; jsc=jsc1 ; jec=jec1
147 
148  if (trim(loaded_fields(index2)%name) /= trim(loaded_fields(index1)%name)) then
149  call mpp_error(fatal, 'time_interp_external_bridge: cannot bridge between different fields.' // &
150  'field1='//trim(loaded_fields(index1)%name)//' field2='//trim(loaded_fields(index2)%name))
151  endif
152 
153  if ((loaded_fields(index1)%numwindows == 1) .and. (loaded_fields(index2)%numwindows == 1)) then
154  nxw = iec-isc+1
155  nyw = jec-jsc+1
156  else
157  if(.not. present(is_in) .or. .not. present(ie_in) .or. .not. present(js_in) .or. .not. present(je_in))then
158  call mpp_error(fatal, 'time_interp_external: is_in, ie_in, js_in and je_in must be present ' // &
159  'when numwindows > 1, field='//trim(loaded_fields(index1)%name))
160  endif
161  nxw = ie_in - is_in + 1
162  nyw = je_in - js_in + 1
163  isc = isc + is_in - 1
164  iec = isc + ie_in - is_in
165  jsc = jsc + js_in - 1
166  jec = jsc + je_in - js_in
167  endif
168 
169  isw = (nx-nxw)/2+1; iew = isw+nxw-1
170  jsw = (ny-nyw)/2+1; jew = jsw+nyw-1
171 
172  if (nx < nxw .or. ny < nyw .or. nz < loaded_fields(index1)%siz(3)) then
173  write(message1,'(i6,2i5)') nx,ny,nz
174  call mpp_error(fatal,'field '//trim(loaded_fields(index1)%name)// &
175  ' Array size mismatch in time_interp_external_bridge.'// &
176  ' Array "data" is too small. shape(time_data)='//message1)
177  endif
178  if (nx < nxw .or. ny < nyw .or. nz < loaded_fields(index2)%siz(3)) then
179  write(message1,'(i6,2i5)') nx,ny,nz
180  call mpp_error(fatal,'field '//trim(loaded_fields(index2)%name)// &
181  ' Array size mismatch in time_interp_external_bridge.'// &
182  ' Array "data" is too small. shape(time_data)='//message1)
183  endif
184 
185  if(PRESENT(mask_out)) then
186  if (size(mask_out,1) /= nx .or. size(mask_out,2) /= ny .or. size(mask_out,3) /= nz) then
187  write(message1,'(i6,2i5)') nx,ny,nz
188  write(message2,'(i6,2i5)') size(mask_out,1),size(mask_out,2),size(mask_out,3)
189  call mpp_error(fatal,'field '//trim(loaded_fields(index1)%name)// &
190  ' array size mismatch in time_interp_external_bridge.'// &
191  ' Shape of array "mask_out" does not match that of array "data".'// &
192  ' shape(time_data)='//message1//' shape(mask_out)='//message2)
193  endif
194  endif
195 
196  if ((loaded_fields(index1)%have_modulo_times) .or. (loaded_fields(index2)%have_modulo_times)) then
197  call mpp_error(fatal, 'time_interp_external_bridge: field '//trim(loaded_fields(index1)%name)// &
198  ' array cannot have modulo time')
199  endif
200  if ((loaded_fields(index1)%modulo_time) .or. (loaded_fields(index2)%modulo_time)) then
201  call mpp_error(fatal, 'time_interp_external_bridge: field '//trim(loaded_fields(index1)%name)// &
202  ' array cannot have modulo time')
203  endif
204 
205  dims1 = get_external_field_size(index1)
206  dims2 = get_external_field_size(index2)
207 
208  t1 = dims1(4)
209  t2 = 1
210 
211  time1 = loaded_fields(index1)%time(t1)
212  time2 = loaded_fields(index2)%time(t2)
213  w2= (time - time1) // (time2 - time1)
214  w1 = 1.0_kindl-w2
215 
216  if (verb) then
217  call get_date(time,yy,mm,dd,hh,minu,ss)
218  write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
219  'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',minu,':',ss
220  write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
221  endif
222 
223  call load_record(loaded_fields(index1),t1,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
224  call load_record(loaded_fields(index2),t2,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
225  i1 = find_buf_index(t1,loaded_fields(index1)%ibuf)
226  i2 = find_buf_index(t2,loaded_fields(index2)%ibuf)
227  if(i1<0.or.i2<0) &
228  call mpp_error(fatal,'time_interp_external_bridge : records were not loaded correctly in memory')
229 
230  if (verb) then
231  write(outunit,*) 'ibuf= ',loaded_fields(index2)%ibuf
232  write(outunit,*) 'i1,i2= ',i1, i2
233  endif
234 
235  if((loaded_fields(index1)%region_type == no_region) .and. (loaded_fields(index2)%region_type == no_region) ) then
236  where(loaded_fields(index1)%mask(isc:iec,jsc:jec,:,i1).and.loaded_fields(index2)%mask(isc:iec,jsc:jec,:,i2))
237  time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index1)%domain_data(isc:iec,jsc:jec,:,i1), kindl)*w1 + &
238  real(loaded_fields(index2)%domain_data(isc:iec,jsc:jec,:,i2), kindl)*w2
239  elsewhere
240  time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index1)%missing, kindl)
241  end where
242  else
243  where(loaded_fields(index1)%mask(isc:iec,jsc:jec,:,i1).and.loaded_fields(index2)%mask(isc:iec,jsc:jec,:,i2))
244  time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index1)%domain_data(isc:iec,jsc:jec,:,i1), kindl)*w1 + &
245  real(loaded_fields(index2)%domain_data(isc:iec,jsc:jec,:,i2), kindl)*w2
246  end where
247  endif
248 
249  if(PRESENT(mask_out)) &
250  mask_out(isw:iew,jsw:jew,:) = &
251  loaded_fields(index1)%mask(isc:iec,jsc:jec,:,i1).and.&
252  loaded_fields(index2)%mask(isc:iec,jsc:jec,:,i2)
253 
254  end subroutine time_interp_external_bridge_3d_
255 
256  subroutine time_interp_external_bridge_0d_(index1, index2, time, time_data, verbose)
257 
258  integer, intent(in) :: index1 !< index of first external field
259  integer, intent(in) :: index2 !< index of second external field
260  type(time_type), intent(in) :: time !< target time for data
261  real(FMS_TI_KIND_), intent(inout) :: time_data !< global or local data array
262  logical, intent(in), optional :: verbose !< flag for debugging
263 
264  integer :: t1, t2 !< temporary to store time index
265  integer :: i1, i2 !< temporary to store time index
266  integer :: yy, mm, dd, hh, minu, ss !< year, month, day, hour, minute, second for date
267  type(time_type) :: time1 !< time type associated with index1 of external field
268  type(time_type) :: time2 !< time type associated with index2 of external field
269  integer :: dims1(4) !< dimensions XYZT of index1
270  integer :: dims2(4) !< dimensions XYZT of index2
271 
272  real(FMS_TI_KIND_) :: w1 !< interp weight for index1
273  real(FMS_TI_KIND_) :: w2 !< interp weight for index2
274  logical :: verb !< verbose
275  integer, parameter :: kindl = fms_ti_kind_
276 
277  verb=.false.
278  if (PRESENT(verbose)) verb=verbose
279  if (debug_this_module) verb = .true.
280 
281  if (index1 < 0 .or. index1 > num_fields) &
282  call mpp_error(fatal,'invalid index1 in call to time_interp_ext' // &
283  ' -- field was not initialized or failed to initialize')
284  if (index2 < 0 .or. index2 > num_fields) &
285  call mpp_error(fatal,'invalid index2 in call to time_interp_ext' // &
286  ' -- field was not initialized or failed to initialize')
287 
288  if ((loaded_fields(index1)%have_modulo_times) .or. (loaded_fields(index2)%have_modulo_times)) then
289  call mpp_error(fatal, 'time_interp_external_bridge: field '//trim(loaded_fields(index1)%name)// &
290  ' array cannot have modulo time')
291  endif
292  if ((loaded_fields(index1)%modulo_time) .or. (loaded_fields(index2)%modulo_time)) then
293  call mpp_error(fatal, 'time_interp_external_bridge: field '//trim(loaded_fields(index1)%name)// &
294  ' array cannot have modulo time')
295  endif
296 
297  dims1 = get_external_field_size(index1)
298  dims2 = get_external_field_size(index2)
299 
300  t1 = dims1(4)
301  t2 = 1
302 
303  time1 = loaded_fields(index1)%time(t1)
304  time2 = loaded_fields(index2)%time(t2)
305  w2= (time - time1) // (time2 - time1)
306  w1 = 1.0_kindl-w2
307 
308  if (verb) then
309  call get_date(time,yy,mm,dd,hh,minu,ss)
310  write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
311  'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',minu,':',ss
312  write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
313  endif
314  call load_record_0d(loaded_fields(index2),t1)
315  call load_record_0d(loaded_fields(index2),t2)
316  i1 = find_buf_index(t1,loaded_fields(index2)%ibuf)
317  i2 = find_buf_index(t2,loaded_fields(index2)%ibuf)
318 
319  if(i1<0.or.i2<0) &
320  call mpp_error(fatal,'time_interp_external : records were not loaded correctly in memory')
321  time_data = real(loaded_fields(index2)%domain_data(1,1,1,i1), fms_ti_kind_)*w1 &
322  + real(loaded_fields(index2)%domain_data(1,1,1,i2), fms_ti_kind_)*w2
323  if (verb) then
324  write(outunit,*) 'ibuf= ',loaded_fields(index2)%ibuf
325  write(outunit,*) 'i1,i2= ',i1, i2
326  endif
327 
328  end subroutine time_interp_external_bridge_0d_
329 
330 
331 ! ============================================================================
subroutine time_interp_external_bridge_1d_(index1, index2, time, data_in, interp, verbose, horz_interp, mask_out, is_in, ie_in, window_id)
1D interpolation for time_interp_external_bridge
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)
3D interpolation for time_interp_external Provide data from external file interpolated to current mod...
subroutine time_interp_external_bridge_0d_(index1, index2, time, time_data, verbose)
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)
2D interpolation for time_interp_external_bridge