| blob | 7c446f4c52e6753a3fc286ebbe57fd7e2b7763fe |
1 !***********************************************************************
2 !* GNU Lesser General Public License
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* FMS is free software: you can redistribute it and/or modify it under
7 !* the terms of the GNU Lesser General Public License as published by
8 !* the Free Software Foundation, either version 3 of the License, or (at
9 !* your option) any later version.
10 !*
11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 !* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 !* for more details.
15 !*
16 !* You should have received a copy of the GNU Lesser General Public
17 !* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18 !***********************************************************************
19 !> @defgroup time_interp_external_mod time_interp_external_mod
20 !> @ingroup time_interp
21 !> @brief Perform I/O and time interpolation of external fields (contained in a file).
22 !> @author M.J. Harrison
23 !!
24 !! Perform I/O and time interpolation for external fields.
25 !! Uses udunits library to calculate calendar dates and
26 !! convert units. Allows for reading data decomposed across
27 !! model horizontal grid using optional domain2d argument
28 !!
29 !! data are defined over data domain for domain2d data
30 !! (halo values are NOT updated by this module)
32 !> @addtogroup time_interp_external_mod
33 !> @{
34 module time_interp_external_mod
35 #ifdef use_deprecated_io
36 #include <fms_platform.h>
37 !
38 !<CONTACT EMAIL="Matthew.Harrison@noaa.gov">M.J. Harrison</CONTACT>
39 !
40 !<REVIEWER EMAIL="hsimmons@iarc.uaf.edu">Harper Simmons</REVIEWER>
41 !
42 !<OVERVIEW>
43 !</OVERVIEW>
45 !<DESCRIPTION>
46 !
47 !</DESCRIPTION>
48 !
49 !<NAMELIST NAME="time_interp_external_nml">
50 ! <DATA NAME="num_io_buffers" TYPE="integer">
51 ! size of record dimension for internal buffer. This is useful for tuning i/o performance
52 ! particularly for large datasets (e.g. daily flux fields)
53 ! </DATA>
54 !</NAMELIST>
56 use fms_mod, only : write_version_number
57 use mpp_mod, only : mpp_error,FATAL,WARNING,mpp_pe, stdout, stdlog, NOTE
58 use mpp_mod, only : input_nml_file
59 use mpp_io_mod, only : mpp_open, mpp_get_atts, mpp_get_info, MPP_NETCDF, MPP_MULTI, MPP_SINGLE,&
60 mpp_get_times, MPP_RDONLY, MPP_ASCII, default_axis,axistype,fieldtype,atttype, &
61 mpp_get_axes, mpp_get_fields, mpp_read, default_field, mpp_close, &
62 mpp_get_tavg_info, validtype, mpp_is_valid, mpp_get_file_name
63 use time_manager_mod, only : time_type, get_date, set_date, operator ( >= ) , operator ( + ) , days_in_month, &
64 operator( - ), operator ( / ) , days_in_year, increment_time, &
65 set_time, get_time, operator( > ), get_calendar_type, NO_CALENDAR
66 use get_cal_time_mod, only : get_cal_time
67 use mpp_domains_mod, only : domain2d, mpp_get_compute_domain, mpp_get_data_domain, &
68 mpp_get_global_domain, NULL_DOMAIN2D
69 use time_interp_mod, only : time_interp, time_interp_init
70 use axis_utils_mod, only : get_axis_cart, get_axis_modulo, get_axis_modulo_times
71 use fms_mod, only : lowercase, open_namelist_file, check_nml_error, close_file
72 use platform_mod, only: r8_kind
73 use horiz_interp_mod, only : horiz_interp, horiz_interp_type
75 implicit none
76 private
78 ! Include variable "version" to be written to log file.
79 #include<file_version.h>
81 integer, parameter, public :: NO_REGION=0, INSIDE_REGION=1, OUTSIDE_REGION=2
82 integer, parameter, private :: modulo_year= 0001
83 integer, parameter, private :: LINEAR_TIME_INTERP = 1 ! not used currently
84 integer, parameter, public :: SUCCESS = 0, ERR_FIELD_NOT_FOUND = 1
85 integer, private :: max_fields = 100, max_files= 40
86 integer, private :: num_fields = 0, num_files=0
87 ! denotes time intervals in file (interpreted from metadata)
88 integer, private :: num_io_buffers = 2 ! set -1 to read all records from disk into memory
89 logical, private :: module_initialized = .false.
90 logical, private :: debug_this_module = .false.
92 public init_external_field, time_interp_external, time_interp_external_init, &
93 time_interp_external_exit, get_external_field_size, get_time_axis, get_external_field_missing
94 public set_override_region, reset_src_data_region, get_external_field_axes
96 private find_buf_index,&
97 set_time_modulo
99 !> @}
101 !> @ingroup time_interp_external_mod
102 type, private :: ext_fieldtype
103 integer :: unit ! keep unit open when not reading all records
104 character(len=128) :: name, units
105 integer :: siz(4), ndim
106 type(domain2d) :: domain
107 type(axistype) :: axes(4)
108 type(time_type), dimension(:), pointer :: time =>NULL() ! midpoint of time interval
109 type(time_type), dimension(:), pointer :: start_time =>NULL(), end_time =>NULL()
110 type(fieldtype) :: field ! mpp_io type
111 type(time_type), dimension(:), pointer :: period =>NULL()
112 logical :: modulo_time ! denote climatological time axis
113 real, dimension(:,:,:,:), pointer :: data =>NULL() ! defined over data domain or global domain
114 logical, dimension(:,:,:,:), pointer :: mask =>NULL() ! defined over data domain or global domain
115 integer, dimension(:), pointer :: ibuf =>NULL() ! record numbers associated with buffers
116 real, dimension(:,:,:,:), pointer :: src_data =>NULL() ! input data buffer
117 type(validtype) :: valid ! data validator
118 integer :: nbuf
119 logical :: domain_present
120 real(DOUBLE_KIND) :: slope, intercept
121 integer :: isc,iec,jsc,jec
122 type(time_type) :: modulo_time_beg, modulo_time_end
123 logical :: have_modulo_times, correct_leap_year_inconsistency
124 integer :: region_type
125 integer :: is_region, ie_region, js_region, je_region
126 integer :: is_src, ie_src, js_src, je_src
127 integer :: tdim
128 integer :: numwindows
129 logical, dimension(:,:), pointer :: need_compute=>NULL()
130 real :: missing ! missing value
131 end type ext_fieldtype
133 !> @ingroup time_interp_external_mod
134 type, private :: filetype
135 character(len=128) :: filename = ''
136 integer :: unit = -1
137 end type filetype
139 !> Provide data from external file interpolated to current model time.
140 !! Data may be local to current processor or global, depending on
141 !! "init_external_field" flags. Uses @ref mpp_io_mod for I/O.
142 !!
143 !! @param index index of external field from previous call to init_external_field
144 !! @param time target time for data
145 !! @param [inout] data global or local data array
146 !! @param interp time_interp_external defined interpolation method (optional). Currently
147 !! this module only supports LINEAR_TIME_INTERP.
148 !! @param verbose verbose flag for debugging (optional).
149 !!
150 !> @ingroup time_interp_external_mod
151 interface time_interp_external
152 module procedure time_interp_external_0d
153 module procedure time_interp_external_2d
154 module procedure time_interp_external_3d
155 end interface
157 !> @addtogroup time_interp_external_mod
158 !> @{
160 integer :: outunit
162 type(ext_fieldtype), save, private, pointer :: field(:) => NULL()
163 type(filetype), save, private, pointer :: opened_files(:) => NULL()
164 !Balaji: really should use field%missing
165 integer, private, parameter :: dk = DOUBLE_KIND
166 real(DOUBLE_KIND), private, parameter :: time_interp_missing=-1e99_dk
167 contains
169 ! <SUBROUTINE NAME="time_interp_external_init">
170 !
171 ! <DESCRIPTION>
172 ! Initialize the time_interp_external module
173 ! </DESCRIPTION>
174 !
175 subroutine time_interp_external_init()
177 integer :: ioun, io_status, logunit, ierr
179 namelist /time_interp_external_nml/ num_io_buffers, debug_this_module, &
180 max_fields, max_files
182 ! open and read namelist
184 if(module_initialized) return
186 logunit = stdlog()
187 outunit = stdout()
188 call write_version_number("TIME_INTERP_EXTERNAL_MOD", version)
190 #ifdef INTERNAL_FILE_NML
191 read (input_nml_file, time_interp_external_nml, iostat=io_status)
192 ierr = check_nml_error(io_status, 'time_interp_external_nml')
193 #else
194 ioun = open_namelist_file ()
195 ierr=1; do while (ierr /= 0)
196 read (ioun, nml=time_interp_external_nml, iostat=io_status, end=10)
197 ierr = check_nml_error(io_status, 'time_interp_external_nml')
198 enddo
199 10 call close_file (ioun)
200 #endif
202 write(logunit,time_interp_external_nml)
203 call realloc_fields(max_fields)
204 call realloc_files(max_files)
206 module_initialized = .true.
208 call time_interp_init()
210 return
212 end subroutine time_interp_external_init
213 ! </SUBROUTINE> NAME="time_interp_external_init"
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>
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,format,threading,domain,desired_units,&
275 verbose,axis_centers,axis_sizes,override,correct_leap_year_inconsistency,&
276 permit_calendar_conversion,use_comp_domain,ierr, nwindows, ignore_axis_atts )
278 character(len=*), intent(in) :: file,fieldname
279 integer, intent(in), optional :: format, threading
280 logical, intent(in), optional :: verbose
281 character(len=*), intent(in), optional :: desired_units
282 type(domain2d), intent(in), optional :: domain
283 type(axistype), intent(inout), optional :: axis_centers(4)
284 integer, intent(inout), optional :: axis_sizes(4)
285 logical, intent(in), optional :: override, correct_leap_year_inconsistency,&
286 permit_calendar_conversion,use_comp_domain
287 integer, intent(out), optional :: ierr
288 integer, intent(in), optional :: nwindows
289 logical, optional :: ignore_axis_atts
290 real :: missing
292 integer :: init_external_field
294 type(fieldtype), dimension(:), allocatable :: flds
295 type(axistype), dimension(:), allocatable :: axes, fld_axes
296 type(axistype) :: time_axis
297 type(atttype), allocatable, dimension(:) :: global_atts
299 real(DOUBLE_KIND) :: slope, intercept
300 integer :: form, thread, fset, unit,ndim,nvar,natt,ntime,i,j
301 integer :: iscomp,iecomp,jscomp,jecomp,isglobal,ieglobal,jsglobal,jeglobal
302 integer :: isdata,iedata,jsdata,jedata, dxsize, dysize,dxsize_max,dysize_max
303 logical :: verb, transpose_xy,use_comp_domain1
304 real, dimension(:), allocatable :: tstamp, tstart, tend, tavg
305 character(len=1) :: cart
306 character(len=1), dimension(4) :: cart_dir
307 character(len=128) :: units, fld_units
308 character(len=128) :: name, msg, calendar_type, timebeg, timeend
309 integer :: siz(4), siz_in(4), gxsize, gysize,gxsize_max, gysize_max
310 type(time_type) :: tdiff
311 integer :: yr, mon, day, hr, minu, sec
312 integer :: len, nfile, nfields_orig, nbuf, nx,ny
313 integer :: numwindows
314 logical :: ignore_axatts
317 if (.not. module_initialized) call mpp_error(FATAL,'Must call time_interp_external_init first')
318 if(present(ierr)) ierr = SUCCESS
319 ignore_axatts=.false.
320 cart_dir(1)='X';cart_dir(2)='Y';cart_dir(3)='Z';cart_dir(4)='T'
321 if(present(ignore_axis_atts)) ignore_axatts = ignore_axis_atts
322 use_comp_domain1 = .false.
323 if(PRESENT(use_comp_domain)) use_comp_domain1 = use_comp_domain
324 form=MPP_NETCDF
325 if (PRESENT(format)) form = format
326 thread = MPP_MULTI
327 if (PRESENT(threading)) thread = threading
328 fset = MPP_SINGLE
329 verb=.false.
330 if (PRESENT(verbose)) verb=verbose
331 if (debug_this_module) verb = .true.
332 numwindows = 1
333 if(present(nwindows)) numwindows = nwindows
335 units = 'same'
336 if (PRESENT(desired_units)) then
337 units = desired_units
338 call mpp_error(FATAL,'==> Unit conversion via time_interp_external &
339 &has been temporarily deprecated. Previous versions of&
340 &this module used udunits_mod to perform unit conversion.&
341 & Udunits_mod is in the process of being replaced since &
342 &there were portability issues associated with this code.&
343 & Please remove the desired_units argument from calls to &
344 &this routine.')
345 endif
346 nfile = 0
347 do i=1,num_files
348 if(trim(opened_files(i)%filename) == trim(file)) then
349 nfile = i
350 exit ! file is already opened
351 endif
352 enddo
353 if(nfile == 0) then
354 call mpp_open(unit,trim(file),MPP_RDONLY,form,threading=thread,&
355 fileset=fset)
356 num_files = num_files + 1
357 if(num_files > max_files) then ! not enough space in the file table, reallocate it
358 !--- z1l: For the case of multiple thread, realoc_files will cause memory leak.
359 !--- If multiple threads are working on file A. One of the thread finished first and
360 !--- begin to work on file B, the realloc_files will cause problem for
361 !--- other threads are working on the file A.
362 ! call realloc_files(2*size(opened_files))
363 call mpp_error(FATAL, "time_interp_external: num_files is greater than max_files, "// &
364 "increase time_interp_external_nml max_files")
365 endif
366 opened_files(num_files)%filename = trim(file)
367 opened_files(num_files)%unit = unit
368 else
369 unit = opened_files(nfile)%unit
370 endif
372 call mpp_get_info(unit,ndim,nvar,natt,ntime)
374 if (ntime < 1) then
375 write(msg,'(a15,a,a58)') 'external field ',trim(fieldname),&
376 ' does not have an associated record dimension (REQUIRED) '
377 call mpp_error(FATAL,trim(msg))
378 endif
379 allocate(global_atts(natt))
380 call mpp_get_atts(unit, global_atts)
381 allocate(axes(ndim))
382 call mpp_get_axes(unit, axes, time_axis)
383 allocate(flds(nvar))
384 call mpp_get_fields(unit,flds)
385 allocate(tstamp(ntime),tstart(ntime),tend(ntime),tavg(ntime))
386 call mpp_get_times(unit,tstamp)
387 transpose_xy = .false.
388 isdata=1; iedata=1; jsdata=1; jedata=1
389 gxsize=1; gysize=1
390 siz_in = 1
392 if (PRESENT(domain)) then
393 call mpp_get_compute_domain(domain,iscomp,iecomp,jscomp,jecomp)
394 nx = iecomp-iscomp+1; ny = jecomp-jscomp+1
395 call mpp_get_data_domain(domain,isdata,iedata,jsdata,jedata,dxsize,dxsize_max,dysize,dysize_max)
396 call mpp_get_global_domain(domain,isglobal,ieglobal,jsglobal,jeglobal,gxsize,gxsize_max,gysize,gysize_max)
397 elseif(use_comp_domain1) then
398 call mpp_error(FATAL,"init_external_field:"//&
399 " use_comp_domain=true but domain is not present")
400 endif
402 init_external_field = -1
403 nfields_orig = num_fields
405 do i=1,nvar
406 call mpp_get_atts(flds(i),name=name,units=fld_units,ndim=ndim,siz=siz_in)
407 call mpp_get_tavg_info(unit,flds(i),flds,tstamp,tstart,tend,tavg)
408 call mpp_get_atts(flds(i),missing=missing)
409 ! why does it convert case of the field name?
410 if (trim(lowercase(name)) /= trim(lowercase(fieldname))) cycle
412 if (verb) write(outunit,*) 'found field ',trim(fieldname), ' in file !!'
413 num_fields = num_fields + 1
414 if(num_fields > max_fields) then
415 !--- z1l: For the case of multiple thread, realoc_fields will cause memory leak.
416 !--- If multiple threads are working on field A. One of the thread finished first and
417 !--- begin to work on field B, the realloc_files will cause problem for
418 !--- other threads are working on the field A.
419 !call realloc_fields(size(field)*2)
420 call mpp_error(FATAL, "time_interp_external: num_fields is greater than max_fields, "// &
421 "increase time_interp_external_nml max_fields")
422 endif
424 init_external_field = num_fields
425 field(num_fields)%unit = unit
426 field(num_fields)%name = trim(name)
427 field(num_fields)%units = trim(fld_units)
428 field(num_fields)%field = flds(i)
429 field(num_fields)%isc = 1
430 field(num_fields)%iec = 1
431 field(num_fields)%jsc = 1
432 field(num_fields)%jec = 1
433 field(num_fields)%region_type = NO_REGION
434 field(num_fields)%is_region = 0
435 field(num_fields)%ie_region = -1
436 field(num_fields)%js_region = 0
437 field(num_fields)%je_region = -1
438 if (PRESENT(domain)) then
439 field(num_fields)%domain_present = .true.
440 field(num_fields)%domain = domain
441 field(num_fields)%isc=iscomp;field(num_fields)%iec = iecomp
442 field(num_fields)%jsc=jscomp;field(num_fields)%jec = jecomp
443 else
444 field(num_fields)%domain_present = .false.
445 endif
447 call mpp_get_atts(flds(i),valid=field(num_fields)%valid )
448 allocate(fld_axes(ndim))
449 call mpp_get_atts(flds(i),axes=fld_axes)
450 if (ndim > 4) call mpp_error(FATAL, &
451 'invalid array rank <=4d fields supported')
452 field(num_fields)%siz = 1
453 field(num_fields)%ndim = ndim
454 field(num_fields)%tdim = 4
455 field(num_fields)%missing = missing
456 do j=1,field(num_fields)%ndim
457 cart = 'N'
458 call get_axis_cart(fld_axes(j), cart)
459 call mpp_get_atts(fld_axes(j),len=len)
460 if (cart == 'N' .and. .not. ignore_axatts) then
461 write(msg,'(a,"/",a)') trim(file),trim(fieldname)
462 call mpp_error(FATAL,'file/field '//trim(msg)// &
463 ' couldnt recognize axis atts in time_interp_external')
464 else if (cart == 'N' .and. ignore_axatts) then
465 cart = cart_dir(j)
466 endif
467 select case (cart)
468 case ('X')
469 if (j.eq.2) transpose_xy = .true.
470 if (.not.PRESENT(domain) .and. .not.PRESENT(override)) then
471 isdata=1;iedata=len
472 iscomp=1;iecomp=len
473 gxsize = len
474 dxsize = len
475 field(num_fields)%isc=iscomp;field(num_fields)%iec=iecomp
476 elseif (PRESENT(override)) then
477 gxsize = len
478 if (PRESENT(axis_sizes)) axis_sizes(1) = len
479 endif
480 field(num_fields)%axes(1) = fld_axes(j)
481 if(use_comp_domain1) then
482 field(num_fields)%siz(1) = nx
483 else
484 field(num_fields)%siz(1) = dxsize
485 endif
486 if (len /= gxsize) then
487 write(msg,'(a,"/",a)') trim(file),trim(fieldname)
488 call mpp_error(FATAL,'time_interp_ext, file/field '//trim(msg)//' x dim doesnt match model')
489 endif
490 case ('Y')
491 field(num_fields)%axes(2) = fld_axes(j)
492 if (.not.PRESENT(domain) .and. .not.PRESENT(override)) then
493 jsdata=1;jedata=len
494 jscomp=1;jecomp=len
495 gysize = len
496 dysize = len
497 field(num_fields)%jsc=jscomp;field(num_fields)%jec=jecomp
498 elseif (PRESENT(override)) then
499 gysize = len
500 if (PRESENT(axis_sizes)) axis_sizes(2) = len
501 endif
502 if(use_comp_domain1) then
503 field(num_fields)%siz(2) = ny
504 else
505 field(num_fields)%siz(2) = dysize
506 endif
507 if (len /= gysize) then
508 write(msg,'(a,"/",a)') trim(file),trim(fieldname)
509 call mpp_error(FATAL,'time_interp_ext, file/field '//trim(msg)//' y dim doesnt match model')
510 endif
511 case ('Z')
512 field(num_fields)%axes(3) = fld_axes(j)
513 field(num_fields)%siz(3) = siz_in(3)
514 case ('T')
515 field(num_fields)%axes(4) = fld_axes(j)
516 field(num_fields)%siz(4) = ntime
517 field(num_fields)%tdim = j
518 end select
519 enddo
520 siz = field(num_fields)%siz
522 if (PRESENT(axis_centers)) then
523 axis_centers = field(num_fields)%axes
524 endif
526 if (PRESENT(axis_sizes) .and. .not.PRESENT(override)) then
527 axis_sizes = field(num_fields)%siz
528 endif
530 deallocate(fld_axes)
531 if (verb) write(outunit,'(a,4i6)') 'field x,y,z,t local size= ',siz
532 if (verb) write(outunit,*) 'field contains data in units = ',trim(field(num_fields)%units)
533 if (transpose_xy) call mpp_error(FATAL,'axis ordering not supported')
534 if (num_io_buffers .le. 1) call mpp_error(FATAL,'time_interp_ext:num_io_buffers should be at least 2')
535 nbuf = min(num_io_buffers,siz(4))
537 field(num_fields)%numwindows = numwindows
538 allocate(field(num_fields)%need_compute(nbuf, numwindows))
539 field(num_fields)%need_compute = .true.
541 allocate(field(num_fields)%data(isdata:iedata,jsdata:jedata,siz(3),nbuf),&
542 field(num_fields)%mask(isdata:iedata,jsdata:jedata,siz(3),nbuf) )
543 field(num_fields)%mask = .false.
544 field(num_fields)%data = 0.0
545 slope=1.0;intercept=0.0
546 ! if (units /= 'same') call convert_units(trim(field(num_fields)%units),trim(units),slope,intercept)
547 ! if (verb.and.units /= 'same') then
548 ! write(outunit,*) 'attempting to convert data to units = ',trim(units)
549 ! write(outunit,'(a,f8.3,a,f8.3)') 'factor = ',slope,' offset= ',intercept
550 ! endif
551 field(num_fields)%slope = slope
552 field(num_fields)%intercept = intercept
553 allocate(field(num_fields)%ibuf(nbuf))
554 field(num_fields)%ibuf = -1
555 field(num_fields)%nbuf = 0 ! initialize buffer number so that first reading fills data(:,:,:,1)
556 if(PRESENT(override)) then
557 field(num_fields)%is_src = 1
558 field(num_fields)%ie_src = gxsize
559 field(num_fields)%js_src = 1
560 field(num_fields)%je_src = gysize
561 allocate(field(num_fields)%src_data(gxsize,gysize,siz(3),nbuf))
562 else
563 field(num_fields)%is_src = isdata
564 field(num_fields)%ie_src = iedata
565 field(num_fields)%js_src = jsdata
566 field(num_fields)%je_src = jedata
567 allocate(field(num_fields)%src_data(isdata:iedata,jsdata:jedata,siz(3),nbuf))
568 endif
570 allocate(field(num_fields)%time(ntime))
571 allocate(field(num_fields)%period(ntime))
572 allocate(field(num_fields)%start_time(ntime))
573 allocate(field(num_fields)%end_time(ntime))
575 call mpp_get_atts(time_axis,units=units,calendar=calendar_type)
576 do j=1,ntime
577 field(num_fields)%time(j) = get_cal_time(tstamp(j),trim(units),trim(calendar_type), &
578 & permit_calendar_conversion)
579 field(num_fields)%start_time(j) = get_cal_time(tstart(j),trim(units),trim(calendar_type), &
580 & permit_calendar_conversion)
581 field(num_fields)%end_time(j) = get_cal_time( tend(j),trim(units),trim(calendar_type), &
582 & permit_calendar_conversion)
583 enddo
585 if (field(num_fields)%modulo_time) then
586 call set_time_modulo(field(num_fields)%Time)
587 call set_time_modulo(field(num_fields)%start_time)
588 call set_time_modulo(field(num_fields)%end_time)
589 endif
591 if(present(correct_leap_year_inconsistency)) then
592 field(num_fields)%correct_leap_year_inconsistency = correct_leap_year_inconsistency
593 else
594 field(num_fields)%correct_leap_year_inconsistency = .false.
595 endif
597 if(get_axis_modulo_times(time_axis, timebeg, timeend)) then
598 if(get_calendar_type() == NO_CALENDAR) then
599 field(num_fields)%modulo_time_beg = set_time(timebeg)
600 field(num_fields)%modulo_time_end = set_time(timeend)
601 else
602 field(num_fields)%modulo_time_beg = set_date(timebeg)
603 field(num_fields)%modulo_time_end = set_date(timeend)
604 endif
605 field(num_fields)%have_modulo_times = .true.
606 else
607 field(num_fields)%have_modulo_times = .false.
608 endif
609 if(ntime == 1) then
610 call mpp_error(NOTE, 'time_interp_external_mod: file '//trim(file)//' has only one time level')
611 else
612 do j= 1, ntime
613 field(num_fields)%period(j) = field(num_fields)%end_time(j)-field(num_fields)%start_time(j)
614 if (field(num_fields)%period(j) > set_time(0,0)) then
615 call get_time(field(num_fields)%period(j), sec, day)
616 sec = sec/2+mod(day,2)*43200
617 day = day/2
618 field(num_fields)%time(j) = field(num_fields)%start_time(j)+&
619 set_time(sec,day)
620 else
621 if (j > 1 .and. j < ntime) then
622 tdiff = field(num_fields)%time(j+1) - field(num_fields)%time(j-1)
623 call get_time(tdiff, sec, day)
624 sec = sec/2+mod(day,2)*43200
625 day = day/2
626 field(num_fields)%period(j) = set_time(sec,day)
627 sec = sec/2+mod(day,2)*43200
628 day = day/2
629 field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
630 field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
631 elseif ( j == 1) then
632 tdiff = field(num_fields)%time(2) - field(num_fields)%time(1)
633 call get_time(tdiff, sec, day)
634 field(num_fields)%period(j) = set_time(sec,day)
635 sec = sec/2+mod(day,2)*43200
636 day = day/2
637 field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
638 field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
639 else
640 tdiff = field(num_fields)%time(ntime) - field(num_fields)%time(ntime-1)
641 call get_time(tdiff, sec, day)
642 field(num_fields)%period(j) = set_time(sec,day)
643 sec = sec/2+mod(day,2)*43200
644 day = day/2
645 field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
646 field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
647 endif
648 endif
649 enddo
650 endif
652 do j=1,ntime-1
653 if (field(num_fields)%time(j) >= field(num_fields)%time(j+1)) then
654 write(msg,'(A,i20)') "times not monotonically increasing. Filename: " &
655 //TRIM(file)//" field: "//TRIM(fieldname)//" timeslice: ", j
656 call mpp_error(FATAL, TRIM(msg))
657 endif
658 enddo
660 field(num_fields)%modulo_time = get_axis_modulo(time_axis)
662 if (verb) then
663 if (field(num_fields)%modulo_time) write(outunit,*) 'data are being treated as modulo in time'
664 do j= 1, ntime
665 write(outunit,*) 'time index, ', j
666 call get_date(field(num_fields)%start_time(j),yr,mon,day,hr,minu,sec)
667 write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
668 'start time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
669 call get_date(field(num_fields)%time(j),yr,mon,day,hr,minu,sec)
670 write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
671 'mid time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
672 call get_date(field(num_fields)%end_time(j),yr,mon,day,hr,minu,sec)
673 write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
674 'end time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
675 enddo
676 end if
678 enddo
680 if (num_fields == nfields_orig) then
681 if (present(ierr)) then
682 ierr = ERR_FIELD_NOT_FOUND
683 else
684 call mpp_error(FATAL,'external field "'//trim(fieldname)//'" not found in file "'//trim(file)//'"')
685 endif
686 endif
688 deallocate(global_atts)
689 deallocate(axes)
690 deallocate(flds)
691 deallocate(tstamp, tstart, tend, tavg)
693 return
695 end function init_external_field
697 !</FUNCTION> NAME="init_external_field"
700 !> @brief 2D time interpolation for @ref time_interp_external
701 subroutine time_interp_external_2d(index, time, data_in, interp, verbose,horz_interp, mask_out, &
702 is_in, ie_in, js_in, je_in, window_id)
704 integer, intent(in) :: index
705 type(time_type), intent(in) :: time
706 real, dimension(:,:), intent(inout) :: data_in
707 integer, intent(in), optional :: interp
708 logical, intent(in), optional :: verbose
709 type(horiz_interp_type),intent(in), optional :: horz_interp
710 logical, dimension(:,:), intent(out), optional :: mask_out ! set to true where output data is valid
711 integer, intent(in), optional :: is_in, ie_in, js_in, je_in
712 integer, intent(in), optional :: window_id
714 real , dimension(size(data_in,1), size(data_in,2), 1) :: data_out
715 logical, dimension(size(data_in,1), size(data_in,2), 1) :: mask3d
717 data_out(:,:,1) = data_in(:,:) ! fill initial values for the portions of array that are not touched by 3d routine
718 call time_interp_external_3d(index, time, data_out, interp, verbose, horz_interp, mask3d, &
719 is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id)
720 data_in(:,:) = data_out(:,:,1)
721 if (PRESENT(mask_out)) mask_out(:,:) = mask3d(:,:,1)
723 return
724 end subroutine time_interp_external_2d
726 !<SUBROUTINE NAME="time_interp_external" >
727 !
728 !<DESCRIPTION>
729 ! Provide data from external file interpolated to current model time.
730 ! Data may be local to current processor or global, depending on
731 ! "init_external_field" flags.
732 !</DESCRIPTION>
733 !
734 !<IN NAME="index" TYPE="integer">
735 ! index of external field from previous call to init_external_field
736 !</IN>
737 !<IN NAME="time" TYPE="time_manager_mod:time_type">
738 ! target time for data
739 !</IN>
740 !<INOUT NAME="data" TYPE="real" DIM="(:,:),(:,:,:)">
741 ! global or local data array
742 !</INOUT>
743 !<IN NAME="interp" TYPE="integer">
744 ! time_interp_external defined interpolation method (optional). Currently this module only supports
745 ! LINEAR_TIME_INTERP.
746 !</IN>
747 !<IN NAME="verbose" TYPE="logical">
748 ! verbose flag for debugging (optional).
749 !</IN>
751 !> @brief 3D time interpolation for @ref time_interp_external
752 subroutine time_interp_external_3d(index, time, data, interp,verbose,horz_interp, mask_out, is_in, ie_in, &
753 & js_in, je_in, window_id)
755 integer, intent(in) :: index
756 type(time_type), intent(in) :: time
757 real, dimension(:,:,:), intent(inout) :: data
758 integer, intent(in), optional :: interp
759 logical, intent(in), optional :: verbose
760 type(horiz_interp_type), intent(in), optional :: horz_interp
761 logical, dimension(:,:,:), intent(out), optional :: mask_out ! set to true where output data is valid
762 integer, intent(in), optional :: is_in, ie_in, js_in, je_in
763 integer, intent(in), optional :: window_id
765 integer :: nx, ny, nz, interp_method, t1, t2
766 integer :: i1, i2, isc, iec, jsc, jec, mod_time
767 integer :: yy, mm, dd, hh, min, ss
768 character(len=256) :: err_msg, filename
770 integer :: isw, iew, jsw, jew, nxw, nyw
771 ! these are boundaries of the updated portion of the "data" argument
772 ! they are calculated using sizes of the "data" and isc,iec,jsc,jsc
773 ! fileds from respective input field, to center the updated portion
774 ! in the output array
776 real :: w1,w2
777 logical :: verb
778 character(len=16) :: message1, message2
780 nx = size(data,1)
781 ny = size(data,2)
782 nz = size(data,3)
784 interp_method = LINEAR_TIME_INTERP
785 if (PRESENT(interp)) interp_method = interp
786 verb=.false.
787 if (PRESENT(verbose)) verb=verbose
788 if (debug_this_module) verb = .true.
790 if (index < 1.or.index > num_fields) &
791 call mpp_error(FATAL, &
792 & 'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize')
794 isc=field(index)%isc;iec=field(index)%iec
795 jsc=field(index)%jsc;jec=field(index)%jec
797 if( field(index)%numwindows == 1 ) then
798 nxw = iec-isc+1
799 nyw = jec-jsc+1
800 else
801 if( .not. present(is_in) .or. .not. present(ie_in) .or. .not. present(js_in) .or. .not. present(je_in) ) then
802 call mpp_error(FATAL, 'time_interp_external: is_in, ie_in, js_in and je_in must be present ' // &
803 'when numwindows > 1, field='//trim(field(index)%name))
804 endif
805 nxw = ie_in - is_in + 1
806 nyw = je_in - js_in + 1
807 isc = isc + is_in - 1
808 iec = isc + ie_in - is_in
809 jsc = jsc + js_in - 1
810 jec = jsc + je_in - js_in
811 endif
813 isw = (nx-nxw)/2+1; iew = isw+nxw-1
814 jsw = (ny-nyw)/2+1; jew = jsw+nyw-1
816 if (nx < nxw .or. ny < nyw .or. nz < field(index)%siz(3)) then
817 write(message1,'(i6,2i5)') nx,ny,nz
818 call mpp_error(FATAL,'field '//trim(field(index)%name)//' Array size mismatch in time_interp_external.'// &
819 ' Array "data" is too small. shape(data)='//message1)
820 endif
821 if(PRESENT(mask_out)) then
822 if (size(mask_out,1) /= nx .or. size(mask_out,2) /= ny .or. size(mask_out,3) /= nz) then
823 write(message1,'(i6,2i5)') nx,ny,nz
824 write(message2,'(i6,2i5)') size(mask_out,1),size(mask_out,2),size(mask_out,3)
825 call mpp_error(FATAL,'field '//trim(field(index)%name)//' array size mismatch in time_interp_external.'// &
826 ' Shape of array "mask_out" does not match that of array "data".'// &
827 ' shape(data)='//message1//' shape(mask_out)='//message2)
828 endif
829 endif
831 if (field(index)%siz(4) == 1) then
832 ! only one record in the file => time-independent field
833 call load_record(field(index),1,horz_interp, is_in, ie_in ,js_in, je_in,window_id)
834 i1 = find_buf_index(1,field(index)%ibuf)
835 if( field(index)%region_type == NO_REGION ) then
836 where(field(index)%mask(isc:iec,jsc:jec,:,i1))
837 data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)
838 elsewhere
839 ! data(isw:iew,jsw:jew,:) = time_interp_missing !field(index)%missing? Balaji
840 data(isw:iew,jsw:jew,:) = field(index)%missing
841 end where
842 else
843 where(field(index)%mask(isc:iec,jsc:jec,:,i1))
844 data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)
845 end where
846 endif
847 if(PRESENT(mask_out)) &
848 mask_out(isw:iew,jsw:jew,:) = field(index)%mask(isc:iec,jsc:jec,:,i1)
849 else
850 if(field(index)%have_modulo_times) then
851 call time_interp(time,field(index)%modulo_time_beg, field(index)%modulo_time_end, field(index)%time(:), &
852 w2, t1, t2, field(index)%correct_leap_year_inconsistency, err_msg=err_msg)
853 if(err_msg .NE. '') then
854 filename = mpp_get_file_name(field(index)%unit)
855 call mpp_error(FATAL,"time_interp_external 1: "//trim(err_msg)//&
856 ",file="//trim(filename)//",field="//trim(field(index)%name) )
857 endif
858 else
859 if(field(index)%modulo_time) then
860 mod_time=1
861 else
862 mod_time=0
863 endif
864 call time_interp(time,field(index)%time(:),w2,t1,t2,modtime=mod_time, err_msg=err_msg)
865 if(err_msg .NE. '') then
866 filename = mpp_get_file_name(field(index)%unit)
867 call mpp_error(FATAL,"time_interp_external 2: "//trim(err_msg)//&
868 ",file="//trim(filename)//",field="//trim(field(index)%name) )
869 endif
870 endif
871 w1 = 1.0-w2
872 if (verb) then
873 call get_date(time,yy,mm,dd,hh,min,ss)
874 write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
875 'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss
876 write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
877 endif
879 call load_record(field(index),t1,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
880 call load_record(field(index),t2,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
881 i1 = find_buf_index(t1,field(index)%ibuf)
882 i2 = find_buf_index(t2,field(index)%ibuf)
883 if(i1<0.or.i2<0) &
884 call mpp_error(FATAL,'time_interp_external : records were not loaded correctly in memory')
886 if (verb) then
887 write(outunit,*) 'ibuf= ',field(index)%ibuf
888 write(outunit,*) 'i1,i2= ',i1, i2
889 endif
891 if( field(index)%region_type == NO_REGION ) then
892 where(field(index)%mask(isc:iec,jsc:jec,:,i1).and.field(index)%mask(isc:iec,jsc:jec,:,i2))
893 data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)*w1 + &
894 field(index)%data(isc:iec,jsc:jec,:,i2)*w2
895 elsewhere
896 ! data(isw:iew,jsw:jew,:) = time_interp_missing !field(index)%missing? Balaji
897 data(isw:iew,jsw:jew,:) = field(index)%missing
898 end where
899 else
900 where(field(index)%mask(isc:iec,jsc:jec,:,i1).and.field(index)%mask(isc:iec,jsc:jec,:,i2))
901 data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)*w1 + &
902 field(index)%data(isc:iec,jsc:jec,:,i2)*w2
903 end where
904 endif
905 if(PRESENT(mask_out)) &
906 mask_out(isw:iew,jsw:jew,:) = &
907 field(index)%mask(isc:iec,jsc:jec,:,i1).and.&
908 field(index)%mask(isc:iec,jsc:jec,:,i2)
909 endif
911 end subroutine time_interp_external_3d
912 !</SUBROUTINE> NAME="time_interp_external"
914 !> @brief Scalar time interpolation for @ref time_interp_external
915 subroutine time_interp_external_0d(index, time, data, verbose)
917 integer, intent(in) :: index
918 type(time_type), intent(in) :: time
919 real, intent(inout) :: data
920 logical, intent(in), optional :: verbose
922 integer :: t1, t2
923 integer :: i1, i2, mod_time
924 integer :: yy, mm, dd, hh, min, ss
925 character(len=256) :: err_msg, filename
927 real :: w1,w2
928 logical :: verb
930 verb=.false.
931 if (PRESENT(verbose)) verb=verbose
932 if (debug_this_module) verb = .true.
934 if (index < 1.or.index > num_fields) &
935 call mpp_error(FATAL, &
936 & 'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize')
938 if (field(index)%siz(4) == 1) then
939 ! only one record in the file => time-independent field
940 call load_record_0d(field(index),1)
941 i1 = find_buf_index(1,field(index)%ibuf)
942 data = field(index)%data(1,1,1,i1)
943 else
944 if(field(index)%have_modulo_times) then
945 call time_interp(time,field(index)%modulo_time_beg, field(index)%modulo_time_end, field(index)%time(:), &
946 w2, t1, t2, field(index)%correct_leap_year_inconsistency, err_msg=err_msg)
947 if(err_msg .NE. '') then
948 filename = mpp_get_file_name(field(index)%unit)
949 call mpp_error(FATAL,"time_interp_external 3:"//trim(err_msg)//&
950 ",file="//trim(filename)//",field="//trim(field(index)%name) )
951 endif
952 else
953 if(field(index)%modulo_time) then
954 mod_time=1
955 else
956 mod_time=0
957 endif
958 call time_interp(time,field(index)%time(:),w2,t1,t2,modtime=mod_time, err_msg=err_msg)
959 if(err_msg .NE. '') then
960 filename = mpp_get_file_name(field(index)%unit)
961 call mpp_error(FATAL,"time_interp_external 4:"//trim(err_msg)// &
962 ",file="//trim(filename)//",field="//trim(field(index)%name) )
963 endif
964 endif
965 w1 = 1.0-w2
966 if (verb) then
967 call get_date(time,yy,mm,dd,hh,min,ss)
968 write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
969 'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss
970 write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
971 endif
972 call load_record_0d(field(index),t1)
973 call load_record_0d(field(index),t2)
974 i1 = find_buf_index(t1,field(index)%ibuf)
975 i2 = find_buf_index(t2,field(index)%ibuf)
977 if(i1<0.or.i2<0) &
978 call mpp_error(FATAL,'time_interp_external : records were not loaded correctly in memory')
979 data = field(index)%data(1,1,1,i1)*w1 + field(index)%data(1,1,1,i2)*w2
980 if (verb) then
981 write(outunit,*) 'ibuf= ',field(index)%ibuf
982 write(outunit,*) 'i1,i2= ',i1, i2
983 endif
984 endif
986 end subroutine time_interp_external_0d
988 subroutine set_time_modulo(Time)
990 type(time_type), intent(inout), dimension(:) :: Time
992 integer :: ntime, n
993 integer :: yr, mon, dy, hr, minu, sec
995 ntime = size(Time(:))
997 do n = 1, ntime
998 call get_date(Time(n), yr, mon, dy, hr, minu, sec)
999 yr = modulo_year
1000 Time(n) = set_date(yr, mon, dy, hr, minu, sec)
1001 enddo
1004 end subroutine set_time_modulo
1006 ! ============================================================================
1007 ! load specified record from file
1008 subroutine load_record(field, rec, interp, is_in, ie_in, js_in, je_in, window_id_in)
1009 type(ext_fieldtype), intent(inout) :: field
1010 integer , intent(in) :: rec ! record number
1011 type(horiz_interp_type), intent(in), optional :: interp
1012 integer, intent(in), optional :: is_in, ie_in, js_in, je_in
1013 integer, intent(in), optional :: window_id_in
1015 ! ---- local vars
1016 integer :: ib ! index in the array of input buffers
1017 integer :: isw,iew,jsw,jew ! boundaries of the domain on each window
1018 integer :: is_region, ie_region, js_region, je_region, i, j
1019 integer :: start(4), nread(4)
1020 logical :: need_compute
1021 real :: mask_in(size(field%src_data,1),size(field%src_data,2),size(field%src_data,3))
1022 real, allocatable :: mask_out(:,:,:)
1023 integer :: window_id
1025 window_id = 1
1026 if( PRESENT(window_id_in) ) window_id = window_id_in
1027 need_compute = .true.
1029 !$OMP CRITICAL
1030 ib = find_buf_index(rec,field%ibuf)
1032 if(ib>0) then
1033 !--- do nothing
1034 need_compute = .false.
1035 else
1036 ! calculate current buffer number in round-robin fasion
1037 field%nbuf = field%nbuf + 1
1038 if(field%nbuf > size(field%data,4).or.field%nbuf <= 0) field%nbuf = 1
1039 ib = field%nbuf
1040 field%ibuf(ib) = rec
1041 field%need_compute(ib,:) = .true.
1043 if (field%domain_present .and. .not.PRESENT(interp)) then
1044 if (debug_this_module) write(outunit,*) 'reading record with domain for field ',trim(field%name)
1045 call mpp_read(field%unit,field%field,field%domain,field%src_data(:,:,:,ib),rec)
1046 else
1047 if (debug_this_module) write(outunit,*) 'reading record without domain for field ',trim(field%name)
1048 start = 1; nread = 1
1049 start(1) = field%is_src; nread(1) = field%ie_src - field%is_src + 1
1050 start(2) = field%js_src; nread(2) = field%je_src - field%js_src + 1
1051 start(3) = 1; nread(3) = size(field%src_data,3)
1052 start(field%tdim) = rec; nread(field%tdim) = 1
1053 call mpp_read(field%unit,field%field,field%src_data(:,:,:,ib),start,nread)
1054 endif
1055 endif
1056 !$OMP END CRITICAL
1057 isw=field%isc;iew=field%iec
1058 jsw=field%jsc;jew=field%jec
1060 if( field%numwindows > 1) then
1061 if( .NOT. PRESENT(is_in) .OR. .NOT. PRESENT(ie_in) .OR. .NOT. PRESENT(js_in) .OR. .NOT. PRESENT(je_in) ) then
1062 call mpp_error(FATAL, &
1063 & 'time_interp_external(load_record): is_in, ie_in, js_in, je_in must be present when numwindows>1')
1064 endif
1065 isw = isw + is_in - 1
1066 iew = isw + ie_in - is_in
1067 jsw = jsw + js_in - 1
1068 jew = jsw + je_in - js_in
1069 endif
1071 ! interpolate to target grid
1073 need_compute = field%need_compute(ib, window_id)
1074 if(need_compute) then
1075 if(PRESENT(interp)) then
1076 is_region = field%is_region; ie_region = field%ie_region
1077 js_region = field%js_region; je_region = field%je_region
1078 mask_in = 0.0
1079 where (mpp_is_valid(field%src_data(:,:,:,ib), field%valid)) mask_in = 1.0
1080 if ( field%region_type .NE. NO_REGION ) then
1081 if( ANY(mask_in == 0.0) ) then
1082 call mpp_error(FATAL, "time_interp_external: mask_in should be all 1 when region_type is not NO_REGION")
1083 endif
1084 if( field%region_type == OUTSIDE_REGION) then
1085 do j = js_region, je_region
1086 do i = is_region, ie_region
1087 mask_in(i,j,:) = 0.0
1088 enddo
1089 enddo
1090 else ! field%region_choice == INSIDE_REGION
1091 do j = 1, size(mask_in,2)
1092 do i = 1, size(mask_in,1)
1093 if( j<js_region .OR. j>je_region .OR. i<is_region .OR. i>ie_region ) mask_in(i,j,:) = 0.0
1094 enddo
1095 enddo
1096 endif
1097 endif
1098 allocate(mask_out(isw:iew,jsw:jew, size(field%src_data,3)))
1099 call horiz_interp(interp,field%src_data(:,:,:,ib),field%data(isw:iew,jsw:jew,:,ib), &
1100 mask_in=mask_in, &
1101 mask_out=mask_out)
1103 field%mask(isw:iew,jsw:jew,:,ib) = mask_out(isw:iew,jsw:jew,:) > 0
1104 deallocate(mask_out)
1105 else
1106 if ( field%region_type .NE. NO_REGION ) then
1107 call mpp_error(FATAL, "time_interp_external: region_type should be NO_REGION when interp is not present")
1108 endif
1109 field%data(isw:iew,jsw:jew,:,ib) = field%src_data(isw:iew,jsw:jew,:,ib)
1110 field%mask(isw:iew,jsw:jew,:,ib) = mpp_is_valid(field%data(isw:iew,jsw:jew,:,ib),field%valid)
1111 endif
1112 ! convert units
1113 where(field%mask(isw:iew,jsw:jew,:,ib)) field%data(isw:iew,jsw:jew,:,ib) = &
1114 field%data(isw:iew,jsw:jew,:,ib)*field%slope + field%intercept
1115 field%need_compute(ib, window_id) = .false.
1116 endif
1118 end subroutine load_record
1121 subroutine load_record_0d(field, rec)
1122 type(ext_fieldtype), intent(inout) :: field
1123 integer , intent(in) :: rec ! record number
1124 ! ---- local vars
1125 integer :: ib ! index in the array of input buffers
1126 integer :: start(4), nread(4)
1128 ib = find_buf_index(rec,field%ibuf)
1130 if(ib>0) then
1131 return
1132 else
1133 ! calculate current buffer number in round-robin fasion
1134 field%nbuf = field%nbuf + 1
1135 if(field%nbuf > size(field%data,4).or.field%nbuf <= 0) field%nbuf = 1
1136 ib = field%nbuf
1137 field%ibuf(ib) = rec
1139 if (debug_this_module) write(outunit,*) 'reading record without domain for field ',trim(field%name)
1140 start = 1; nread = 1
1141 start(3) = 1; nread(3) = size(field%src_data,3)
1142 start(field%tdim) = rec; nread(field%tdim) = 1
1143 call mpp_read(field%unit,field%field,field%src_data(:,:,:,ib),start,nread)
1144 if ( field%region_type .NE. NO_REGION ) then
1145 call mpp_error(FATAL, "time_interp_external: region_type should be NO_REGION when field is scalar")
1146 endif
1147 field%data(1,1,:,ib) = field%src_data(1,1,:,ib)
1148 field%mask(1,1,:,ib) = mpp_is_valid(field%data(1,1,:,ib),field%valid)
1149 ! convert units
1150 where(field%mask(1,1,:,ib)) field%data(1,1,:,ib) = &
1151 field%data(1,1,:,ib)*field%slope + field%intercept
1152 endif
1154 end subroutine load_record_0d
1156 ! ============================================================================
1157 subroutine reset_src_data_region(index, is, ie, js, je)
1158 integer, intent(in) :: index
1159 integer, intent(in) :: is, ie, js, je
1160 integer :: nk, nbuf
1162 if( is == field(index)%is_src .AND. ie == field(index)%ie_src .AND. &
1163 js == field(index)%js_src .AND. ie == field(index)%je_src ) return
1165 if( .NOT. ASSOCIATED(field(index)%src_data) ) call mpp_error(FATAL, &
1166 "time_interp_external: field(index)%src_data is not associated")
1167 nk = size(field(index)%src_data,3)
1168 nbuf = size(field(index)%src_data,4)
1169 deallocate(field(index)%src_data)
1170 allocate(field(index)%src_data(is:ie,js:je,nk,nbuf))
1171 field(index)%is_src = is
1172 field(index)%ie_src = ie
1173 field(index)%js_src = js
1174 field(index)%je_src = je
1177 end subroutine reset_src_data_region
1179 ! ============================================================================
1180 subroutine set_override_region(index, region_type, is_region, ie_region, js_region, je_region)
1181 integer, intent(in) :: index, region_type
1182 integer, intent(in) :: is_region, ie_region, js_region, je_region
1184 field(index)%region_type = region_type
1185 field(index)%is_region = is_region
1186 field(index)%ie_region = ie_region
1187 field(index)%js_region = js_region
1188 field(index)%je_region = je_region
1190 return
1192 end subroutine set_override_region
1194 ! ============================================================================
1195 ! reallocates array of fields, increasing its size
1196 subroutine realloc_files(n)
1197 integer, intent(in) :: n ! new size
1199 type(filetype), pointer :: ptr(:)
1200 integer :: i
1202 if (associated(opened_files)) then
1203 if (n <= size(opened_files)) return ! do nothing, if requested size no more than current
1204 endif
1206 allocate(ptr(n))
1207 do i = 1, size(ptr)
1208 ptr(i)%filename = ''
1209 ptr(i)%unit = -1
1210 enddo
1212 if (associated(opened_files))then
1213 ptr(1:size(opened_files)) = opened_files(:)
1214 deallocate(opened_files)
1215 endif
1216 opened_files => ptr
1218 end subroutine realloc_files
1220 ! ============================================================================
1221 ! reallocates array of fields,increasing its size
1222 subroutine realloc_fields(n)
1223 integer, intent(in) :: n ! new size
1225 type(ext_fieldtype), pointer :: ptr(:)
1226 integer :: i, ier
1228 if (associated(field)) then
1229 if (n <= size(field)) return ! do nothing if requested size no more then current
1230 endif
1232 allocate(ptr(n))
1233 do i=1,size(ptr)
1234 ptr(i)%unit=-1
1235 ptr(i)%name=''
1236 ptr(i)%units=''
1237 ptr(i)%siz=-1
1238 ptr(i)%ndim=-1
1239 ptr(i)%domain = NULL_DOMAIN2D
1240 ptr(i)%axes(:) = default_axis
1241 if (ASSOCIATED(ptr(i)%time)) DEALLOCATE(ptr(i)%time, stat=ier)
1242 if (ASSOCIATED(ptr(i)%start_time)) DEALLOCATE(ptr(i)%start_time, stat=ier)
1243 if (ASSOCIATED(ptr(i)%end_time)) DEALLOCATE(ptr(i)%end_time, stat=ier)
1244 ptr(i)%field = default_field
1245 if (ASSOCIATED(ptr(i)%period)) DEALLOCATE(ptr(i)%period, stat=ier)
1246 ptr(i)%modulo_time=.false.
1247 if (ASSOCIATED(ptr(i)%data)) DEALLOCATE(ptr(i)%data, stat=ier)
1248 if (ASSOCIATED(ptr(i)%ibuf)) DEALLOCATE(ptr(i)%ibuf, stat=ier)
1249 if (ASSOCIATED(ptr(i)%src_data)) DEALLOCATE(ptr(i)%src_data, stat=ier)
1250 ptr(i)%nbuf=-1
1251 ptr(i)%domain_present=.false.
1252 ptr(i)%slope=1.0
1253 ptr(i)%intercept=0.0
1254 ptr(i)%isc=-1;ptr(i)%iec=-1
1255 ptr(i)%jsc=-1;ptr(i)%jec=-1
1256 enddo
1257 if (associated(field)) then
1258 ptr(1:size(field)) = field(:)
1259 deallocate(field)
1260 endif
1261 field=>ptr
1263 end subroutine realloc_fields
1266 function find_buf_index(indx,buf)
1267 integer :: indx
1268 integer, dimension(:) :: buf
1269 integer :: find_buf_index
1271 integer :: nbuf, i
1273 nbuf = size(buf(:))
1275 find_buf_index = -1
1277 do i=1,nbuf
1278 if (buf(i) == indx) then
1279 find_buf_index = i
1280 exit
1281 endif
1282 enddo
1284 end function find_buf_index
1286 !<FUNCTION NAME="get_external_field_size" TYPE="integer" DIM="(4)">
1287 !
1288 !<DESCRIPTION>
1289 ! return size of field after call to init_external_field.
1290 ! Ordering is X/Y/Z/T.
1291 ! This call only makes sense for non-distributed reads.
1292 !</DESCRIPTION>
1293 !
1294 !<IN NAME="index" TYPE="integer">
1295 ! returned from previous call to init_external_field.
1296 !</IN>
1298 function get_external_field_size(index)
1300 integer :: index
1301 integer :: get_external_field_size(4)
1303 if (index .lt. 1 .or. index .gt. num_fields) &
1304 call mpp_error(FATAL,'invalid index in call to get_external_field_size')
1307 get_external_field_size(1) = field(index)%siz(1)
1308 get_external_field_size(2) = field(index)%siz(2)
1309 get_external_field_size(3) = field(index)%siz(3)
1310 get_external_field_size(4) = field(index)%siz(4)
1312 end function get_external_field_size
1313 !</FUNCTION> NAME="get_external_field_size"
1316 !<FUNCTION NAME="get_external_field_missing" TYPE="real">
1317 !
1318 !<DESCRIPTION>
1319 ! return missing value
1320 !</DESCRIPTION>
1321 !
1322 !<IN NAME="index" TYPE="integer">
1323 ! returned from previous call to init_external_field.
1324 !</IN>
1326 function get_external_field_missing(index)
1328 integer :: index
1329 real :: get_external_field_missing
1331 if (index .lt. 1 .or. index .gt. num_fields) &
1332 call mpp_error(FATAL,'invalid index in call to get_external_field_size')
1335 ! call mpp_get_atts(field(index)%field,missing=missing)
1336 get_external_field_missing = field(index)%missing
1338 end function get_external_field_missing
1339 !</FUNCTION> NAME="get_external_field_missing"
1341 !<FUNCTION NAME="get_external_field_axes" TYPE="axistype" DIM="(4)">
1342 !
1343 !<DESCRIPTION>
1344 ! return field axes after call to init_external_field.
1345 ! Ordering is X/Y/Z/T.
1346 !</DESCRIPTION>
1347 !
1348 !<IN NAME="index" TYPE="integer">
1349 ! returned from previous call to init_external_field.
1350 !</IN>
1353 function get_external_field_axes(index)
1355 integer :: index
1356 type(axistype), dimension(4) :: get_external_field_axes
1358 if (index .lt. 1 .or. index .gt. num_fields) &
1359 call mpp_error(FATAL,'invalid index in call to get_external_field_size')
1362 get_external_field_axes(1) = field(index)%axes(1)
1363 get_external_field_axes(2) = field(index)%axes(2)
1364 get_external_field_axes(3) = field(index)%axes(3)
1365 get_external_field_axes(4) = field(index)%axes(4)
1367 end function get_external_field_axes
1368 !</FUNCTION> NAME="get_external_field_axes"
1370 ! ===========================================================================
1371 subroutine get_time_axis(index, time)
1372 integer , intent(in) :: index ! field id
1373 type(time_type), intent(out) :: time(:) ! array of time values to be filled
1375 integer :: n ! size of the data to be assigned
1377 if (index < 1.or.index > num_fields) &
1378 call mpp_error(FATAL,'invalid index in call to get_time_axis')
1380 n = min(size(time),size(field(index)%time))
1382 time(1:n) = field(index)%time(1:n)
1383 end subroutine
1385 !<SUBROUTINE NAME="time_interp_external_exit">
1386 !
1387 !<DESCRIPTION>
1388 ! exit time_interp_external_mod. Close all open files and
1389 ! release storage
1390 !</DESCRIPTION>
1392 subroutine time_interp_external_exit()
1394 integer :: i,j
1395 !
1396 ! release storage arrays
1397 !
1398 do i=1,num_fields
1399 deallocate(field(i)%time,field(i)%start_time,field(i)%end_time,&
1400 field(i)%period,field(i)%data,field(i)%mask,field(i)%ibuf)
1401 if (ASSOCIATED(field(i)%src_data)) deallocate(field(i)%src_data)
1402 do j=1,4
1403 field(i)%axes(j) = default_axis
1404 enddo
1405 field(i)%domain = NULL_DOMAIN2D
1406 field(i)%field = default_field
1407 field(i)%nbuf = 0
1408 field(i)%slope = 0.
1409 field(i)%intercept = 0.
1410 enddo
1412 deallocate(field)
1413 deallocate(opened_files)
1415 num_fields = 0
1417 module_initialized = .false.
1419 end subroutine time_interp_external_exit
1420 !</SUBROUTINE> NAME="time_interp_external_exit"
1421 #endif
1422 end module time_interp_external_mod
1423 !> @}
1424 ! close documentation grouping