forked from NOAA-GFDL/FMS
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtime_interp_external.F90
1420 lines (1234 loc) · 57.2 KB
/
time_interp_external.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!***********************************************************************
!* GNU Lesser General Public License
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* FMS is free software: you can redistribute it and/or modify it under
!* the terms of the GNU Lesser General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or (at
!* your option) any later version.
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
!* for more details.
!*
!* You should have received a copy of the GNU Lesser General Public
!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
!***********************************************************************
!> @defgroup time_interp_external_mod time_interp_external_mod
!> @ingroup time_interp
!> @brief Perform I/O and time interpolation of external fields (contained in a file).
!> @author M.J. Harrison
!!
!! Perform I/O and time interpolation for external fields.
!! Uses udunits library to calculate calendar dates and
!! convert units. Allows for reading data decomposed across
!! model horizontal grid using optional domain2d argument
!!
!! data are defined over data domain for domain2d data
!! (halo values are NOT updated by this module)
!> @file
!> @brief File for @ref time_interp_external_mod
!> @addtogroup time_interp_external_mod
!> @{
module time_interp_external_mod
#include <fms_platform.h>
!
!<CONTACT EMAIL="[email protected]">M.J. Harrison</CONTACT>
!
!<REVIEWER EMAIL="[email protected]">Harper Simmons</REVIEWER>
!
!<OVERVIEW>
!</OVERVIEW>
!<DESCRIPTION>
!
!</DESCRIPTION>
!
!<NAMELIST NAME="time_interp_external_nml">
! <DATA NAME="num_io_buffers" TYPE="integer">
! size of record dimension for internal buffer. This is useful for tuning i/o performance
! particularly for large datasets (e.g. daily flux fields)
! </DATA>
!</NAMELIST>
use fms_mod, only : write_version_number
use mpp_mod, only : mpp_error,FATAL,WARNING,mpp_pe, stdout, stdlog, NOTE
use mpp_mod, only : input_nml_file
use mpp_io_mod, only : mpp_open, mpp_get_atts, mpp_get_info, MPP_NETCDF, MPP_MULTI, MPP_SINGLE,&
mpp_get_times, MPP_RDONLY, MPP_ASCII, default_axis,axistype,fieldtype,atttype, &
mpp_get_axes, mpp_get_fields, mpp_read, default_field, mpp_close, &
mpp_get_tavg_info, validtype, mpp_is_valid, mpp_get_file_name
use time_manager_mod, only : time_type, get_date, set_date, operator ( >= ) , operator ( + ) , days_in_month, &
operator( - ), operator ( / ) , days_in_year, increment_time, &
set_time, get_time, operator( > ), get_calendar_type, NO_CALENDAR
use get_cal_time_mod, only : get_cal_time
use mpp_domains_mod, only : domain2d, mpp_get_compute_domain, mpp_get_data_domain, &
mpp_get_global_domain, NULL_DOMAIN2D
use time_interp_mod, only : time_interp, time_interp_init
use axis_utils_mod, only : get_axis_cart, get_axis_modulo, get_axis_modulo_times
use fms_mod, only : lowercase, open_namelist_file, check_nml_error, close_file
use platform_mod, only: r8_kind
use horiz_interp_mod, only : horiz_interp, horiz_interp_type
implicit none
private
! Include variable "version" to be written to log file.
#include<file_version.h>
integer, parameter, public :: NO_REGION=0, INSIDE_REGION=1, OUTSIDE_REGION=2
integer, parameter, private :: modulo_year= 0001
integer, parameter, private :: LINEAR_TIME_INTERP = 1 ! not used currently
integer, parameter, public :: SUCCESS = 0, ERR_FIELD_NOT_FOUND = 1
integer, private :: max_fields = 100, max_files= 40
integer, private :: num_fields = 0, num_files=0
! denotes time intervals in file (interpreted from metadata)
integer, private :: num_io_buffers = 2 ! set -1 to read all records from disk into memory
logical, private :: module_initialized = .false.
logical, private :: debug_this_module = .false.
public init_external_field, time_interp_external, time_interp_external_init, &
time_interp_external_exit, get_external_field_size, get_time_axis, get_external_field_missing
public set_override_region, reset_src_data_region, get_external_field_axes
private find_buf_index,&
set_time_modulo
!> @}
!> @ingroup time_interp_external_mod
type, private :: ext_fieldtype
integer :: unit ! keep unit open when not reading all records
character(len=128) :: name, units
integer :: siz(4), ndim
type(domain2d) :: domain
type(axistype) :: axes(4)
type(time_type), dimension(:), pointer :: time =>NULL() ! midpoint of time interval
type(time_type), dimension(:), pointer :: start_time =>NULL(), end_time =>NULL()
type(fieldtype) :: field ! mpp_io type
type(time_type), dimension(:), pointer :: period =>NULL()
logical :: modulo_time ! denote climatological time axis
real, dimension(:,:,:,:), pointer :: data =>NULL() ! defined over data domain or global domain
logical, dimension(:,:,:,:), pointer :: mask =>NULL() ! defined over data domain or global domain
integer, dimension(:), pointer :: ibuf =>NULL() ! record numbers associated with buffers
real, dimension(:,:,:,:), pointer :: src_data =>NULL() ! input data buffer
type(validtype) :: valid ! data validator
integer :: nbuf
logical :: domain_present
real(DOUBLE_KIND) :: slope, intercept
integer :: isc,iec,jsc,jec
type(time_type) :: modulo_time_beg, modulo_time_end
logical :: have_modulo_times, correct_leap_year_inconsistency
integer :: region_type
integer :: is_region, ie_region, js_region, je_region
integer :: is_src, ie_src, js_src, je_src
integer :: tdim
integer :: numwindows
logical, dimension(:,:), pointer :: need_compute=>NULL()
real :: missing ! missing value
end type ext_fieldtype
!> @ingroup time_interp_external_mod
type, private :: filetype
character(len=128) :: filename = ''
integer :: unit = -1
end type filetype
!> Provide data from external file interpolated to current model time.
!! Data may be local to current processor or global, depending on
!! "init_external_field" flags. Uses @ref mpp_io_mod for I/O.
!!
!! @param index index of external field from previous call to init_external_field
!! @param time target time for data
!! @param [inout] data global or local data array
!! @param interp time_interp_external defined interpolation method (optional). Currently
!! this module only supports LINEAR_TIME_INTERP.
!! @param verbose verbose flag for debugging (optional).
!!
!> @ingroup time_interp_external_mod
interface time_interp_external
module procedure time_interp_external_0d
module procedure time_interp_external_2d
module procedure time_interp_external_3d
end interface
!> @addtogroup time_interp_external_mod
!> @{
integer :: outunit
type(ext_fieldtype), save, private, pointer :: field(:) => NULL()
type(filetype), save, private, pointer :: opened_files(:) => NULL()
!Balaji: really should use field%missing
integer, private, parameter :: dk = DOUBLE_KIND
real(DOUBLE_KIND), private, parameter :: time_interp_missing=-1e99_dk
contains
! <SUBROUTINE NAME="time_interp_external_init">
!
! <DESCRIPTION>
! Initialize the time_interp_external module
! </DESCRIPTION>
!
subroutine time_interp_external_init()
integer :: ioun, io_status, logunit, ierr
namelist /time_interp_external_nml/ num_io_buffers, debug_this_module, &
max_fields, max_files
! open and read namelist
if(module_initialized) return
logunit = stdlog()
outunit = stdout()
call write_version_number("TIME_INTERP_EXTERNAL_MOD", version)
#ifdef INTERNAL_FILE_NML
read (input_nml_file, time_interp_external_nml, iostat=io_status)
ierr = check_nml_error(io_status, 'time_interp_external_nml')
#else
ioun = open_namelist_file ()
ierr=1; do while (ierr /= 0)
read (ioun, nml=time_interp_external_nml, iostat=io_status, end=10)
ierr = check_nml_error(io_status, 'time_interp_external_nml')
enddo
10 call close_file (ioun)
#endif
write(logunit,time_interp_external_nml)
call realloc_fields(max_fields)
call realloc_files(max_files)
module_initialized = .true.
call time_interp_init()
return
end subroutine time_interp_external_init
! </SUBROUTINE> NAME="time_interp_external_init"
!<FUNCTION NAME="init_external_field" TYPE="integer">
!
!<DESCRIPTION>
! initialize an external field. Buffer "num_io_buffers" (default=2) in memory to reduce memory allocations.
! distributed reads are supported using the optional "domain" flag.
! Units conversion via the optional "desired_units" flag using udunits_mod.
!
! Return integer id of field for future calls to time_interp_external.
!
! </DESCRIPTION>
!
!<IN NAME="file" TYPE="character(len=*)">
! filename
!</IN>
!<IN NAME="fieldname" TYPE="character(len=*)">
! fieldname (in file)
!</IN>
!<IN NAME="format" TYPE="integer">
! mpp_io flag for format of file (optional). Currently only "MPP_NETCDF" supported
!</IN>
!<IN NAME="threading" TYPE="integer">
! mpp_io flag for threading (optional). "MPP_SINGLE" means root pe reads global field and distributes to other PEs
! "MPP_MULTI" means all PEs read data
!</IN>
!<IN NAME="domain" TYPE="mpp_domains_mod:domain2d">
! domain flag (optional)
!</IN>
!<IN NAME="desired_units" TYPE="character(len=*)">
! Target units for data (optional), e.g. convert from deg_K to deg_C.
! Failure to convert using udunits will result in failure of this module.
!</IN>
!<IN NAME="verbose" TYPE="logical">
! verbose flag for debugging (optional).
!</IN>
!<INOUT NAME="axis_centers" TYPE="axistype" DIM="(4)">
! MPP_IO axistype array for grid centers ordered X-Y-Z-T (optional).
!</INOUT>
!<INOUT NAME="axis_sizes" TYPE="integer" DIM="(4)">
! array of axis lengths ordered X-Y-Z-T (optional).
!</INOUT>
!> Initialize an external field. Buffer "num_io_buffers" (default=2) in memory to reduce memory allocations.
!! distributed reads are supported using the optional "domain" flag.
!! Units conversion via the optional "desired_units" flag using udunits_mod.
!!
!> @return integer id of field for future calls to time_interp_external.
!> @param file filename
!> @param fieldname fieldname (in file)
!> @param format mpp_io flag for format of file(optional). Currently only "MPP_NETCDF" supported
!> @param threading mpp_io flag for threading (optional). "MPP_SINGLE" means root pe reads
!! global field and distributes to other PEs. "MPP_MULTI" means all PEs read data
!> @param domain domain flag (optional)
!> @param desired_units Target units for data (optional), e.g. convert from deg_K to deg_C.
!! Failure to convert using udunits will result in failure of this module.
!> @param verbose verbose flag for debugging (optional).
!> @param [out] axis_names List of axis names (optional).
!> @param [inout] axis_sizes array of axis lengths ordered X-Y-Z-T (optional).
function init_external_field(file,fieldname,format,threading,domain,desired_units,&
verbose,axis_centers,axis_sizes,override,correct_leap_year_inconsistency,&
permit_calendar_conversion,use_comp_domain,ierr, nwindows, ignore_axis_atts )
character(len=*), intent(in) :: file,fieldname
integer, intent(in), optional :: format, threading
logical, intent(in), optional :: verbose
character(len=*), intent(in), optional :: desired_units
type(domain2d), intent(in), optional :: domain
type(axistype), intent(inout), optional :: axis_centers(4)
integer, intent(inout), optional :: axis_sizes(4)
logical, intent(in), optional :: override, correct_leap_year_inconsistency,&
permit_calendar_conversion,use_comp_domain
integer, intent(out), optional :: ierr
integer, intent(in), optional :: nwindows
logical, optional :: ignore_axis_atts
real :: missing
integer :: init_external_field
type(fieldtype), dimension(:), allocatable :: flds
type(axistype), dimension(:), allocatable :: axes, fld_axes
type(axistype) :: time_axis
type(atttype), allocatable, dimension(:) :: global_atts
real(DOUBLE_KIND) :: slope, intercept
integer :: form, thread, fset, unit,ndim,nvar,natt,ntime,i,j
integer :: iscomp,iecomp,jscomp,jecomp,isglobal,ieglobal,jsglobal,jeglobal
integer :: isdata,iedata,jsdata,jedata, dxsize, dysize,dxsize_max,dysize_max
logical :: verb, transpose_xy,use_comp_domain1
real, dimension(:), allocatable :: tstamp, tstart, tend, tavg
character(len=1) :: cart
character(len=1), dimension(4) :: cart_dir
character(len=128) :: units, fld_units
character(len=128) :: name, msg, calendar_type, timebeg, timeend
integer :: siz(4), siz_in(4), gxsize, gysize,gxsize_max, gysize_max
type(time_type) :: tdiff
integer :: yr, mon, day, hr, minu, sec
integer :: len, nfile, nfields_orig, nbuf, nx,ny
integer :: numwindows
logical :: ignore_axatts
if (.not. module_initialized) call mpp_error(FATAL,'Must call time_interp_external_init first')
if(present(ierr)) ierr = SUCCESS
ignore_axatts=.false.
cart_dir(1)='X';cart_dir(2)='Y';cart_dir(3)='Z';cart_dir(4)='T'
if(present(ignore_axis_atts)) ignore_axatts = ignore_axis_atts
use_comp_domain1 = .false.
if(PRESENT(use_comp_domain)) use_comp_domain1 = use_comp_domain
form=MPP_NETCDF
if (PRESENT(format)) form = format
thread = MPP_MULTI
if (PRESENT(threading)) thread = threading
fset = MPP_SINGLE
verb=.false.
if (PRESENT(verbose)) verb=verbose
if (debug_this_module) verb = .true.
numwindows = 1
if(present(nwindows)) numwindows = nwindows
units = 'same'
if (PRESENT(desired_units)) then
units = desired_units
call mpp_error(FATAL,'==> Unit conversion via time_interp_external &
&has been temporarily deprecated. Previous versions of&
&this module used udunits_mod to perform unit conversion.&
& Udunits_mod is in the process of being replaced since &
&there were portability issues associated with this code.&
& Please remove the desired_units argument from calls to &
&this routine.')
endif
nfile = 0
do i=1,num_files
if(trim(opened_files(i)%filename) == trim(file)) then
nfile = i
exit ! file is already opened
endif
enddo
if(nfile == 0) then
call mpp_open(unit,trim(file),MPP_RDONLY,form,threading=thread,&
fileset=fset)
num_files = num_files + 1
if(num_files > max_files) then ! not enough space in the file table, reallocate it
!--- z1l: For the case of multiple thread, realoc_files will cause memory leak.
!--- If multiple threads are working on file A. One of the thread finished first and
!--- begin to work on file B, the realloc_files will cause problem for
!--- other threads are working on the file A.
! call realloc_files(2*size(opened_files))
call mpp_error(FATAL, "time_interp_external: num_files is greater than max_files, "// &
"increase time_interp_external_nml max_files")
endif
opened_files(num_files)%filename = trim(file)
opened_files(num_files)%unit = unit
else
unit = opened_files(nfile)%unit
endif
call mpp_get_info(unit,ndim,nvar,natt,ntime)
if (ntime < 1) then
write(msg,'(a15,a,a58)') 'external field ',trim(fieldname),&
' does not have an associated record dimension (REQUIRED) '
call mpp_error(FATAL,trim(msg))
endif
allocate(global_atts(natt))
call mpp_get_atts(unit, global_atts)
allocate(axes(ndim))
call mpp_get_axes(unit, axes, time_axis)
allocate(flds(nvar))
call mpp_get_fields(unit,flds)
allocate(tstamp(ntime),tstart(ntime),tend(ntime),tavg(ntime))
call mpp_get_times(unit,tstamp)
transpose_xy = .false.
isdata=1; iedata=1; jsdata=1; jedata=1
gxsize=1; gysize=1
siz_in = 1
if (PRESENT(domain)) then
call mpp_get_compute_domain(domain,iscomp,iecomp,jscomp,jecomp)
nx = iecomp-iscomp+1; ny = jecomp-jscomp+1
call mpp_get_data_domain(domain,isdata,iedata,jsdata,jedata,dxsize,dxsize_max,dysize,dysize_max)
call mpp_get_global_domain(domain,isglobal,ieglobal,jsglobal,jeglobal,gxsize,gxsize_max,gysize,gysize_max)
elseif(use_comp_domain1) then
call mpp_error(FATAL,"init_external_field:"//&
" use_comp_domain=true but domain is not present")
endif
init_external_field = -1
nfields_orig = num_fields
do i=1,nvar
call mpp_get_atts(flds(i),name=name,units=fld_units,ndim=ndim,siz=siz_in)
call mpp_get_tavg_info(unit,flds(i),flds,tstamp,tstart,tend,tavg)
call mpp_get_atts(flds(i),missing=missing)
! why does it convert case of the field name?
if (trim(lowercase(name)) /= trim(lowercase(fieldname))) cycle
if (verb) write(outunit,*) 'found field ',trim(fieldname), ' in file !!'
num_fields = num_fields + 1
if(num_fields > max_fields) then
!--- z1l: For the case of multiple thread, realoc_fields will cause memory leak.
!--- If multiple threads are working on field A. One of the thread finished first and
!--- begin to work on field B, the realloc_files will cause problem for
!--- other threads are working on the field A.
!call realloc_fields(size(field)*2)
call mpp_error(FATAL, "time_interp_external: num_fields is greater than max_fields, "// &
"increase time_interp_external_nml max_fields")
endif
init_external_field = num_fields
field(num_fields)%unit = unit
field(num_fields)%name = trim(name)
field(num_fields)%units = trim(fld_units)
field(num_fields)%field = flds(i)
field(num_fields)%isc = 1
field(num_fields)%iec = 1
field(num_fields)%jsc = 1
field(num_fields)%jec = 1
field(num_fields)%region_type = NO_REGION
field(num_fields)%is_region = 0
field(num_fields)%ie_region = -1
field(num_fields)%js_region = 0
field(num_fields)%je_region = -1
if (PRESENT(domain)) then
field(num_fields)%domain_present = .true.
field(num_fields)%domain = domain
field(num_fields)%isc=iscomp;field(num_fields)%iec = iecomp
field(num_fields)%jsc=jscomp;field(num_fields)%jec = jecomp
else
field(num_fields)%domain_present = .false.
endif
call mpp_get_atts(flds(i),valid=field(num_fields)%valid )
allocate(fld_axes(ndim))
call mpp_get_atts(flds(i),axes=fld_axes)
if (ndim > 4) call mpp_error(FATAL, &
'invalid array rank <=4d fields supported')
field(num_fields)%siz = 1
field(num_fields)%ndim = ndim
field(num_fields)%tdim = 4
field(num_fields)%missing = missing
do j=1,field(num_fields)%ndim
cart = 'N'
call get_axis_cart(fld_axes(j), cart)
call mpp_get_atts(fld_axes(j),len=len)
if (cart == 'N' .and. .not. ignore_axatts) then
write(msg,'(a,"/",a)') trim(file),trim(fieldname)
call mpp_error(FATAL,'file/field '//trim(msg)// &
' couldnt recognize axis atts in time_interp_external')
else if (cart == 'N' .and. ignore_axatts) then
cart = cart_dir(j)
endif
select case (cart)
case ('X')
if (j.eq.2) transpose_xy = .true.
if (.not.PRESENT(domain) .and. .not.PRESENT(override)) then
isdata=1;iedata=len
iscomp=1;iecomp=len
gxsize = len
dxsize = len
field(num_fields)%isc=iscomp;field(num_fields)%iec=iecomp
elseif (PRESENT(override)) then
gxsize = len
if (PRESENT(axis_sizes)) axis_sizes(1) = len
endif
field(num_fields)%axes(1) = fld_axes(j)
if(use_comp_domain1) then
field(num_fields)%siz(1) = nx
else
field(num_fields)%siz(1) = dxsize
endif
if (len /= gxsize) then
write(msg,'(a,"/",a)') trim(file),trim(fieldname)
call mpp_error(FATAL,'time_interp_ext, file/field '//trim(msg)//' x dim doesnt match model')
endif
case ('Y')
field(num_fields)%axes(2) = fld_axes(j)
if (.not.PRESENT(domain) .and. .not.PRESENT(override)) then
jsdata=1;jedata=len
jscomp=1;jecomp=len
gysize = len
dysize = len
field(num_fields)%jsc=jscomp;field(num_fields)%jec=jecomp
elseif (PRESENT(override)) then
gysize = len
if (PRESENT(axis_sizes)) axis_sizes(2) = len
endif
if(use_comp_domain1) then
field(num_fields)%siz(2) = ny
else
field(num_fields)%siz(2) = dysize
endif
if (len /= gysize) then
write(msg,'(a,"/",a)') trim(file),trim(fieldname)
call mpp_error(FATAL,'time_interp_ext, file/field '//trim(msg)//' y dim doesnt match model')
endif
case ('Z')
field(num_fields)%axes(3) = fld_axes(j)
field(num_fields)%siz(3) = siz_in(3)
case ('T')
field(num_fields)%axes(4) = fld_axes(j)
field(num_fields)%siz(4) = ntime
field(num_fields)%tdim = j
end select
enddo
siz = field(num_fields)%siz
if (PRESENT(axis_centers)) then
axis_centers = field(num_fields)%axes
endif
if (PRESENT(axis_sizes) .and. .not.PRESENT(override)) then
axis_sizes = field(num_fields)%siz
endif
deallocate(fld_axes)
if (verb) write(outunit,'(a,4i6)') 'field x,y,z,t local size= ',siz
if (verb) write(outunit,*) 'field contains data in units = ',trim(field(num_fields)%units)
if (transpose_xy) call mpp_error(FATAL,'axis ordering not supported')
if (num_io_buffers .le. 1) call mpp_error(FATAL,'time_interp_ext:num_io_buffers should be at least 2')
nbuf = min(num_io_buffers,siz(4))
field(num_fields)%numwindows = numwindows
allocate(field(num_fields)%need_compute(nbuf, numwindows))
field(num_fields)%need_compute = .true.
allocate(field(num_fields)%data(isdata:iedata,jsdata:jedata,siz(3),nbuf),&
field(num_fields)%mask(isdata:iedata,jsdata:jedata,siz(3),nbuf) )
field(num_fields)%mask = .false.
field(num_fields)%data = 0.0
slope=1.0;intercept=0.0
! if (units /= 'same') call convert_units(trim(field(num_fields)%units),trim(units),slope,intercept)
! if (verb.and.units /= 'same') then
! write(outunit,*) 'attempting to convert data to units = ',trim(units)
! write(outunit,'(a,f8.3,a,f8.3)') 'factor = ',slope,' offset= ',intercept
! endif
field(num_fields)%slope = slope
field(num_fields)%intercept = intercept
allocate(field(num_fields)%ibuf(nbuf))
field(num_fields)%ibuf = -1
field(num_fields)%nbuf = 0 ! initialize buffer number so that first reading fills data(:,:,:,1)
if(PRESENT(override)) then
field(num_fields)%is_src = 1
field(num_fields)%ie_src = gxsize
field(num_fields)%js_src = 1
field(num_fields)%je_src = gysize
allocate(field(num_fields)%src_data(gxsize,gysize,siz(3),nbuf))
else
field(num_fields)%is_src = isdata
field(num_fields)%ie_src = iedata
field(num_fields)%js_src = jsdata
field(num_fields)%je_src = jedata
allocate(field(num_fields)%src_data(isdata:iedata,jsdata:jedata,siz(3),nbuf))
endif
allocate(field(num_fields)%time(ntime))
allocate(field(num_fields)%period(ntime))
allocate(field(num_fields)%start_time(ntime))
allocate(field(num_fields)%end_time(ntime))
call mpp_get_atts(time_axis,units=units,calendar=calendar_type)
do j=1,ntime
field(num_fields)%time(j) = get_cal_time(tstamp(j),trim(units),trim(calendar_type),permit_calendar_conversion)
field(num_fields)%start_time(j) = get_cal_time(tstart(j),trim(units),trim(calendar_type),permit_calendar_conversion)
field(num_fields)%end_time(j) = get_cal_time( tend(j),trim(units),trim(calendar_type),permit_calendar_conversion)
enddo
if (field(num_fields)%modulo_time) then
call set_time_modulo(field(num_fields)%Time)
call set_time_modulo(field(num_fields)%start_time)
call set_time_modulo(field(num_fields)%end_time)
endif
if(present(correct_leap_year_inconsistency)) then
field(num_fields)%correct_leap_year_inconsistency = correct_leap_year_inconsistency
else
field(num_fields)%correct_leap_year_inconsistency = .false.
endif
if(get_axis_modulo_times(time_axis, timebeg, timeend)) then
if(get_calendar_type() == NO_CALENDAR) then
field(num_fields)%modulo_time_beg = set_time(timebeg)
field(num_fields)%modulo_time_end = set_time(timeend)
else
field(num_fields)%modulo_time_beg = set_date(timebeg)
field(num_fields)%modulo_time_end = set_date(timeend)
endif
field(num_fields)%have_modulo_times = .true.
else
field(num_fields)%have_modulo_times = .false.
endif
if(ntime == 1) then
call mpp_error(NOTE, 'time_interp_external_mod: file '//trim(file)//' has only one time level')
else
do j= 1, ntime
field(num_fields)%period(j) = field(num_fields)%end_time(j)-field(num_fields)%start_time(j)
if (field(num_fields)%period(j) > set_time(0,0)) then
call get_time(field(num_fields)%period(j), sec, day)
sec = sec/2+mod(day,2)*43200
day = day/2
field(num_fields)%time(j) = field(num_fields)%start_time(j)+&
set_time(sec,day)
else
if (j > 1 .and. j < ntime) then
tdiff = field(num_fields)%time(j+1) - field(num_fields)%time(j-1)
call get_time(tdiff, sec, day)
sec = sec/2+mod(day,2)*43200
day = day/2
field(num_fields)%period(j) = set_time(sec,day)
sec = sec/2+mod(day,2)*43200
day = day/2
field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
elseif ( j == 1) then
tdiff = field(num_fields)%time(2) - field(num_fields)%time(1)
call get_time(tdiff, sec, day)
field(num_fields)%period(j) = set_time(sec,day)
sec = sec/2+mod(day,2)*43200
day = day/2
field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
else
tdiff = field(num_fields)%time(ntime) - field(num_fields)%time(ntime-1)
call get_time(tdiff, sec, day)
field(num_fields)%period(j) = set_time(sec,day)
sec = sec/2+mod(day,2)*43200
day = day/2
field(num_fields)%start_time(j) = field(num_fields)%time(j) - set_time(sec,day)
field(num_fields)%end_time(j) = field(num_fields)%time(j) + set_time(sec,day)
endif
endif
enddo
endif
do j=1,ntime-1
if (field(num_fields)%time(j) >= field(num_fields)%time(j+1)) then
write(msg,'(A,i20)') "times not monotonically increasing. Filename: " &
//TRIM(file)//" field: "//TRIM(fieldname)//" timeslice: ", j
call mpp_error(FATAL, TRIM(msg))
endif
enddo
field(num_fields)%modulo_time = get_axis_modulo(time_axis)
if (verb) then
if (field(num_fields)%modulo_time) write(outunit,*) 'data are being treated as modulo in time'
do j= 1, ntime
write(outunit,*) 'time index, ', j
call get_date(field(num_fields)%start_time(j),yr,mon,day,hr,minu,sec)
write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
'start time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
call get_date(field(num_fields)%time(j),yr,mon,day,hr,minu,sec)
write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
'mid time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
call get_date(field(num_fields)%end_time(j),yr,mon,day,hr,minu,sec)
write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
'end time: yyyy/mm/dd hh:mm:ss= ',yr,'/',mon,'/',day,hr,':',minu,':',sec
enddo
end if
enddo
if (num_fields == nfields_orig) then
if (present(ierr)) then
ierr = ERR_FIELD_NOT_FOUND
else
call mpp_error(FATAL,'external field "'//trim(fieldname)//'" not found in file "'//trim(file)//'"')
endif
endif
deallocate(global_atts)
deallocate(axes)
deallocate(flds)
deallocate(tstamp, tstart, tend, tavg)
return
end function init_external_field
!</FUNCTION> NAME="init_external_field"
!> @brief 2D time interpolation for @ref time_interp_external
subroutine time_interp_external_2d(index, time, data_in, interp, verbose,horz_interp, mask_out, &
is_in, ie_in, js_in, je_in, window_id)
integer, intent(in) :: index
type(time_type), intent(in) :: time
real, dimension(:,:), intent(inout) :: data_in
integer, intent(in), optional :: interp
logical, intent(in), optional :: verbose
type(horiz_interp_type),intent(in), optional :: horz_interp
logical, dimension(:,:), intent(out), optional :: mask_out ! set to true where output data is valid
integer, intent(in), optional :: is_in, ie_in, js_in, je_in
integer, intent(in), optional :: window_id
real , dimension(size(data_in,1), size(data_in,2), 1) :: data_out
logical, dimension(size(data_in,1), size(data_in,2), 1) :: mask3d
data_out(:,:,1) = data_in(:,:) ! fill initial values for the portions of array that are not touched by 3d routine
call time_interp_external_3d(index, time, data_out, interp, verbose, horz_interp, mask3d, &
is_in=is_in,ie_in=ie_in,js_in=js_in,je_in=je_in,window_id=window_id)
data_in(:,:) = data_out(:,:,1)
if (PRESENT(mask_out)) mask_out(:,:) = mask3d(:,:,1)
return
end subroutine time_interp_external_2d
!<SUBROUTINE NAME="time_interp_external" >
!
!<DESCRIPTION>
! Provide data from external file interpolated to current model time.
! Data may be local to current processor or global, depending on
! "init_external_field" flags.
!</DESCRIPTION>
!
!<IN NAME="index" TYPE="integer">
! index of external field from previous call to init_external_field
!</IN>
!<IN NAME="time" TYPE="time_manager_mod:time_type">
! target time for data
!</IN>
!<INOUT NAME="data" TYPE="real" DIM="(:,:),(:,:,:)">
! global or local data array
!</INOUT>
!<IN NAME="interp" TYPE="integer">
! time_interp_external defined interpolation method (optional). Currently this module only supports
! LINEAR_TIME_INTERP.
!</IN>
!<IN NAME="verbose" TYPE="logical">
! verbose flag for debugging (optional).
!</IN>
!> @brief 3D time interpolation for @ref time_interp_external
subroutine time_interp_external_3d(index, time, data, interp,verbose,horz_interp, mask_out, is_in, ie_in, js_in, je_in, window_id)
integer, intent(in) :: index
type(time_type), intent(in) :: time
real, dimension(:,:,:), intent(inout) :: data
integer, intent(in), optional :: interp
logical, intent(in), optional :: verbose
type(horiz_interp_type), intent(in), optional :: horz_interp
logical, dimension(:,:,:), intent(out), optional :: mask_out ! set to true where output data is valid
integer, intent(in), optional :: is_in, ie_in, js_in, je_in
integer, intent(in), optional :: window_id
integer :: nx, ny, nz, interp_method, t1, t2
integer :: i1, i2, isc, iec, jsc, jec, mod_time
integer :: yy, mm, dd, hh, min, ss
character(len=256) :: err_msg, filename
integer :: isw, iew, jsw, jew, nxw, nyw
! these are boundaries of the updated portion of the "data" argument
! they are calculated using sizes of the "data" and isc,iec,jsc,jsc
! fileds from respective input field, to center the updated portion
! in the output array
real :: w1,w2
logical :: verb
character(len=16) :: message1, message2
nx = size(data,1)
ny = size(data,2)
nz = size(data,3)
interp_method = LINEAR_TIME_INTERP
if (PRESENT(interp)) interp_method = interp
verb=.false.
if (PRESENT(verbose)) verb=verbose
if (debug_this_module) verb = .true.
if (index < 1.or.index > num_fields) &
call mpp_error(FATAL,'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize')
isc=field(index)%isc;iec=field(index)%iec
jsc=field(index)%jsc;jec=field(index)%jec
if( field(index)%numwindows == 1 ) then
nxw = iec-isc+1
nyw = jec-jsc+1
else
if( .not. present(is_in) .or. .not. present(ie_in) .or. .not. present(js_in) .or. .not. present(je_in) ) then
call mpp_error(FATAL, 'time_interp_external: is_in, ie_in, js_in and je_in must be present ' // &
'when numwindows > 1, field='//trim(field(index)%name))
endif
nxw = ie_in - is_in + 1
nyw = je_in - js_in + 1
isc = isc + is_in - 1
iec = isc + ie_in - is_in
jsc = jsc + js_in - 1
jec = jsc + je_in - js_in
endif
isw = (nx-nxw)/2+1; iew = isw+nxw-1
jsw = (ny-nyw)/2+1; jew = jsw+nyw-1
if (nx < nxw .or. ny < nyw .or. nz < field(index)%siz(3)) then
write(message1,'(i6,2i5)') nx,ny,nz
call mpp_error(FATAL,'field '//trim(field(index)%name)//' Array size mismatch in time_interp_external.'// &
' Array "data" is too small. shape(data)='//message1)
endif
if(PRESENT(mask_out)) then
if (size(mask_out,1) /= nx .or. size(mask_out,2) /= ny .or. size(mask_out,3) /= nz) then
write(message1,'(i6,2i5)') nx,ny,nz
write(message2,'(i6,2i5)') size(mask_out,1),size(mask_out,2),size(mask_out,3)
call mpp_error(FATAL,'field '//trim(field(index)%name)//' array size mismatch in time_interp_external.'// &
' Shape of array "mask_out" does not match that of array "data".'// &
' shape(data)='//message1//' shape(mask_out)='//message2)
endif
endif
if (field(index)%siz(4) == 1) then
! only one record in the file => time-independent field
call load_record(field(index),1,horz_interp, is_in, ie_in ,js_in, je_in,window_id)
i1 = find_buf_index(1,field(index)%ibuf)
if( field(index)%region_type == NO_REGION ) then
where(field(index)%mask(isc:iec,jsc:jec,:,i1))
data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)
elsewhere
! data(isw:iew,jsw:jew,:) = time_interp_missing !field(index)%missing? Balaji
data(isw:iew,jsw:jew,:) = field(index)%missing
end where
else
where(field(index)%mask(isc:iec,jsc:jec,:,i1))
data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)
end where
endif
if(PRESENT(mask_out)) &
mask_out(isw:iew,jsw:jew,:) = field(index)%mask(isc:iec,jsc:jec,:,i1)
else
if(field(index)%have_modulo_times) then
call time_interp(time,field(index)%modulo_time_beg, field(index)%modulo_time_end, field(index)%time(:), &
w2, t1, t2, field(index)%correct_leap_year_inconsistency, err_msg=err_msg)
if(err_msg .NE. '') then
filename = mpp_get_file_name(field(index)%unit)
call mpp_error(FATAL,"time_interp_external 1: "//trim(err_msg)//&
",file="//trim(filename)//",field="//trim(field(index)%name) )
endif
else
if(field(index)%modulo_time) then
mod_time=1
else
mod_time=0
endif
call time_interp(time,field(index)%time(:),w2,t1,t2,modtime=mod_time, err_msg=err_msg)
if(err_msg .NE. '') then
filename = mpp_get_file_name(field(index)%unit)
call mpp_error(FATAL,"time_interp_external 2: "//trim(err_msg)//&
",file="//trim(filename)//",field="//trim(field(index)%name) )
endif
endif
w1 = 1.0-w2
if (verb) then
call get_date(time,yy,mm,dd,hh,min,ss)
write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss
write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
endif
call load_record(field(index),t1,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
call load_record(field(index),t2,horz_interp, is_in, ie_in ,js_in, je_in, window_id)
i1 = find_buf_index(t1,field(index)%ibuf)
i2 = find_buf_index(t2,field(index)%ibuf)
if(i1<0.or.i2<0) &
call mpp_error(FATAL,'time_interp_external : records were not loaded correctly in memory')
if (verb) then
write(outunit,*) 'ibuf= ',field(index)%ibuf
write(outunit,*) 'i1,i2= ',i1, i2
endif
if( field(index)%region_type == NO_REGION ) then
where(field(index)%mask(isc:iec,jsc:jec,:,i1).and.field(index)%mask(isc:iec,jsc:jec,:,i2))
data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)*w1 + &
field(index)%data(isc:iec,jsc:jec,:,i2)*w2
elsewhere
! data(isw:iew,jsw:jew,:) = time_interp_missing !field(index)%missing? Balaji
data(isw:iew,jsw:jew,:) = field(index)%missing
end where
else
where(field(index)%mask(isc:iec,jsc:jec,:,i1).and.field(index)%mask(isc:iec,jsc:jec,:,i2))
data(isw:iew,jsw:jew,:) = field(index)%data(isc:iec,jsc:jec,:,i1)*w1 + &
field(index)%data(isc:iec,jsc:jec,:,i2)*w2
end where
endif
if(PRESENT(mask_out)) &
mask_out(isw:iew,jsw:jew,:) = &
field(index)%mask(isc:iec,jsc:jec,:,i1).and.&
field(index)%mask(isc:iec,jsc:jec,:,i2)
endif
end subroutine time_interp_external_3d
!</SUBROUTINE> NAME="time_interp_external"
!> @brief Scalar time interpolation for @ref time_interp_external
subroutine time_interp_external_0d(index, time, data, verbose)
integer, intent(in) :: index
type(time_type), intent(in) :: time
real, intent(inout) :: data
logical, intent(in), optional :: verbose
integer :: t1, t2
integer :: i1, i2, mod_time
integer :: yy, mm, dd, hh, min, ss
character(len=256) :: err_msg, filename
real :: w1,w2
logical :: verb
verb=.false.
if (PRESENT(verbose)) verb=verbose
if (debug_this_module) verb = .true.
if (index < 1.or.index > num_fields) &
call mpp_error(FATAL,'invalid index in call to time_interp_ext -- field was not initialized or failed to initialize')
if (field(index)%siz(4) == 1) then
! only one record in the file => time-independent field
call load_record_0d(field(index),1)
i1 = find_buf_index(1,field(index)%ibuf)
data = field(index)%data(1,1,1,i1)
else
if(field(index)%have_modulo_times) then
call time_interp(time,field(index)%modulo_time_beg, field(index)%modulo_time_end, field(index)%time(:), &
w2, t1, t2, field(index)%correct_leap_year_inconsistency, err_msg=err_msg)
if(err_msg .NE. '') then
filename = mpp_get_file_name(field(index)%unit)
call mpp_error(FATAL,"time_interp_external 3:"//trim(err_msg)//&
",file="//trim(filename)//",field="//trim(field(index)%name) )
endif
else
if(field(index)%modulo_time) then
mod_time=1
else
mod_time=0
endif
call time_interp(time,field(index)%time(:),w2,t1,t2,modtime=mod_time, err_msg=err_msg)
if(err_msg .NE. '') then
filename = mpp_get_file_name(field(index)%unit)
call mpp_error(FATAL,"time_interp_external 4:"//trim(err_msg)// &
",file="//trim(filename)//",field="//trim(field(index)%name) )
endif
endif
w1 = 1.0-w2
if (verb) then
call get_date(time,yy,mm,dd,hh,min,ss)
write(outunit,'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
'target time yyyy/mm/dd hh:mm:ss= ',yy,'/',mm,'/',dd,hh,':',min,':',ss
write(outunit,*) 't1, t2, w1, w2= ', t1, t2, w1, w2
endif
call load_record_0d(field(index),t1)
call load_record_0d(field(index),t2)
i1 = find_buf_index(t1,field(index)%ibuf)
i2 = find_buf_index(t2,field(index)%ibuf)
if(i1<0.or.i2<0) &
call mpp_error(FATAL,'time_interp_external : records were not loaded correctly in memory')
data = field(index)%data(1,1,1,i1)*w1 + field(index)%data(1,1,1,i2)*w2
if (verb) then
write(outunit,*) 'ibuf= ',field(index)%ibuf
write(outunit,*) 'i1,i2= ',i1, i2
endif
endif
end subroutine time_interp_external_0d
subroutine set_time_modulo(Time)
type(time_type), intent(inout), dimension(:) :: Time
integer :: ntime, n
integer :: yr, mon, dy, hr, minu, sec
ntime = size(Time(:))
do n = 1, ntime
call get_date(Time(n), yr, mon, dy, hr, minu, sec)
yr = modulo_year
Time(n) = set_date(yr, mon, dy, hr, minu, sec)
enddo
end subroutine set_time_modulo