FMS  2025.04
Flexible Modeling System
time_interp_external2.F90
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 !> @defgroup time_interp_external2_mod time_interp_external2_mod
19 !> @ingroup time_interp
20 !> @brief Perform I/O and time interpolation of external fields (contained in a file), using
21 !! fms2_io.
22 !!
23 !> @author M.J. Harrison
24 !!
25 !> Perform I/O and time interpolation for external fields.
26 !! Uses udunits library to calculate calendar dates and
27 !! convert units. Allows for reading data decomposed across
28 !! model horizontal grid using optional domain2d argument
29 !!
30 !! data are defined over data domain for domain2d data
31 !! (halo values are NOT updated by this module)
32 
33 !> @addtogroup time_interp_external2_mod
34 !> @{
35 module time_interp_external2_mod
36 
37 !<NAMELIST NAME="time_interp_external_nml">
38 ! <DATA NAME="num_io_buffers" TYPE="integer">
39 ! size of record dimension for internal buffer. This is useful for tuning i/o performance
40 ! particularly for large datasets (e.g. daily flux fields)
41 ! </DATA>
42 !</NAMELIST>
43 
44  use platform_mod, only : r8_kind, r4_kind
45  use fms_mod, only : write_version_number
46  use mpp_mod, only : mpp_error,fatal,warning,mpp_pe, stdout, stdlog, note
47  use mpp_mod, only : input_nml_file, mpp_npes, mpp_root_pe, mpp_broadcast, mpp_get_current_pelist
48  use time_manager_mod, only : time_type, get_date, set_date, operator ( >= ) , operator ( + ) , days_in_month, &
49  operator( - ), operator ( / ), operator ( // ) , days_in_year, increment_time, &
50  set_time, get_time, operator( > ), get_calendar_type, no_calendar
51  use get_cal_time_mod, only : get_cal_time
52  use mpp_domains_mod, only : domain2d, mpp_get_compute_domain, mpp_get_data_domain, &
53  mpp_get_global_domain, null_domain2d
54  use time_interp_mod, only : time_interp, time_interp_init
55  use axis_utils2_mod, only : get_axis_cart, get_axis_modulo, get_axis_modulo_times
56  use fms_mod, only : lowercase, check_nml_error
57  use platform_mod, only: r8_kind, fms_path_len, fms_file_len
58  use horiz_interp_mod, only : horiz_interp, horiz_interp_type
59  use horiz_interp_type_mod, only : conserve
60  use fms2_io_mod, only : valid_t, fmsnetcdfdomainfile_t, open_file, get_unlimited_dimension_name, &
61  variable_att_exists, fmsnetcdffile_t, &
62  variable_exists, get_valid, get_variable_num_dimensions, read_data, &
63  is_valid, close_file, get_dimension_size, get_variable_dimension_names, &
64  get_variable_size, get_time_calendar, get_variable_missing, get_variable_units
65 
66  implicit none
67  private
68 
69 ! Include variable "version" to be written to log file.
70 #include<file_version.h>
71 
72  integer, parameter, public :: NO_REGION=0, inside_region=1, outside_region=2
73  integer, parameter, private :: modulo_year= 0001
74  integer, parameter, private :: LINEAR_TIME_INTERP = 1 ! not used currently
75  integer, parameter, public :: SUCCESS = 0, err_field_not_found = 1
76  real(r8_kind), parameter, private :: DEFAULT_MISSING_VALUE = -1e20_r8_kind
77  integer, private :: max_fields = 100, max_files= 40
78  integer, private :: num_fields = 0, num_files=0
79  ! denotes time intervals in file (interpreted from metadata)
80  integer, private :: num_io_buffers = 2 ! set -1 to read all records from disk into memory
81  logical, private :: module_initialized = .false.
82  logical, private :: debug_this_module = .false.
83 
86  public set_override_region, reset_src_data_region
87  public get_external_fileobj
89 
90  private find_buf_index,&
91  set_time_modulo
92  !> @}
93 
94  !> Represents external fields
95  !> @ingroup time_interp_external2_mod
96  type, private :: ext_fieldtype
97  type(fmsnetcdffile_t), pointer :: fileobj=>null() !< keep unit open when not reading all records
98  character(len=128) :: name, units
99  integer :: siz(4), ndim
100  character(len=32) :: axisname(4)
101  type(domain2d) :: domain
102  type(time_type), dimension(:), pointer :: time =>null() !< midpoint of time interval
103  type(time_type), dimension(:), pointer :: start_time =>null(), end_time =>null()
104  type(time_type), dimension(:), pointer :: period =>null()
105  logical :: modulo_time !< denote climatological time axis
106  real(r8_kind), dimension(:,:,:,:), pointer :: domain_data =>null() !< defined over data domain or global domain
107  logical, dimension(:,:,:,:), pointer :: mask =>null() !< defined over data domain or global domain
108  integer, dimension(:), pointer :: ibuf =>null() !< record numbers associated with buffers
109  real(r8_kind), dimension(:,:,:,:), pointer :: src_data =>null() !< input data buffer
110  type(valid_t) :: valid ! data validator
111  integer :: nbuf
112  logical :: domain_present
113  real(r8_kind) :: slope, intercept
114  integer :: isc,iec,jsc,jec
115  type(time_type) :: modulo_time_beg, modulo_time_end
116  logical :: have_modulo_times, correct_leap_year_inconsistency
117  integer :: region_type
118  integer :: is_region, ie_region, js_region, je_region
119  integer :: is_src, ie_src, js_src, je_src
120  integer :: tdim
121  integer :: numwindows
122  logical, dimension(:,:), pointer :: need_compute=>null()
123  real(r8_kind) :: missing ! missing value
124  end type ext_fieldtype
125 
126  !> Holds filename and file object
127  !> @ingroup time_interp_external2_mod
128  type, private :: filetype
129  character(len=FMS_FILE_LEN) :: filename = ''
130  type(FmsNetcdfFile_t), pointer :: fileobj => null()
131  end type filetype
132 
133  !> Provide data from external file interpolated to current model time.
134  !! Data may be local to current processor or global, depending on
135  !! "init_external_field" flags. Uses @ref fms2_io for I/O.
136  !!
137  !! @param index index of external field from previous call to init_external_field
138  !! @param time target time for data
139  !! @param [inout] data global or local data array
140  !! @param interp time_interp_external defined interpolation method (optional). Currently
141  !! this module only supports LINEAR_TIME_INTERP.
142  !! @param verbose verbose flag for debugging (optional).
143  !!
144  !> @ingroup time_interp_external2_mod
146  module procedure time_interp_external_0d_r4
147  module procedure time_interp_external_1d_r4
148  module procedure time_interp_external_2d_r4
149  module procedure time_interp_external_3d_r4
150  module procedure time_interp_external_0d_r8
151  module procedure time_interp_external_1d_r8
152  module procedure time_interp_external_2d_r8
153  module procedure time_interp_external_3d_r8
154  end interface
155 
157  module procedure time_interp_external_bridge_0d_r4
158  module procedure time_interp_external_bridge_1d_r4
159  module procedure time_interp_external_bridge_2d_r4
160  module procedure time_interp_external_bridge_3d_r4
161  module procedure time_interp_external_bridge_0d_r8
162  module procedure time_interp_external_bridge_1d_r8
163  module procedure time_interp_external_bridge_2d_r8
164  module procedure time_interp_external_bridge_3d_r8
165  end interface
166 
167  !> @addtogroup time_interp_external2_mod
168  !> @{
169 
170  integer :: outunit
171 
172  type(ext_fieldtype), save, private, pointer :: loaded_fields(:) => null()
173  type(filetype), save, private, pointer :: opened_files(:) => null()
174 !Balaji: really should use field%missing
175  real(r8_kind), private, parameter :: time_interp_missing=-1e99_r8_kind
176  contains
177 
178 ! <SUBROUTINE NAME="time_interp_external_init">
179 !
180 ! <DESCRIPTION>
181 ! Initialize the time_interp_external module
182 ! </DESCRIPTION>
183 !
184  !> @brief Initialize the @ref time_interp_external_mod module
186 
187  integer :: io_status, logunit, ierr
188 
189  namelist /time_interp_external_nml/ num_io_buffers, debug_this_module, &
190  max_fields, max_files
191 
192  ! open and read namelist
193 
194  if(module_initialized) return
195 
196  logunit = stdlog()
197  outunit = stdout()
198  call write_version_number("TIME_INTERP_EXTERNAL_MOD", version)
199 
200  read (input_nml_file, time_interp_external_nml, iostat=io_status)
201  ierr = check_nml_error(io_status, 'time_interp_external_nml')
202 
203  write(logunit,time_interp_external_nml)
204  call realloc_fields(max_fields)
205  call realloc_files(max_files)
206 
207  module_initialized = .true.
208 
209  call time_interp_init()
210 
211  return
212 
213  end subroutine time_interp_external_init
214 
215 
216 !<FUNCTION NAME="init_external_field" TYPE="integer">
217 !
218 !<DESCRIPTION>
219 ! initialize an external field. Buffer "num_io_buffers" (default=2) in memory to reduce memory allocations.
220 ! distributed reads are supported using the optional "domain" flag.
221 ! Units conversion via the optional "desired_units" flag using udunits_mod.
222 !
223 ! Return integer id of field for future calls to time_interp_external.
224 !
225 ! </DESCRIPTION>
226 !
227 !<IN NAME="file" TYPE="character(len=*)">
228 ! filename
229 !</IN>
230 !<IN NAME="fieldname" TYPE="character(len=*)">
231 ! fieldname (in file)
232 !</IN>
233 !<IN NAME="format" TYPE="integer">
234 ! mpp_io flag for format of file (optional). Currently only "MPP_NETCDF" supported
235 !</IN>
236 !<IN NAME="threading" TYPE="integer">
237 ! mpp_io flag for threading (optional). "MPP_SINGLE" means root pe reads global field and distributes to other PEs
238 ! "MPP_MULTI" means all PEs read data
239 !</IN>
240 !<IN NAME="domain" TYPE="mpp_domains_mod:domain2d">
241 ! domain flag (optional)
242 !</IN>
243 !<IN NAME="desired_units" TYPE="character(len=*)">
244 ! Target units for data (optional), e.g. convert from deg_K to deg_C.
245 ! Failure to convert using udunits will result in failure of this module.
246 !</IN>
247 !<IN NAME="verbose" TYPE="logical">
248 ! verbose flag for debugging (optional).
249 !</IN>
250 !<INOUT NAME="axis_centers" TYPE="axistype" DIM="(4)">
251 ! MPP_IO axistype array for grid centers ordered X-Y-Z-T (optional).
252 !</INOUT>
253 !<INOUT NAME="axis_sizes" TYPE="integer" DIM="(4)">
254 ! array of axis lengths ordered X-Y-Z-T (optional).
255 !</INOUT>
256 
257 
258  !> Initialize an external field. Buffer "num_io_buffers" (default=2) in memory to reduce memory allocations.
259  !! distributed reads are supported using the optional "domain" flag.
260  !! Units conversion via the optional "desired_units" flag using udunits_mod.
261  !!
262  !> @return integer id of field for future calls to time_interp_external.
263  !> @param file filename
264  !> @param fieldname fieldname (in file)
265  !> @param format mpp_io flag for format of file(optional). Currently only "MPP_NETCDF" supported
266  !> @param threading mpp_io flag for threading (optional). "MPP_SINGLE" means root pe reads
267  !! global field and distributes to other PEs. "MPP_MULTI" means all PEs read data
268  !> @param domain domain flag (optional)
269  !> @param desired_units Target units for data (optional), e.g. convert from deg_K to deg_C.
270  !! Failure to convert using udunits will result in failure of this module.
271  !> @param verbose verbose flag for debugging (optional).
272  !> @param [out] axis_names List of axis names (optional).
273  !> @param [inout] axis_sizes array of axis lengths ordered X-Y-Z-T (optional).
274  function init_external_field(file,fieldname,domain,desired_units,&
275  verbose,axis_names, axis_sizes,override,correct_leap_year_inconsistency,&
276  permit_calendar_conversion,use_comp_domain,ierr, nwindows, ignore_axis_atts, ongrid )
277 
278  character(len=*), intent(in) :: file,fieldname
279  logical, intent(in), optional :: verbose
280  character(len=*), intent(in), optional :: desired_units
281  type(domain2d), intent(in), optional :: domain
282  integer, intent(inout), optional :: axis_sizes(4)
283  character(len=*), intent(out), optional :: axis_names(4)
284  logical, intent(in), optional :: override, correct_leap_year_inconsistency,&
285  permit_calendar_conversion,use_comp_domain
286  integer, intent(out), optional :: ierr
287  integer, intent(in), optional :: nwindows
288  logical, optional :: ignore_axis_atts
289  logical, optional :: ongrid !< Optional flag indicating if the data is ongrid
290 
291  logical :: ongrid_local !< Flag indicating if the data is ongrid
292 
293  integer :: init_external_field
294 
295  real(r8_kind) :: slope, intercept
296  integer :: ndim,ntime,i,j
297  integer :: iscomp,iecomp,jscomp,jecomp,isglobal,ieglobal,jsglobal,jeglobal
298  integer :: isdata,iedata,jsdata,jedata, dxsize, dysize,dxsize_max,dysize_max
299  logical :: verb, transpose_xy,use_comp_domain1
300  real(r8_kind), dimension(:), allocatable :: tstamp, tstart, tend, tavg
301  character(len=1) :: cart
302  character(len=1), dimension(4) :: cart_dir
303  character(len=128) :: units, fld_units
304  character(len=128) :: msg, calendar_type, timebeg, timeend
305  character(len=128) :: timename, timeunits
306  character(len=128), allocatable :: axisname(:)
307  integer, allocatable :: axislen(:)
308  integer :: siz(4), siz_in(4), gxsize, gysize,gxsize_max, gysize_max
309  type(time_type) :: tdiff
310  integer :: yr, mon, day, hr, minu, sec
311  integer :: len, nfile, nfields_orig, nbuf, nx,ny
312  integer :: numwindows
313  logical :: ignore_axatts
314  logical :: have_modulo_time
315  type(fmsnetcdffile_t), pointer :: fileobj=>null()
316  integer, dimension(:), allocatable :: pes !< List of ranks in the current pelist
317 
318  if (.not. module_initialized) call mpp_error(fatal,'Must call time_interp_external_init first')
319  if(present(ierr)) ierr = success
320  ignore_axatts=.false.
321  cart_dir(1)='X';cart_dir(2)='Y';cart_dir(3)='Z';cart_dir(4)='T'
322  if(present(ignore_axis_atts)) ignore_axatts = ignore_axis_atts
323  use_comp_domain1 = .false.
324  if(PRESENT(use_comp_domain)) use_comp_domain1 = use_comp_domain
325  verb=.false.
326  if (PRESENT(verbose)) verb=verbose
327  if (debug_this_module) verb = .true.
328  numwindows = 1
329  if(present(nwindows)) numwindows = nwindows
330 
331  units = 'same'
332  if (PRESENT(desired_units)) then
333  units = desired_units
334  call mpp_error(fatal,'==> Unit conversion via time_interp_external &
335  &has been temporarily deprecated. Previous versions of&
336  &this module used udunits_mod to perform unit conversion.&
337  & Udunits_mod is in the process of being replaced since &
338  &there were portability issues associated with this code.&
339  & Please remove the desired_units argument from calls to &
340  &this routine.')
341  endif
342  nfile = 0
343  do i=1,num_files
344  if(trim(opened_files(i)%filename) == trim(file)) then
345  nfile = i
346  exit ! file is already opened
347  endif
348  enddo
349  if(nfile == 0) then
350  num_files = num_files + 1
351  if(num_files > max_files) then ! not enough space in the file table, reallocate it
352  !--- z1l: For the case of multiple thread, realoc_files will cause memory leak.
353  !--- If multiple threads are working on file A. One of the thread finished first and
354  !--- begin to work on file B, the realloc_files will cause problem for
355  !--- other threads are working on the file A.
356  ! call realloc_files(2*size(opened_files))
357  call mpp_error(fatal, "time_interp_external: num_files is greater than max_files, "// &
358  "increase time_interp_external_nml max_files")
359  endif
360  opened_files(num_files)%filename = trim(file)
361  allocate(opened_files(num_files)%fileobj)
362  fileobj => opened_files(num_files)%fileobj
363 
364  if(.not. open_file(opened_files(num_files)%fileobj, trim(file), 'read')) &
365  call mpp_error(fatal, 'time_interp_external_mod: Error in opening file '//trim(file))
366  else
367  fileobj => opened_files(nfile)%fileobj
368  endif
369 
370  call get_unlimited_dimension_name(fileobj, timename)
371  call get_dimension_size(fileobj, timename, ntime)
372 
373  if (ntime < 1) then
374  write(msg,'(a15,a,a58)') 'external field ',trim(fieldname),&
375  ' does not have an associated record dimension (REQUIRED) '
376  call mpp_error(fatal,trim(msg))
377  endif
378 
379  !--- get time calendar_type
380  call get_time_calendar(fileobj, timename, calendar_type)
381 
382  !--- get time units
383  call get_variable_units(fileobj, timename, timeunits)
384 
385  !--- get timebeg and timeend
386  have_modulo_time = get_axis_modulo_times(fileobj, timename, timebeg, timeend)
387 
388  allocate(pes(mpp_npes()))
389  call mpp_get_current_pelist(pes)
390  allocate(tstamp(ntime),tstart(ntime),tend(ntime),tavg(ntime))
391 
392  !< Only root reads the unlimited dimension and broadcasts it to the other ranks
393  if (mpp_root_pe() .eq. mpp_pe()) call read_data(fileobj, timename, tstamp)
394  call mpp_broadcast(tstamp, size(tstamp), mpp_root_pe(), pelist=pes)
395  deallocate(pes)
396 
397  transpose_xy = .false.
398  isdata=1; iedata=1; jsdata=1; jedata=1
399  gxsize=1; gysize=1
400  siz_in = 1
401 
402  if (PRESENT(domain)) then
403  call mpp_get_compute_domain(domain,iscomp,iecomp,jscomp,jecomp)
404  nx = iecomp-iscomp+1; ny = jecomp-jscomp+1
405  call mpp_get_data_domain(domain,isdata,iedata,jsdata,jedata,dxsize,dxsize_max,dysize,dysize_max)
406  call mpp_get_global_domain(domain,isglobal,ieglobal,jsglobal,jeglobal,gxsize,gxsize_max,gysize,gysize_max)
407  ongrid_local = .false.
408  if (present(ongrid)) ongrid_local = ongrid
409  !> If this is an ongrid case, set is[e]js[e]data to be equal to the compute domain.
410  !! This is what it is used to allocate space for the data!
411  if (ongrid_local) then
412  isdata=iscomp
413  iedata=iecomp
414  jsdata=jscomp
415  jedata=jecomp
416  endif
417  elseif(use_comp_domain1) then
418  call mpp_error(fatal,"init_external_field:"//&
419  " use_comp_domain=true but domain is not present")
420  endif
421 
423  nfields_orig = num_fields
424 
425  if (.not. variable_exists(fileobj, fieldname) ) then
426  if (present(ierr)) then
427  ierr = err_field_not_found
428  return
429  else
430  call mpp_error(fatal,'external field "'//trim(fieldname)//'" not found in file "'//trim(file)//'"')
431  endif
432  endif
433 
434  tavg = -1.0_r8_kind
435  tstart = tstamp
436  tend = tstamp
437  if(variable_att_exists(fileobj, fieldname, 'time_avg_info')) then
438  if(variable_exists(fileobj, 'average_T1')) call read_data(fileobj, 'average_T1', tstart)
439  if(variable_exists(fileobj, 'average_T2')) call read_data(fileobj, 'average_T2', tend)
440  if(variable_exists(fileobj, 'average_DT')) call read_data(fileobj, 'average_DT', tavg)
441  endif
442 
443  num_fields = num_fields + 1
444  if(num_fields > max_fields) then
445  !--- z1l: For the case of multiple thread, realoc_fields will cause memory leak.
446  !--- If multiple threads are working on field A. One of the thread finished first and
447  !--- begin to work on field B, the realloc_files will cause problem for
448  !--- other threads are working on the field A.
449  !call realloc_fields(size(field)*2)
450  call mpp_error(fatal, "time_interp_external: num_fields is greater than max_fields, "// &
451  "increase time_interp_external_nml max_fields")
452  endif
453 
454  !--- get field attribute
455  call get_variable_units(fileobj, fieldname, fld_units)
456 
457  init_external_field = num_fields
458  loaded_fields(num_fields)%fileobj => fileobj
459  loaded_fields(num_fields)%name = trim(fieldname)
460  loaded_fields(num_fields)%units = trim(fld_units)
461  loaded_fields(num_fields)%isc = 1
462  loaded_fields(num_fields)%iec = 1
463  loaded_fields(num_fields)%jsc = 1
464  loaded_fields(num_fields)%jec = 1
465  loaded_fields(num_fields)%region_type = no_region
466  loaded_fields(num_fields)%is_region = 0
467  loaded_fields(num_fields)%ie_region = -1
468  loaded_fields(num_fields)%js_region = 0
469  loaded_fields(num_fields)%je_region = -1
470  if (PRESENT(domain)) then
471  loaded_fields(num_fields)%domain_present = .true.
472  loaded_fields(num_fields)%domain = domain
473  loaded_fields(num_fields)%isc=iscomp;loaded_fields(num_fields)%iec = iecomp
474  loaded_fields(num_fields)%jsc=jscomp;loaded_fields(num_fields)%jec = jecomp
475  else
476  loaded_fields(num_fields)%domain_present = .false.
477  endif
478 
479  loaded_fields(num_fields)%valid = get_valid(fileobj, fieldname)
480  ndim = get_variable_num_dimensions(fileobj, fieldname)
481  if (ndim > 4) call mpp_error(fatal, &
482  'invalid array rank <=4d fields supported')
483 
484  loaded_fields(num_fields)%siz = 1
485  loaded_fields(num_fields)%ndim = ndim
486  loaded_fields(num_fields)%tdim = 4
487  !--- get field missing value
488  loaded_fields(num_fields)%missing = get_variable_missing(fileobj, fieldname)
489 
490  allocate(axisname(ndim), axislen(ndim))
491 
492  call get_variable_dimension_names(fileobj, fieldname, axisname)
493  call get_variable_size(fileobj, fieldname, axislen)
494  do j=1,loaded_fields(num_fields)%ndim
495  call get_axis_cart(fileobj, axisname(j), cart)
496  len = axislen(j)
497  if (cart == 'N' .and. .not. ignore_axatts) then
498  write(msg,'(a,"/",a)') trim(file),trim(fieldname)
499  call mpp_error(fatal,'file/field '//trim(msg)// &
500  ' couldnt recognize axis atts in time_interp_external')
501  else if (cart == 'N' .and. ignore_axatts) then
502  cart = cart_dir(j)
503  endif
504  select case (cart)
505  case ('X')
506  if (j.eq.2) transpose_xy = .true.
507  if (.not.PRESENT(domain) .and. .not.PRESENT(override)) then
508  isdata=1;iedata=len
509  iscomp=1;iecomp=len
510  gxsize = len
511  dxsize = len
512  loaded_fields(num_fields)%isc=iscomp;loaded_fields(num_fields)%iec=iecomp
513  elseif (PRESENT(override)) then
514  gxsize = len
515  if (PRESENT(axis_sizes)) axis_sizes(1) = len
516  endif
517  loaded_fields(num_fields)%axisname(1) = axisname(j)
518  if(use_comp_domain1) then
519  loaded_fields(num_fields)%siz(1) = nx
520  else
521  loaded_fields(num_fields)%siz(1) = dxsize
522  endif
523  if (len /= gxsize) then
524  write(msg,'(a,"/",a)') trim(file),trim(fieldname)
525  call mpp_error(fatal,'time_interp_ext, file/field '//trim(msg)//' x dim doesnt match model')
526  endif
527  case ('Y')
528  loaded_fields(num_fields)%axisname(2) = axisname(j)
529  if (.not.PRESENT(domain) .and. .not.PRESENT(override)) then
530  jsdata=1;jedata=len
531  jscomp=1;jecomp=len
532  gysize = len
533  dysize = len
534  loaded_fields(num_fields)%jsc=jscomp;loaded_fields(num_fields)%jec=jecomp
535  elseif (PRESENT(override)) then
536  gysize = len
537  if (PRESENT(axis_sizes)) axis_sizes(2) = len
538  endif
539  if(use_comp_domain1) then
540  loaded_fields(num_fields)%siz(2) = ny
541  else
542  loaded_fields(num_fields)%siz(2) = dysize
543  endif
544  if (len /= gysize) then
545  write(msg,'(a,"/",a)') trim(file),trim(fieldname)
546  call mpp_error(fatal,'time_interp_ext, file/field '//trim(msg)//' y dim doesnt match model')
547  endif
548  case ('Z')
549  loaded_fields(num_fields)%axisname(3) = axisname(j)
550  loaded_fields(num_fields)%siz(3) = len
551  case ('T')
552  loaded_fields(num_fields)%axisname(4) = axisname(j)
553  loaded_fields(num_fields)%siz(4) = ntime
554  loaded_fields(num_fields)%tdim = j
555  end select
556  enddo
557  siz = loaded_fields(num_fields)%siz
558  if(PRESENT(axis_names)) axis_names = loaded_fields(num_fields)%axisname
559  if (PRESENT(axis_sizes) .and. .not.PRESENT(override)) then
560  axis_sizes = loaded_fields(num_fields)%siz
561  endif
562 
563  if (verb) write(outunit,'(a,4i6)') 'field x,y,z,t local size= ',siz
564  if (verb) write(outunit,*) 'field contains data in units = ',trim(loaded_fields(num_fields)%units)
565  if (transpose_xy) call mpp_error(fatal,'axis ordering not supported')
566  if (num_io_buffers .le. 1) call mpp_error(fatal,'time_interp_ext:num_io_buffers should be at least 2')
567  nbuf = min(num_io_buffers,siz(4))
568 
569  loaded_fields(num_fields)%numwindows = numwindows
570  allocate(loaded_fields(num_fields)%need_compute(nbuf, numwindows))
571  loaded_fields(num_fields)%need_compute = .true.
572 
573  allocate(loaded_fields(num_fields)%domain_data(isdata:iedata,jsdata:jedata,siz(3),nbuf),&
574  loaded_fields(num_fields)%mask(isdata:iedata,jsdata:jedata,siz(3),nbuf) )
575  loaded_fields(num_fields)%mask = .false.
576  loaded_fields(num_fields)%domain_data = 0.0_r8_kind
577  slope=1.0_r8_kind;intercept=0.0_r8_kind
578 ! if (units /= 'same') call convert_units(trim(field(num_fields)%units),trim(units),slope,intercept)
579 ! if (verb.and.units /= 'same') then
580 ! write(outunit,*) 'attempting to convert data to units = ',trim(units)
581 ! write(outunit,'(a,f8.3,a,f8.3)') 'factor = ',slope,' offset= ',intercept
582 ! endif
583  loaded_fields(num_fields)%slope = slope
584  loaded_fields(num_fields)%intercept = intercept
585  allocate(loaded_fields(num_fields)%ibuf(nbuf))
586  loaded_fields(num_fields)%ibuf = -1
587  loaded_fields(num_fields)%nbuf = 0 ! initialize buffer number so that first reading fills data(:,:,:,1)
588  if(PRESENT(override)) then
589  loaded_fields(num_fields)%is_src = 1
590  loaded_fields(num_fields)%ie_src = gxsize
591  loaded_fields(num_fields)%js_src = 1
592  loaded_fields(num_fields)%je_src = gysize
593  allocate(loaded_fields(num_fields)%src_data(gxsize,gysize,siz(3),nbuf))
594  else
595  loaded_fields(num_fields)%is_src = isdata
596  loaded_fields(num_fields)%ie_src = iedata
597  loaded_fields(num_fields)%js_src = jsdata
598  loaded_fields(num_fields)%je_src = jedata
599  allocate(loaded_fields(num_fields)%src_data(isdata:iedata,jsdata:jedata,siz(3),nbuf))
600  endif
601 
602  allocate(loaded_fields(num_fields)%time(ntime))
603  allocate(loaded_fields(num_fields)%period(ntime))
604  allocate(loaded_fields(num_fields)%start_time(ntime))
605  allocate(loaded_fields(num_fields)%end_time(ntime))
606 
607  do j=1,ntime
608  loaded_fields(num_fields)%time(j) = get_cal_time(tstamp(j),trim(timeunits),trim(calendar_type), &
609  & permit_calendar_conversion)
610  loaded_fields(num_fields)%start_time(j) = get_cal_time(tstart(j),trim(timeunits),trim(calendar_type), &
611  & permit_calendar_conversion)
612  loaded_fields(num_fields)%end_time(j) = get_cal_time( tend(j),trim(timeunits),trim(calendar_type), &
613  & permit_calendar_conversion)
614  enddo
615 
616  if (loaded_fields(num_fields)%modulo_time) then
617  call set_time_modulo(loaded_fields(num_fields)%Time)
618  call set_time_modulo(loaded_fields(num_fields)%start_time)
619  call set_time_modulo(loaded_fields(num_fields)%end_time)
620  endif
621 
622  if(present(correct_leap_year_inconsistency)) then
623  loaded_fields(num_fields)%correct_leap_year_inconsistency = correct_leap_year_inconsistency
624  else
625  loaded_fields(num_fields)%correct_leap_year_inconsistency = .false.
626  endif
627 
628  if(have_modulo_time) then
629  if(get_calendar_type() == no_calendar) then
630  loaded_fields(num_fields)%modulo_time_beg = set_time(timebeg)
631  loaded_fields(num_fields)%modulo_time_end = set_time(timeend)
632  else
633  loaded_fields(num_fields)%modulo_time_beg = set_date(timebeg)
634  loaded_fields(num_fields)%modulo_time_end = set_date(timeend)
635  endif
636  loaded_fields(num_fields)%have_modulo_times = .true.
637  else
638  loaded_fields(num_fields)%have_modulo_times = .false.
639  endif
640  if(ntime == 1) then
641  call mpp_error(note, 'time_interp_external_mod: file '//trim(file)//' has only one time level')
642  else
643  do j= 1, ntime
644  loaded_fields(num_fields)%period(j) = loaded_fields(num_fields)%end_time(j) &
645  - loaded_fields(num_fields)%start_time(j)
646  if (loaded_fields(num_fields)%period(j) > set_time(0,0)) then
647  call get_time(loaded_fields(num_fields)%period(j), sec, day)
648  sec = sec/2+mod(day,2)*43200
649  day = day/2
650  loaded_fields(num_fields)%time(j) = loaded_fields(num_fields)%start_time(j)+&
651  set_time(sec,day)
652  else
653  if (j > 1 .and. j < ntime) then
654  tdiff = loaded_fields(num_fields)%time(j+1) - loaded_fields(num_fields)%time(j-1)
655  call get_time(tdiff, sec, day)
656  sec = sec/2+mod(day,2)*43200
657  day = day/2
658  loaded_fields(num_fields)%period(j) = set_time(sec,day)
659  sec = sec/2+mod(day,2)*43200
660  day = day/2
661  loaded_fields(num_fields)%start_time(j) = loaded_fields(num_fields)%time(j) - set_time(sec,day)
662  loaded_fields(num_fields)%end_time(j) = loaded_fields(num_fields)%time(j) + set_time(sec,day)
663  elseif ( j == 1) then
664  tdiff = loaded_fields(num_fields)%time(2) - loaded_fields(num_fields)%time(1)
665  call get_time(tdiff, sec, day)
666  loaded_fields(num_fields)%period(j) = set_time(sec,day)
667  sec = sec/2+mod(day,2)*43200
668  day = day/2
669  loaded_fields(num_fields)%start_time(j) = loaded_fields(num_fields)%time(j) - set_time(sec,day)
670  loaded_fields(num_fields)%end_time(j) = loaded_fields(num_fields)%time(j) + set_time(sec,day)
671  else
672  tdiff = loaded_fields(num_fields)%time(ntime) - loaded_fields(num_fields)%time(ntime-1)
673  call get_time(tdiff, sec, day)
674  loaded_fields(num_fields)%period(j) = set_time(sec,day)
675  sec = sec/2+mod(day,2)*43200
676  day = day/2
677  loaded_fields(num_fields)%start_time(j) = loaded_fields(num_fields)%time(j) - set_time(sec,day)
678  loaded_fields(num_fields)%end_time(j) = loaded_fields(num_fields)%time(j) + set_time(sec,day)
679  endif
680  endif
681  enddo
682  endif
683 
684  do j=1,ntime-1
685  if (loaded_fields(num_fields)%time(j) >= loaded_fields(num_fields)%time(j+1)) then
686  write(msg,'(A,i20)') "times not monotonically increasing. Filename: " &
687  //trim(file)//" field: "//trim(fieldname)//" timeslice: ", j
688  call mpp_error(fatal, trim(msg))
689  endif
690  enddo
691 
692  loaded_fields(num_fields)%modulo_time = get_axis_modulo(fileobj, timename)
693 
694  if (verb) then
695  if (loaded_fields(num_fields)%modulo_time) write(outunit,*) 'data are being treated as modulo in time'
696  do j= 1, ntime
697  write(outunit,*) 'time index, ', j
698  call get_date(loaded_fields(num_fields)%start_time(j),yr,mon,day,hr,minu,sec)
699  write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
700  'start time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
701  call get_date(loaded_fields(num_fields)%time(j),yr,mon,day,hr,minu,sec)
702  write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
703  'mid time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
704  call get_date(loaded_fields(num_fields)%end_time(j),yr,mon,day,hr,minu,sec)
705  write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
706  'end time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
707  enddo
708  end if
709 
710  deallocate(axisname, axislen)
711  deallocate(tstamp, tstart, tend, tavg)
712 
713  return
714 
715  end function init_external_field
716 
717 
718  function get_external_fileobj(filename, fileobj)
719  character(len=*), intent(in) :: filename
720  type(fmsnetcdffile_t), intent(out) :: fileobj
721  logical :: get_external_fileobj
722  integer :: i
723 
724  get_external_fileobj = .false.
725  do i=1,num_files
726  if(trim(opened_files(i)%filename) == trim(filename)) then
727  fileobj = opened_files(i)%fileobj
728  get_external_fileobj = .true.
729  exit ! file is already opened
730  endif
731  enddo
732 
733  end function get_external_fileobj
734 
735 
736  subroutine set_time_modulo(Time)
737 
738  type(time_type), intent(inout), dimension(:) :: time
739 
740  integer :: ntime, n
741  integer :: yr, mon, dy, hr, minu, sec
742 
743  ntime = size(time(:))
744 
745  do n = 1, ntime
746  call get_date(time(n), yr, mon, dy, hr, minu, sec)
747  yr = modulo_year
748  time(n) = set_date(yr, mon, dy, hr, minu, sec)
749  enddo
750 
751 
752  end subroutine set_time_modulo
753 
754 ! ============================================================================
755 !> load specified record from file
756 !> Always loads in r8, casts down for horiz_interp if interp argument is already allocated for r4's.
757 subroutine load_record(field, rec, interp, is_in, ie_in, js_in, je_in, window_id_in)
758  type(ext_fieldtype), intent(inout) :: field
759  integer , intent(in) :: rec ! record number
760  type(horiz_interp_type), intent(in), optional :: interp
761  integer, intent(in), optional :: is_in, ie_in, js_in, je_in
762  integer, intent(in), optional :: window_id_in
763 
764  ! ---- local vars
765  integer :: ib ! index in the array of input buffers
766  integer :: isw,iew,jsw,jew ! boundaries of the domain on each window
767  integer :: is_region, ie_region, js_region, je_region, i, j
768  integer :: start(4), nread(4)
769  logical :: need_compute
770  integer :: window_id
771  real(r8_kind) :: mask_in(size(field%src_data,1),size(field%src_data,2),size(field%src_data,3))
772  real(r8_kind), allocatable :: mask_out(:,:,:)
773  real(r4_kind), allocatable :: hi_tmp_data(:,:,:,:) !< used to hold a copy of field%domain_data if using r4_kind
774  real(r4_kind), allocatable :: hi_tmp_msk_out(:,:,:) !< used return the field mask if using r4_kind
775 
776  window_id = 1
777  if( PRESENT(window_id_in) ) window_id = window_id_in
778  need_compute = .true.
779 
780 !$OMP CRITICAL
781  ib = find_buf_index(rec,field%ibuf)
782 
783  if(ib>0) then
784  !--- do nothing
785  need_compute = .false.
786  else
787  ! calculate current buffer number in round-robin fasion
788  field%nbuf = field%nbuf + 1
789  if(field%nbuf > size(field%domain_data,4).or.field%nbuf <= 0) field%nbuf = 1
790  ib = field%nbuf
791  field%ibuf(ib) = rec
792  field%need_compute(ib,:) = .true.
793 
794  if (debug_this_module) write(outunit,*) 'reading record without domain for field ',trim(field%name)
795  start = 1; nread = 1
796  start(1) = field%is_src; nread(1) = field%ie_src - field%is_src + 1
797  start(2) = field%js_src; nread(2) = field%je_src - field%js_src + 1
798  start(3) = 1; nread(3) = size(field%src_data,3)
799  start(field%tdim) = rec; nread(field%tdim) = 1
800  call read_data(field%fileobj,field%name,field%src_data(:,:,:,ib:ib),corner=start,edge_lengths=nread)
801  endif
802 !$OMP END CRITICAL
803  isw=field%isc;iew=field%iec
804  jsw=field%jsc;jew=field%jec
805 
806  if( field%numwindows > 1) then
807  if( .NOT. PRESENT(is_in) .OR. .NOT. PRESENT(ie_in) .OR. .NOT. PRESENT(js_in) .OR. .NOT. PRESENT(je_in) ) then
808  call mpp_error(fatal, &
809  & 'time_interp_external(load_record): is_in, ie_in, js_in, je_in must be present when numwindows>1')
810  endif
811  isw = isw + is_in - 1
812  iew = isw + ie_in - is_in
813  jsw = jsw + js_in - 1
814  jew = jsw + je_in - js_in
815  endif
816 
817  ! interpolate to target grid
818 
819  need_compute = field%need_compute(ib, window_id)
820  if(need_compute) then
821  if(PRESENT(interp)) then
822  is_region = field%is_region; ie_region = field%ie_region
823  js_region = field%js_region; je_region = field%je_region
824  mask_in = 0.0_r8_kind
825  where (is_valid(field%src_data(:,:,:,ib), field%valid)) mask_in = 1.0_r8_kind
826  if ( field%region_type .NE. no_region ) then
827  if( any(mask_in == 0.0_r8_kind) ) then
828  call mpp_error(fatal, "time_interp_external: mask_in should be all 1 when region_type is not NO_REGION")
829  end if
830  if(interp%interp_method == conserve) then
831  call mpp_error(warning, &
832  "conservative interpolation when region_type is not NO_REGION. Make sure interp matches the region")
833  endif
834  if( field%region_type == outside_region) then
835  do j = js_region, je_region
836  do i = is_region, ie_region
837  mask_in(i,j,:) = 0.0_r8_kind
838  enddo
839  enddo
840  else ! field%region_choice == INSIDE_REGION
841  do j = 1, size(mask_in,2)
842  do i = 1, size(mask_in,1)
843  if( j<js_region .OR. j>je_region .OR. i<is_region .OR. i>ie_region ) mask_in(i,j,:) = 0.0_r8_kind
844  enddo
845  enddo
846  endif
847  endif
848 
849  !! added for mixed mode. Data is always read in as r8 (via ext_fieldtype). if existing horiz_interp_type was
850  !! initialized in r4, needs to cast down in order to match up with saved values in horiz_interp_type.
851  !! creates some temporary arrays since intent(out) vars can't get passed in directly
852  if (interp%horizInterpReals4_type%is_allocated) then
853  ! allocate (there may be a better way to do this, had issues with gnu)
854  allocate(hi_tmp_msk_out(isw:iew,jsw:jew, SIZE(field%src_data,3)))
855  allocate(hi_tmp_data(lbound(field%domain_data,1):ubound(field%domain_data,1), &
856  lbound(field%domain_data,2):ubound(field%domain_data,2), &
857  lbound(field%domain_data,3):ubound(field%domain_data,3), &
858  lbound(field%domain_data,4):ubound(field%domain_data,4)))
859  ! copy over to r4
860  hi_tmp_data = real(field%domain_data, r4_kind)
861  ! do interpolation
862  if(interp%interp_method == conserve) then
863  call horiz_interp(interp, real(field%src_data(:,:,:,ib),r4_kind), hi_tmp_data(isw:iew,jsw:jew,:,ib))
864  hi_tmp_msk_out = 1.0_r4_kind
865  else
866  call horiz_interp(interp, real(field%src_data(:,:,:,ib),r4_kind), hi_tmp_data(isw:iew,jsw:jew,:,ib), &
867  mask_in=real(mask_in,r4_kind), mask_out=hi_tmp_msk_out)
868  end if
869  ! assign any output
870  field%domain_data = real(hi_tmp_data, r8_kind)
871  field%mask(isw:iew,jsw:jew,:,ib) = hi_tmp_msk_out(isw:iew,jsw:jew,:) > 0.0_r4_kind
872 
873  if(allocated(hi_tmp_data)) deallocate(hi_tmp_data)
874  if(allocated(hi_tmp_msk_out)) deallocate(hi_tmp_msk_out)
875  else
876  allocate(mask_out(isw:iew,jsw:jew, size(field%src_data,3)))
877  if(interp%interp_method == conserve) then
878  call horiz_interp(interp, field%src_data(:,:,:,ib), field%domain_data(isw:iew,jsw:jew,:,ib))
879  mask_out = 1.0_r8_kind
880  else
881  call horiz_interp(interp, field%src_data(:,:,:,ib),field%domain_data(isw:iew,jsw:jew,:,ib), &
882  mask_in=mask_in, mask_out=mask_out)
883  end if
884  field%mask(isw:iew,jsw:jew,:,ib) = mask_out(isw:iew,jsw:jew,:) > 0.0_r8_kind
885  deallocate(mask_out)
886  endif
887 
888  else
889  if ( field%region_type .NE. no_region ) then
890  call mpp_error(fatal, "time_interp_external: region_type should be NO_REGION when interp is not present")
891  endif
892  field%domain_data(isw:iew,jsw:jew,:,ib) = field%src_data(isw:iew,jsw:jew,:,ib)
893  field%mask(isw:iew,jsw:jew,:,ib) = is_valid(field%domain_data(isw:iew,jsw:jew,:,ib),field%valid)
894  endif
895  ! convert units
896  where(field%mask(isw:iew,jsw:jew,:,ib)) field%domain_data(isw:iew,jsw:jew,:,ib) = &
897  field%domain_data(isw:iew,jsw:jew,:,ib)*field%slope + field%intercept
898  field%need_compute(ib, window_id) = .false.
899 
900  endif
901 
902 end subroutine load_record
903 
904 !> Given a initialized ext_fieldtype and record number, loads the given index into field%src_data
905 subroutine load_record_0d(field, rec)
906  type(ext_fieldtype), intent(inout) :: field
907  integer , intent(in) :: rec ! record number
908  ! ---- local vars
909  integer :: ib ! index in the array of input buffers
910  integer :: start(4), nread(4)
911 
912  ib = find_buf_index(rec,field%ibuf)
913 
914  if(ib>0) then
915  return
916  else
917  ! calculate current buffer number in round-robin fasion
918  field%nbuf = field%nbuf + 1
919  if(field%nbuf > size(field%domain_data,4).or.field%nbuf <= 0) field%nbuf = 1
920  ib = field%nbuf
921  field%ibuf(ib) = rec
922 
923  if (debug_this_module) write(outunit,*) 'reading record without domain for field ',trim(field%name)
924  start = 1; nread = 1
925  start(3) = 1; nread(3) = size(field%src_data,3)
926  start(field%tdim) = rec; nread(field%tdim) = 1
927  call read_data(field%fileobj,field%name,field%src_data(:,:,:,ib),corner=start,edge_lengths=nread)
928  if ( field%region_type .NE. no_region ) then
929  call mpp_error(fatal, "time_interp_external: region_type should be NO_REGION when field is scalar")
930  endif
931  field%domain_data(1,1,:,ib) = field%src_data(1,1,:,ib)
932  field%mask(1,1,:,ib) = is_valid(field%domain_data(1,1,:,ib),field%valid)
933  ! convert units
934  where(field%mask(1,1,:,ib)) field%domain_data(1,1,:,ib) = &
935  field%domain_data(1,1,:,ib)*field%slope + field%intercept
936  endif
937 
938 end subroutine load_record_0d
939 
940 !> Reallocates src_data for field from module level loaded_fields array
941 subroutine reset_src_data_region(index, is, ie, js, je)
942  integer, intent(in) :: index
943  integer, intent(in) :: is, ie, js, je
944  integer :: nk, nbuf
945 
946  if( is == loaded_fields(index)%is_src .AND. ie == loaded_fields(index)%ie_src .AND. &
947  js == loaded_fields(index)%js_src .AND. ie == loaded_fields(index)%je_src ) return
948 
949  if( .NOT. ASSOCIATED(loaded_fields(index)%src_data) ) call mpp_error(fatal, &
950  "time_interp_external: field(index)%src_data is not associated")
951  nk = size(loaded_fields(index)%src_data,3)
952  nbuf = size(loaded_fields(index)%src_data,4)
953  deallocate(loaded_fields(index)%src_data)
954  allocate(loaded_fields(index)%src_data(is:ie,js:je,nk,nbuf))
955  loaded_fields(index)%is_src = is
956  loaded_fields(index)%ie_src = ie
957  loaded_fields(index)%js_src = js
958  loaded_fields(index)%je_src = je
959 
960 
961 end subroutine reset_src_data_region
962 
963 subroutine set_override_region(index, region_type, is_region, ie_region, js_region, je_region)
964  integer, intent(in) :: index, region_type
965  integer, intent(in) :: is_region, ie_region, js_region, je_region
966 
967  loaded_fields(index)%region_type = region_type
968  loaded_fields(index)%is_region = is_region
969  loaded_fields(index)%ie_region = ie_region
970  loaded_fields(index)%js_region = js_region
971  loaded_fields(index)%je_region = je_region
972 
973  return
974 
975 end subroutine set_override_region
976 
977 !> reallocates array of fields, increasing its size
978 subroutine realloc_files(n)
979  integer, intent(in) :: n ! new size
980 
981  type(filetype), pointer :: ptr(:)
982  integer :: i
983 
984  if (associated(opened_files)) then
985  if (n <= size(opened_files)) return ! do nothing, if requested size no more than current
986  endif
987 
988  allocate(ptr(n))
989  do i = 1, size(ptr)
990  ptr(i)%filename = ''
991  if(Associated(ptr(i)%fileobj)) then
992  call close_file(ptr(i)%fileobj)
993  deallocate(ptr(i)%fileobj)
994  endif
995  enddo
996 
997  if (associated(opened_files))then
998  ptr(1:size(opened_files)) = opened_files(:)
999  deallocate(opened_files)
1000  endif
1001  opened_files => ptr
1002 
1003 end subroutine realloc_files
1004 
1005 !> reallocates array of fields,increasing its size
1006 subroutine realloc_fields(n)
1007  integer, intent(in) :: n ! new size
1008 
1009  type(ext_fieldtype), pointer :: ptr(:)
1010  integer :: i, ier
1011 
1012  if (associated(loaded_fields)) then
1013  if (n <= size(loaded_fields)) return ! do nothing if requested size no more then current
1014  endif
1015 
1016  allocate(ptr(n))
1017  do i=1,size(ptr)
1018  ptr(i)%fileobj => null()
1019  ptr(i)%name=''
1020  ptr(i)%units=''
1021  ptr(i)%siz=-1
1022  ptr(i)%ndim=-1
1023  ptr(i)%domain = null_domain2d
1024  if (ASSOCIATED(ptr(i)%time)) DEALLOCATE(ptr(i)%time, stat=ier)
1025  if (ASSOCIATED(ptr(i)%start_time)) DEALLOCATE(ptr(i)%start_time, stat=ier)
1026  if (ASSOCIATED(ptr(i)%end_time)) DEALLOCATE(ptr(i)%end_time, stat=ier)
1027  if (ASSOCIATED(ptr(i)%period)) DEALLOCATE(ptr(i)%period, stat=ier)
1028  ptr(i)%modulo_time=.false.
1029  if (ASSOCIATED(ptr(i)%domain_data)) DEALLOCATE(ptr(i)%domain_data, stat=ier)
1030  if (ASSOCIATED(ptr(i)%ibuf)) DEALLOCATE(ptr(i)%ibuf, stat=ier)
1031  if (ASSOCIATED(ptr(i)%src_data)) DEALLOCATE(ptr(i)%src_data, stat=ier)
1032  ptr(i)%nbuf=-1
1033  ptr(i)%domain_present=.false.
1034  ptr(i)%slope=1.0_r8_kind
1035  ptr(i)%intercept=0.0_r8_kind
1036  ptr(i)%isc=-1;ptr(i)%iec=-1
1037  ptr(i)%jsc=-1;ptr(i)%jec=-1
1038  enddo
1039  if (associated(loaded_fields)) then
1040  ptr(1:size(loaded_fields)) = loaded_fields(:)
1041  deallocate(loaded_fields)
1042  endif
1043  loaded_fields=>ptr
1044 
1045 end subroutine realloc_fields
1046 
1047 !> simple linear search for given value in given list
1048 !! TODO should use better search if this list is bigger
1049 function find_buf_index(indx,buf)
1050  integer :: indx
1051  integer, dimension(:) :: buf
1052  integer :: find_buf_index
1053 
1054  integer :: nbuf, i
1055 
1056  nbuf = size(buf(:))
1057 
1058  find_buf_index = -1
1059 
1060  do i=1,nbuf
1061  if (buf(i) == indx) then
1062  find_buf_index = i
1063  exit
1064  endif
1065  enddo
1066 
1067 end function find_buf_index
1068 
1069 !> Returns size of field after call to init_external_field.
1070 !! Ordering is X/Y/Z/T.
1071 !! This call only makes sense for non-distributed reads.
1073 
1074  integer :: index !< returned from previous call to init_external_field.
1075  integer :: get_external_field_size(4)
1076 
1077  if (index .lt. 1 .or. index .gt. num_fields) &
1078  call mpp_error(fatal,'invalid index in call to get_external_field_size')
1079 
1080 
1081  get_external_field_size(1) = loaded_fields(index)%siz(1)
1082  get_external_field_size(2) = loaded_fields(index)%siz(2)
1083  get_external_field_size(3) = loaded_fields(index)%siz(3)
1084  get_external_field_size(4) = loaded_fields(index)%siz(4)
1085 
1086 end function get_external_field_size
1087 
1088 !> return missing value
1090 
1091  integer :: index !< returned from previous call to init_external_field.
1092  real(r8_kind) :: get_external_field_missing
1093 
1094  if (index .lt. 1 .or. index .gt. num_fields) &
1095  call mpp_error(fatal,'invalid index in call to get_external_field_size')
1096 
1097 
1098  get_external_field_missing = loaded_fields(index)%missing
1099 
1100 end function get_external_field_missing
1101 
1102 subroutine get_time_axis(index, time)
1103  integer , intent(in) :: index !< field id
1104  type(time_type), intent(out) :: time(:) !< array of time values to be filled
1105 
1106  integer :: n !< size of the data to be assigned
1107 
1108  if (index < 1.or.index > num_fields) &
1109  call mpp_error(fatal,'invalid index in call to get_time_axis')
1110 
1111  n = min(size(time),size(loaded_fields(index)%time))
1112 
1113  time(1:n) = loaded_fields(index)%time(1:n)
1114 end subroutine
1115 
1116 
1117 !> exit time_interp_external_mod. Close all open files and
1118 !! release storage
1120 
1121  integer :: i
1122  !
1123  ! release storage arrays
1124  !
1125  do i=1,num_fields
1126  deallocate(loaded_fields(i)%time,loaded_fields(i)%start_time,loaded_fields(i)%end_time,&
1127  loaded_fields(i)%period,loaded_fields(i)%domain_data,loaded_fields(i)%mask,loaded_fields(i)%ibuf)
1128  if (ASSOCIATED(loaded_fields(i)%src_data)) deallocate(loaded_fields(i)%src_data)
1129  loaded_fields(i)%domain = null_domain2d
1130  loaded_fields(i)%nbuf = 0
1131  loaded_fields(i)%slope = 0.0_r8_kind
1132  loaded_fields(i)%intercept = 0.0_r8_kind
1133  enddo
1134 
1135  !-- close all the files opended
1136  do i = 1, num_files
1137  call close_file(opened_files(i)%fileobj)
1138  deallocate(opened_files(i)%fileobj)
1139  enddo
1140 
1141  deallocate(loaded_fields)
1142  deallocate(opened_files)
1143 
1144  num_fields = 0
1145 
1146  module_initialized = .false.
1147 
1148 end subroutine time_interp_external_exit
1149 
1150 #include "time_interp_external2_r4.fh"
1151 #include "time_interp_external2_r8.fh"
1152 
1153 #include "time_interp_external2_bridge_r4.fh"
1154 #include "time_interp_external2_bridge_r8.fh"
1155 
1156 end module time_interp_external2_mod
1157 !> @}
1158 ! close documentation grouping
subroutine, public get_axis_cart(fileobj, axisname, cart)
Returns X,Y,Z or T cartesian attribute.
logical function, public get_axis_modulo(fileobj, axisname)
Checks if 'modulo' variable exists for a given axis.
logical function, public get_axis_modulo_times(fileobj, axisname, tbeg, tend)
Close a netcdf or domain file opened with open_file or open_virtual_file.
Definition: fms2_io.F90:165
Opens a given netcdf or domain file.
Definition: fms2_io.F90:121
Read data from a defined field in a file.
Definition: fms2_io.F90:291
integer function, public check_nml_error(IOSTAT, NML_NAME)
Checks the iostat argument that is returned after reading a namelist and determines if the error code...
Definition: fms.F90:523
subroutine, public write_version_number(version, tag, unit)
Prints to the log file (or a specified unit) the version id string and tag name.
Definition: fms.F90:701
Added for mixed precision support. Updates force time_manager math to be done with kind=8 reals _wrap...
Subroutine for performing the horizontal interpolation between two grids.
Holds data pointers and metadata for horizontal interpolations, passed between the horiz_interp modul...
These routines retrieve the axis specifications associated with the compute domains....
These routines retrieve the axis specifications associated with the data domains. The domain is a der...
These routines retrieve the axis specifications associated with the global domains....
The domain2D type contains all the necessary information to define the global, compute and data domai...
integer function stdout()
This function returns the current standard fortran unit numbers for output.
Definition: mpp_util.inc:42
integer function stdlog()
This function returns the current standard fortran unit numbers for log messages. Log messages,...
Definition: mpp_util.inc:58
integer function mpp_npes()
Returns processor count for current pelist.
Definition: mpp_util.inc:420
integer function mpp_pe()
Returns processor ID.
Definition: mpp_util.inc:406
Perform parallel broadcasts.
Definition: mpp.F90:1104
Error handler.
Definition: mpp.F90:381
subroutine, public time_interp_external_exit()
exit time_interp_external_mod. Close all open files and release storage
integer function, dimension(4), public get_external_field_size(index)
Returns size of field after call to init_external_field. Ordering is X/Y/Z/T. This call only makes se...
subroutine, public reset_src_data_region(index, is, ie, js, je)
Reallocates src_data for field from module level loaded_fields array.
subroutine load_record(field, rec, interp, is_in, ie_in, js_in, je_in, window_id_in)
load specified record from file Always loads in r8, casts down for horiz_interp if interp argument is...
subroutine, public time_interp_external_init()
Initialize the time_interp_external_mod module.
real(r8_kind) function, public get_external_field_missing(index)
return missing value
subroutine realloc_files(n)
reallocates array of fields, increasing its size
integer function, public init_external_field(file, fieldname, domain, desired_units, verbose, axis_names, axis_sizes, override, correct_leap_year_inconsistency, permit_calendar_conversion, use_comp_domain, ierr, nwindows, ignore_axis_atts, ongrid)
Initialize an external field. Buffer "num_io_buffers" (default=2) in memory to reduce memory allocati...
integer function, private find_buf_index(indx, buf)
simple linear search for given value in given list TODO should use better search if this list is bigg...
subroutine, public get_time_axis(index, time)
subroutine realloc_fields(n)
reallocates array of fields,increasing its size
subroutine load_record_0d(field, rec)
Given a initialized ext_fieldtype and record number, loads the given index into fieldsrc_data.
Provide data from external file interpolated to current model time. Data may be local to current proc...
Holds filename and file object.
Returns a weight and dates or indices for interpolating between two dates. The interface fraction_of_...
subroutine, public get_time(Time, seconds, days, ticks, err_msg)
Returns days and seconds ( < 86400 ) corresponding to a time. err_msg should be checked for any error...
subroutine, public get_date(time, year, month, day, hour, minute, second, tick, err_msg)
Gets the date for different calendar types. Given a time_interval, returns the corresponding date und...
integer function, public days_in_year(Time)
Returns the number of days in the calendar year corresponding to the date represented by time for the...
type(time_type) function, public increment_time(Time, seconds, days, ticks, err_msg, allow_neg_inc)
Increments a time by seconds and days.
integer function, public days_in_month(Time, err_msg)
Given a time, computes the corresponding date given the selected date time mapping algorithm.
integer function, public get_calendar_type()
Returns default calendar type for mapping from time to date.
Given an input date in year, month, days, etc., creates a time_type that represents this time interva...
Given some number of seconds and days, returns the corresponding time_type.
Type to represent amounts of time. Implemented as seconds and days to allow for larger intervals.