FMS  2025.04
Flexible Modeling System
time_interp_external2.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
23  subroutine time_interp_external_1d_(index, time, data_in, interp, verbose,horz_interp, mask_out, &
24  is_in, ie_in, window_id)
25 
26  integer, intent(in) :: index
27  type(time_type), intent(in) :: time
28  real(FMS_TI_KIND_), dimension(:), intent(inout) :: data_in
29  integer, intent(in), optional :: interp
30  logical, intent(in), optional :: verbose
31  type(horiz_interp_type),intent(in), optional :: horz_interp
32  logical, dimension(:), intent(out), optional :: mask_out !< set to true where output data is valid
33  integer, intent(in), optional :: is_in, ie_in
34  integer, intent(in), optional :: window_id
35 
36  real(FMS_TI_KIND_), dimension(size(data_in,1), 1, 1) :: data_out
37  logical, dimension(size(data_in), 1, 1) :: mask3d
38 
39  data_out(:,1,1) = data_in(:) ! fill initial values for the portions of array that are not touched by 3d routine
40  call time_interp_external(index, time, data_out, interp, verbose, horz_interp, mask3d, &
41  is_in=is_in,ie_in=ie_in,js_in=1,je_in=1,window_id=window_id)
42  data_in(:) = data_out(:,1,1)
43  if (PRESENT(mask_out)) mask_out(:) = mask3d(:,1,1)
44 
45  end subroutine time_interp_external_1d_
46 
47  !> @brief 2D interpolation for @ref time_interp_external
48  subroutine time_interp_external_2d_(index, time, data_in, interp, verbose,horz_interp, mask_out, &
49  is_in, ie_in, js_in, je_in, window_id)
50 
51  integer, intent(in) :: index
52  type(time_type), intent(in) :: time
53  real(FMS_TI_KIND_), dimension(:,:), intent(inout) :: data_in
54  integer, intent(in), optional :: interp
55  logical, intent(in), optional :: verbose
56  type(horiz_interp_type),intent(in), optional :: horz_interp
57  logical, dimension(:,:), intent(out), optional :: mask_out !< set to true where output data is valid
58  integer, intent(in), optional :: is_in, ie_in, js_in, je_in
59  integer, intent(in), optional :: window_id
60 
61  real(FMS_TI_KIND_), dimension(size(data_in,1), size(data_in,2), 1) :: data_out
62  logical, dimension(size(data_in,1), size(data_in,2), 1) :: mask3d
63 
64  data_out(:,:,1) = data_in(:,:) ! fill initial values for the portions of array that are not touched by 3d routine
65  call time_interp_external(index, time, data_out, interp, verbose, horz_interp, mask3d, &
66  is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id)
67  data_in(:,:) = data_out(:,:,1)
68  if (PRESENT(mask_out)) mask_out(:,:) = mask3d(:,:,1)
69 
70  return
71  end subroutine time_interp_external_2d_
72 
73 
74 
75  !> 3D interpolation for @ref time_interp_external
76  !! Provide data from external file interpolated to current model time.
77  !! Data may be local to current processor or global, depending on
78  !! "init_external_field" flags.
79  subroutine time_interp_external_3d_(index, time, time_data, interp,verbose,horz_interp, mask_out, is_in, ie_in, &
80  & js_in, je_in, window_id)
81 
82  integer, intent(in) :: index !< index of external field from previous call
83  !! to init_external_field
84  type(time_type), intent(in) :: time !< target time for data
85  real(FMS_TI_KIND_), dimension(:,:,:), intent(inout) :: time_data !< global or local data array
86  integer, intent(in), optional :: interp
87  logical, intent(in), optional :: verbose !< flag for debugging
88  type(horiz_interp_type), intent(in), optional :: horz_interp
89  logical, dimension(:,:,:), intent(out), optional :: mask_out !< set to true where output data is valid
90  integer, intent(in), optional :: is_in, ie_in, js_in, je_in
91  integer, intent(in), optional :: window_id
92 
93  integer :: nx, ny, nz, interp_method, t1, t2
94  integer :: i1, i2, isc, iec, jsc, jec, mod_time
95  integer :: yy, mm, dd, hh, min, ss
96  character(len=256) :: err_msg
97  character(len=FMS_PATH_LEN) :: filename
98 
99  integer :: isw, iew, jsw, jew, nxw, nyw
100  ! these are boundaries of the updated portion of the "data" argument
101  ! they are calculated using sizes of the "data" and isc,iec,jsc,jsc
102  ! fileds from respective input field, to center the updated portion
103  ! in the output array
104 
105  real(FMS_TI_KIND_) :: w1,w2
106  logical :: verb
107  character(len=16) :: message1, message2
108  integer, parameter :: kindl = fms_ti_kind_
109 
110  nx = size(time_data,1)
111  ny = size(time_data,2)
112  nz = size(time_data,3)
113 
114  interp_method = linear_time_interp
115  if (PRESENT(interp)) interp_method = interp
116  verb=.false.
117  if (PRESENT(verbose)) verb=verbose
118  if (debug_this_module) verb = .true.
119 
120  if (index < 1.or.index > num_fields) &
121  call mpp_error(fatal, &
122  & 'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize')
123 
124  isc=loaded_fields(index)%isc;iec=loaded_fields(index)%iec
125  jsc=loaded_fields(index)%jsc;jec=loaded_fields(index)%jec
126 
127  if( loaded_fields(index)%numwindows == 1 ) then
128  nxw = iec-isc+1
129  nyw = jec-jsc+1
130  else
131  if(.not. present(is_in) .or. .not. present(ie_in) .or. .not. present(js_in) .or. .not. present(je_in))then
132  call mpp_error(fatal, 'time_interp_external: is_in, ie_in, js_in and je_in must be present ' // &
133  'when numwindows > 1, field='//trim(loaded_fields(index)%name))
134  endif
135  nxw = ie_in - is_in + 1
136  nyw = je_in - js_in + 1
137  isc = isc + is_in - 1
138  iec = isc + ie_in - is_in
139  jsc = jsc + js_in - 1
140  jec = jsc + je_in - js_in
141  endif
142 
143  isw = (nx-nxw)/2+1; iew = isw+nxw-1
144  jsw = (ny-nyw)/2+1; jew = jsw+nyw-1
145 
146  if (nx < nxw .or. ny < nyw .or. nz < loaded_fields(index)%siz(3)) then
147  write(message1,'(i6,2i5)') nx,ny,nz
148  call mpp_error(fatal,'field '//trim(loaded_fields(index)%name)// &
149  ' Array size mismatch in time_interp_external. Array "data" is too small. shape(data)=' &
150  //message1)
151  endif
152  if(PRESENT(mask_out)) then
153  if (size(mask_out,1) /= nx .or. size(mask_out,2) /= ny .or. size(mask_out,3) /= nz) then
154  write(message1,'(i6,2i5)') nx,ny,nz
155  write(message2,'(i6,2i5)') size(mask_out,1),size(mask_out,2),size(mask_out,3)
156  call mpp_error(fatal,'field '//trim(loaded_fields(index)%name)// &
157  ' array size mismatch in time_interp_external.'// &
158  ' Shape of array "mask_out" does not match that of array "data".'// &
159  ' shape(data)='//message1//' shape(mask_out)='//message2)
160  endif
161  endif
162 
163  ! only one record in the file => time-independent field
164  if (loaded_fields(index)%siz(4) == 1) then
165  call load_record(loaded_fields(index),1,horz_interp, is_in, ie_in ,js_in, je_in,window_id)
166  i1 = find_buf_index(1,loaded_fields(index)%ibuf)
167  if( loaded_fields(index)%region_type == no_region ) then
168  where(loaded_fields(index)%mask(isc:iec,jsc:jec,:,i1))
169  time_data(isw:iew,jsw:jew,:) = real( loaded_fields(index)%domain_data(isc:iec,jsc:jec,:,i1),&
170  fms_ti_kind_)
171  elsewhere
172  ! time_data(isw:iew,jsw:jew,:) = time_interp_missing !field(index)%missing? Balaji
173  time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index)%missing, fms_ti_kind_)
174  end where
175  else
176  where(loaded_fields(index)%mask(isc:iec,jsc:jec,:,i1))
177  time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index)%domain_data(isc:iec,jsc:jec,:,i1),&
178  fms_ti_kind_)
179  end where
180  endif
181  if(PRESENT(mask_out)) mask_out(isw:iew,jsw:jew,:) = loaded_fields(index)%mask(isc:iec,jsc:jec,:,i1)
182  ! otherwise do interpolation
183  else
184  ! using time_interp_modulo
185  if(loaded_fields(index)%have_modulo_times) then
186  call time_interp(time,loaded_fields(index)%modulo_time_beg, loaded_fields(index)%modulo_time_end, &
187  loaded_fields(index)%time(:), w2, t1, t2, &
188  loaded_fields(index)%correct_leap_year_inconsistency, err_msg=err_msg)
189  if(err_msg .NE. '') then
190  filename = trim(loaded_fields(index)%fileobj%path)
191  call mpp_error(fatal,"time_interp_external 1: "//trim(err_msg)//&
192  ",file="//trim(filename)//",field="//trim(loaded_fields(index)%name) )
193  endif
194  ! using time_interp_list
195  else
196  if(loaded_fields(index)%modulo_time) then
197  mod_time=1
198  else
199  mod_time=0
200  endif
201  call time_interp(time,loaded_fields(index)%time(:),w2,t1,t2,modtime=mod_time, err_msg=err_msg)
202  if(err_msg .NE. '') then
203  filename = trim(loaded_fields(index)%fileobj%path)
204  call mpp_error(fatal,"time_interp_external 2: "//trim(err_msg)//&
205  ",file="//trim(filename)//",field="//trim(loaded_fields(index)%name) )
206  endif
207  endif
208  w1 = 1.0_kindl -w2
209  if (verb) then
210  call get_date(time,yy,mm,dd,hh,min,ss)
211  write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
212  'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss
213  write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
214  endif
215 
216  call load_record(loaded_fields(index),t1,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
217  call load_record(loaded_fields(index),t2,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
218  i1 = find_buf_index(t1,loaded_fields(index)%ibuf)
219  i2 = find_buf_index(t2,loaded_fields(index)%ibuf)
220  if(i1<0.or.i2<0) &
221  call mpp_error(fatal,'time_interp_external : records were not loaded correctly in memory')
222 
223  if (verb) then
224  write(outunit,*) 'ibuf= ',loaded_fields(index)%ibuf
225  write(outunit,*) 'i1,i2= ',i1, i2
226  endif
227 
228  if( loaded_fields(index)%region_type == no_region ) then
229  where(loaded_fields(index)%mask(isc:iec,jsc:jec,:,i1) .and. &
230  loaded_fields(index)%mask(isc:iec,jsc:jec,:,i2))
231  time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index)%domain_data(isc:iec,jsc:jec,:,i1), kindl)&
232  * w1 + real(loaded_fields(index)%domain_data(isc:iec,jsc:jec,:,i2), kindl) * w2
233  elsewhere
234  ! time_data(isw:iew,jsw:jew,:) = time_interp_missing !loaded_fields(index)%missing? Balaji
235  time_data(isw:iew,jsw:jew,:) = real(loaded_fields(index)%missing, kindl)
236  end where
237  else
238  where(loaded_fields(index)%mask(isc:iec,jsc:jec,:,i1) .and. &
239  loaded_fields(index)%mask(isc:iec,jsc:jec,:,i2))
240  time_data(isw:iew,jsw:jew,:) = real( loaded_fields(index)%domain_data(isc:iec,jsc:jec,:,i1), kindl)&
241  * w1 + real(loaded_fields(index)%domain_data(isc:iec,jsc:jec,:,i2), kindl) * w2
242  end where
243  endif
244  if(PRESENT(mask_out)) &
245  mask_out(isw:iew,jsw:jew,:) = &
246  loaded_fields(index)%mask(isc:iec,jsc:jec,:,i1).and.&
247  loaded_fields(index)%mask(isc:iec,jsc:jec,:,i2)
248  endif
249 
250  end subroutine time_interp_external_3d_
251 !</SUBROUTINE> NAME="time_interp_external"
252 
253  !> @brief Scalar interpolation for @ref time_interp_external
254  subroutine time_interp_external_0d_(index, time, time_data, verbose)
255 
256  integer, intent(in) :: index
257  type(time_type), intent(in) :: time
258  real(FMS_TI_KIND_), intent(inout) :: time_data
259  logical, intent(in), optional :: verbose
260 
261  integer :: t1, t2
262  integer :: i1, i2, mod_time
263  integer :: yy, mm, dd, hh, min, ss
264  character(len=256) :: err_msg
265  character(len=FMS_PATH_LEN) :: filename
266 
267  real(FMS_TI_KIND_) :: w1,w2
268  logical :: verb
269  integer, parameter :: kindl = fms_ti_kind_
270 
271  verb=.false.
272  if (PRESENT(verbose)) verb=verbose
273  if (debug_this_module) verb = .true.
274 
275  if (index < 1.or.index > num_fields) &
276  call mpp_error(fatal, &
277  & 'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize')
278 
279  if (loaded_fields(index)%siz(4) == 1) then
280  ! only one record in the file => time-independent loaded_fields
281  call load_record_0d(loaded_fields(index),1)
282  i1 = find_buf_index(1,loaded_fields(index)%ibuf)
283  time_data = real(loaded_fields(index)%domain_data(1,1,1,i1), fms_ti_kind_)
284  else
285  if(loaded_fields(index)%have_modulo_times) then
286  call time_interp(time,loaded_fields(index)%modulo_time_beg, loaded_fields(index)%modulo_time_end, &
287  loaded_fields(index)%time(:), w2, t1, t2, &
288  loaded_fields(index)%correct_leap_year_inconsistency, err_msg=err_msg)
289  if(err_msg .NE. '') then
290  filename = trim(loaded_fields(index)%fileobj%path)
291  call mpp_error(fatal,"time_interp_external 3:"//trim(err_msg)//&
292  ",file="//trim(filename)//",field="//trim(loaded_fields(index)%name) )
293  endif
294  else
295  if(loaded_fields(index)%modulo_time) then
296  mod_time=1
297  else
298  mod_time=0
299  endif
300  call time_interp(time,loaded_fields(index)%time(:),w2,t1,t2,modtime=mod_time, err_msg=err_msg)
301  if(err_msg .NE. '') then
302  filename = trim(loaded_fields(index)%fileobj%path)
303  call mpp_error(fatal,"time_interp_external 4:"//trim(err_msg)// &
304  ",file="//trim(filename)//",field="//trim(loaded_fields(index)%name) )
305  endif
306  endif
307  w1 = 1.0_kindl-w2
308  if (verb) then
309  call get_date(time,yy,mm,dd,hh,min,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,':',min,':',ss
312  write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
313  endif
314  call load_record_0d(loaded_fields(index),t1)
315  call load_record_0d(loaded_fields(index),t2)
316  i1 = find_buf_index(t1,loaded_fields(index)%ibuf)
317  i2 = find_buf_index(t2,loaded_fields(index)%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(index)%domain_data(1,1,1,i1), fms_ti_kind_)*w1 &
322  & + real(loaded_fields(index)%domain_data(1,1,1,i2), fms_ti_kind_)*w2
323  if (verb) then
324  write(outunit,*) 'ibuf= ',loaded_fields(index)%ibuf
325  write(outunit,*) 'i1,i2= ',i1, i2
326  endif
327  endif
328 
329  end subroutine time_interp_external_0d_
330 
331 ! ============================================================================
subroutine time_interp_external_2d_(index, 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
subroutine time_interp_external_3d_(index, 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_0d_(index, time, time_data, verbose)
Scalar interpolation for time_interp_external.
subroutine time_interp_external_1d_(index, time, data_in, interp, verbose, horz_interp, mask_out, is_in, ie_in, window_id)
1D interpolation for time_interp_external