35 module time_interp_external2_mod
 
   44   use platform_mod, 
only : r8_kind, r4_kind
 
   54   use time_interp_mod, 
only : 
time_interp, time_interp_init
 
   57   use platform_mod, 
only: r8_kind, fms_path_len, fms_file_len
 
   58   use horiz_interp_mod, 
only : 
horiz_interp, horiz_interp_type
 
   59   use horiz_interp_type_mod, 
only : conserve
 
   60   use fms2_io_mod,      
only : valid_t, fmsnetcdfdomainfile_t, 
open_file, get_unlimited_dimension_name, &
 
   61                                variable_att_exists, fmsnetcdffile_t, &
 
   62                                variable_exists, get_valid, get_variable_num_dimensions, 
read_data, &
 
   63                                is_valid, 
close_file, get_dimension_size, get_variable_dimension_names, &
 
   64                                get_variable_size, get_time_calendar, get_variable_missing, get_variable_units
 
   70 #include<file_version.h> 
   72   integer, 
parameter, 
public  :: NO_REGION=0, inside_region=1, outside_region=2
 
   73   integer, 
parameter, 
private :: modulo_year= 0001
 
   74   integer, 
parameter, 
private :: LINEAR_TIME_INTERP = 1 
 
   75   integer, 
parameter, 
public  :: SUCCESS = 0, err_field_not_found = 1
 
   76   real(r8_kind),    
parameter, 
private :: DEFAULT_MISSING_VALUE = -1e20_r8_kind
 
   77   integer,            
private :: max_fields = 100, max_files= 40
 
   78   integer, 
private :: num_fields = 0, num_files=0
 
   80   integer, 
private :: num_io_buffers = 2 
 
   81   logical, 
private :: module_initialized = .false.
 
   82   logical, 
private :: debug_this_module = .false.
 
   87   public get_external_fileobj
 
   97         type(fmsnetcdffile_t), 
pointer :: fileobj=>null() 
 
   98         character(len=128) :: name, units
 
   99         integer :: siz(4), ndim
 
  100         character(len=32) :: axisname(4)
 
  103         type(
time_type), 
dimension(:), 
pointer :: start_time =>null(), end_time =>null()
 
  104         type(
time_type), 
dimension(:), 
pointer :: period =>null()
 
  105         logical :: modulo_time
 
  106         real(r8_kind), 
dimension(:,:,:,:), 
pointer :: domain_data =>null() 
 
  107         logical, 
dimension(:,:,:,:), 
pointer :: mask =>null() 
 
  108         integer, 
dimension(:), 
pointer :: ibuf  =>null() 
 
  109         real(r8_kind), 
dimension(:,:,:,:), 
pointer :: src_data  =>null() 
 
  110         type(valid_t) :: valid 
 
  112         logical :: domain_present
 
  113         real(r8_kind) :: slope, intercept
 
  114         integer :: isc,iec,jsc,jec
 
  115         type(
time_type) :: modulo_time_beg, modulo_time_end
 
  116         logical :: have_modulo_times, correct_leap_year_inconsistency
 
  117         integer :: region_type
 
  118         integer :: is_region, ie_region, js_region, je_region
 
  119         integer :: is_src, ie_src, js_src, je_src
 
  121         integer :: numwindows
 
  122         logical, 
dimension(:,:), 
pointer :: need_compute=>null()
 
  123         real(r8_kind)    :: missing 
 
  129         character(len=FMS_FILE_LEN) :: filename = 
'' 
  130         type(FmsNetcdfFile_t), 
pointer :: fileobj => null()
 
  146       module procedure time_interp_external_0d_r4
 
  147       module procedure time_interp_external_1d_r4
 
  148       module procedure time_interp_external_2d_r4
 
  149       module procedure time_interp_external_3d_r4
 
  150       module procedure time_interp_external_0d_r8
 
  151       module procedure time_interp_external_1d_r8
 
  152       module procedure time_interp_external_2d_r8
 
  153       module procedure time_interp_external_3d_r8
 
  157      module procedure time_interp_external_bridge_0d_r4
 
  158      module procedure time_interp_external_bridge_1d_r4
 
  159      module procedure time_interp_external_bridge_2d_r4
 
  160      module procedure time_interp_external_bridge_3d_r4
 
  161      module procedure time_interp_external_bridge_0d_r8
 
  162      module procedure time_interp_external_bridge_1d_r8
 
  163      module procedure time_interp_external_bridge_2d_r8
 
  164      module procedure time_interp_external_bridge_3d_r8
 
  172   type(ext_fieldtype), 
save, 
private, 
pointer :: loaded_fields(:) => null()
 
  173   type(filetype),      
save, 
private, 
pointer :: opened_files(:) => null()
 
  175   real(r8_kind), 
private, 
parameter :: time_interp_missing=-1e99_r8_kind
 
  187       integer :: io_status, logunit, ierr
 
  189       namelist /time_interp_external_nml/ num_io_buffers, debug_this_module, &
 
  190                                           max_fields, max_files
 
  194       if(module_initialized) 
return 
  198       call write_version_number(
"TIME_INTERP_EXTERNAL_MOD", version)
 
  200       read (input_nml_file, time_interp_external_nml, iostat=io_status)
 
  203       write(logunit,time_interp_external_nml)
 
  207       module_initialized = .true.
 
  209       call time_interp_init()
 
  275          verbose,axis_names, axis_sizes,override,correct_leap_year_inconsistency,&
 
  276          permit_calendar_conversion,use_comp_domain,ierr, nwindows, ignore_axis_atts, ongrid )
 
  278       character(len=*), 
intent(in)            :: file,fieldname
 
  279       logical, 
intent(in), 
optional           :: verbose
 
  280       character(len=*), 
intent(in), 
optional  :: desired_units
 
  281       type(
domain2d), 
intent(in), 
optional    :: domain
 
  282       integer, 
intent(inout), 
optional        :: axis_sizes(4)
 
  283       character(len=*), 
intent(out), 
optional :: axis_names(4)
 
  284       logical, 
intent(in), 
optional           :: override, correct_leap_year_inconsistency,&
 
  285                                                  permit_calendar_conversion,use_comp_domain
 
  286       integer,          
intent(out), 
optional :: ierr
 
  287       integer,          
intent(in),  
optional :: nwindows
 
  288       logical, 
optional                       :: ignore_axis_atts
 
  289       logical, 
optional                       :: ongrid
 
  291       logical :: ongrid_local
 
  295       real(r8_kind) :: slope, intercept
 
  296       integer :: ndim,ntime,i,j
 
  297       integer :: iscomp,iecomp,jscomp,jecomp,isglobal,ieglobal,jsglobal,jeglobal
 
  298       integer :: isdata,iedata,jsdata,jedata, dxsize, dysize,dxsize_max,dysize_max
 
  299       logical :: verb, transpose_xy,use_comp_domain1
 
  300       real(r8_kind), 
dimension(:), 
allocatable :: tstamp, tstart, tend, tavg
 
  301       character(len=1) :: cart
 
  302       character(len=1), 
dimension(4) :: cart_dir
 
  303       character(len=128) :: units, fld_units
 
  304       character(len=128) :: msg, calendar_type, timebeg, timeend
 
  305       character(len=128) :: timename, timeunits
 
  306       character(len=128), 
allocatable :: axisname(:)
 
  307       integer,            
allocatable :: axislen(:)
 
  308       integer :: siz(4), siz_in(4), gxsize, gysize,gxsize_max, gysize_max
 
  310       integer :: yr, mon, day, hr, minu, sec
 
  311       integer :: len, nfile, nfields_orig, nbuf, nx,ny
 
  312       integer :: numwindows
 
  313       logical :: ignore_axatts
 
  314       logical :: have_modulo_time
 
  315       type(fmsnetcdffile_t), 
pointer :: fileobj=>null()
 
  316       integer, 
dimension(:), 
allocatable :: pes
 
  318       if (.not. module_initialized) 
call mpp_error(fatal,
'Must call time_interp_external_init first')
 
  319       if(
present(ierr)) ierr = success
 
  320       ignore_axatts=.false.
 
  321       cart_dir(1)=
'X';cart_dir(2)=
'Y';cart_dir(3)=
'Z';cart_dir(4)=
'T' 
  322       if(
present(ignore_axis_atts)) ignore_axatts = ignore_axis_atts
 
  323       use_comp_domain1 = .false.
 
  324       if(
PRESENT(use_comp_domain)) use_comp_domain1 = use_comp_domain
 
  326       if (
PRESENT(verbose)) verb=verbose
 
  327       if (debug_this_module) verb = .true.
 
  329       if(
present(nwindows)) numwindows = nwindows
 
  332       if (
PRESENT(desired_units)) 
then 
  333           units = desired_units
 
  334           call mpp_error(fatal,
'==> Unit conversion via time_interp_external & 
  335                &has been temporarily deprecated.  Previous versions of& 
  336                &this module used udunits_mod to perform unit conversion.& 
  337                &  Udunits_mod is in the process of being replaced since & 
  338                &there were portability issues associated with this code.& 
  339                & Please remove the desired_units argument from calls to & 
  344          if(trim(opened_files(i)%filename) == trim(file)) 
then 
  350          num_files = num_files + 1
 
  351          if(num_files > max_files) 
then  
  357             call mpp_error(fatal, 
"time_interp_external: num_files is greater than max_files, "// &
 
  358                                       "increase time_interp_external_nml max_files")
 
  360          opened_files(num_files)%filename = trim(file)
 
  361          allocate(opened_files(num_files)%fileobj)
 
  362          fileobj => opened_files(num_files)%fileobj
 
  364          if(.not. 
open_file(opened_files(num_files)%fileobj, trim(file), 
'read')) &
 
  365              call mpp_error(fatal, 
'time_interp_external_mod: Error in opening file '//trim(file))
 
  367          fileobj => opened_files(nfile)%fileobj
 
  370       call get_unlimited_dimension_name(fileobj, timename)
 
  371       call get_dimension_size(fileobj, timename, ntime)
 
  374           write(msg,
'(a15,a,a58)') 
'external field ',trim(fieldname),&
 
  375            ' does not have an associated record dimension (REQUIRED) ' 
  380       call get_time_calendar(fileobj, timename, calendar_type)
 
  383       call get_variable_units(fileobj, timename, timeunits)
 
  389       call mpp_get_current_pelist(pes)
 
  390       allocate(tstamp(ntime),tstart(ntime),tend(ntime),tavg(ntime))
 
  393       if (mpp_root_pe() .eq. 
mpp_pe()) 
call read_data(fileobj, timename, tstamp)
 
  394       call mpp_broadcast(tstamp, 
size(tstamp), mpp_root_pe(), pelist=pes)
 
  397       transpose_xy = .false.
 
  398       isdata=1; iedata=1; jsdata=1; jedata=1
 
  402       if (
PRESENT(domain)) 
then 
  404          nx = iecomp-iscomp+1; ny = jecomp-jscomp+1
 
  405          call mpp_get_data_domain(domain,isdata,iedata,jsdata,jedata,dxsize,dxsize_max,dysize,dysize_max)
 
  406          call mpp_get_global_domain(domain,isglobal,ieglobal,jsglobal,jeglobal,gxsize,gxsize_max,gysize,gysize_max)
 
  407          ongrid_local = .false.
 
  408          if (
present(ongrid)) ongrid_local = ongrid
 
  411          if (ongrid_local) 
then 
  417       elseif(use_comp_domain1) 
then 
  418          call mpp_error(fatal,
"init_external_field:"//&
 
  419               " use_comp_domain=true but domain is not present")
 
  423       nfields_orig = num_fields
 
  425       if (.not. variable_exists(fileobj, fieldname) ) 
then 
  426         if (
present(ierr)) 
then 
  427            ierr = err_field_not_found
 
  430            call mpp_error(fatal,
'external field "'//trim(fieldname)//
'" not found in file "'//trim(file)//
'"')
 
  437       if(variable_att_exists(fileobj, fieldname, 
'time_avg_info')) 
then 
  438         if(variable_exists(fileobj, 
'average_T1')) 
call read_data(fileobj, 
'average_T1', tstart)
 
  439         if(variable_exists(fileobj, 
'average_T2')) 
call read_data(fileobj, 
'average_T2', tend)
 
  440         if(variable_exists(fileobj, 
'average_DT')) 
call read_data(fileobj, 
'average_DT', tavg)
 
  443       num_fields = num_fields + 1
 
  444       if(num_fields > max_fields) 
then 
  450          call mpp_error(fatal, 
"time_interp_external: num_fields is greater than max_fields, "// &
 
  451                                    "increase time_interp_external_nml max_fields")
 
  455       call get_variable_units(fileobj, fieldname, fld_units)
 
  458       loaded_fields(num_fields)%fileobj => fileobj
 
  459       loaded_fields(num_fields)%name = trim(fieldname)
 
  460       loaded_fields(num_fields)%units = trim(fld_units)
 
  461       loaded_fields(num_fields)%isc = 1
 
  462       loaded_fields(num_fields)%iec = 1
 
  463       loaded_fields(num_fields)%jsc = 1
 
  464       loaded_fields(num_fields)%jec = 1
 
  465       loaded_fields(num_fields)%region_type = no_region
 
  466       loaded_fields(num_fields)%is_region   = 0
 
  467       loaded_fields(num_fields)%ie_region   = -1
 
  468       loaded_fields(num_fields)%js_region   = 0
 
  469       loaded_fields(num_fields)%je_region   = -1
 
  470       if (
PRESENT(domain)) 
then 
  471          loaded_fields(num_fields)%domain_present = .true.
 
  472          loaded_fields(num_fields)%domain = domain
 
  473          loaded_fields(num_fields)%isc=iscomp;loaded_fields(num_fields)%iec = iecomp
 
  474          loaded_fields(num_fields)%jsc=jscomp;loaded_fields(num_fields)%jec = jecomp
 
  476          loaded_fields(num_fields)%domain_present = .false.
 
  479       loaded_fields(num_fields)%valid = get_valid(fileobj, fieldname)
 
  480       ndim = get_variable_num_dimensions(fileobj, fieldname)
 
  482            'invalid array rank <=4d fields supported')
 
  484       loaded_fields(num_fields)%siz = 1
 
  485       loaded_fields(num_fields)%ndim = ndim
 
  486       loaded_fields(num_fields)%tdim = 4
 
  488       loaded_fields(num_fields)%missing = get_variable_missing(fileobj, fieldname)
 
  490       allocate(axisname(ndim), axislen(ndim))
 
  492       call get_variable_dimension_names(fileobj, fieldname, axisname)
 
  493       call get_variable_size(fileobj, fieldname, axislen)
 
  494       do j=1,loaded_fields(num_fields)%ndim
 
  497          if (cart == 
'N' .and. .not. ignore_axatts) 
then 
  498             write(msg,
'(a,"/",a)')  trim(file),trim(fieldname)
 
  499             call mpp_error(fatal,
'file/field '//trim(msg)// &
 
  500                  ' couldnt recognize axis atts in time_interp_external')
 
  501          else if (cart == 
'N' .and. ignore_axatts) 
then 
  506             if (j.eq.2) transpose_xy = .true.
 
  507             if (.not.
PRESENT(domain) .and. .not.
PRESENT(override)) 
then 
  512                loaded_fields(num_fields)%isc=iscomp;loaded_fields(num_fields)%iec=iecomp
 
  513             elseif (
PRESENT(override)) 
then 
  515                if (
PRESENT(axis_sizes)) axis_sizes(1) = len
 
  517             loaded_fields(num_fields)%axisname(1) = axisname(j)
 
  518             if(use_comp_domain1) 
then 
  519                loaded_fields(num_fields)%siz(1) = nx
 
  521                loaded_fields(num_fields)%siz(1) = dxsize
 
  523             if (len /= gxsize) 
then 
  524                write(msg,
'(a,"/",a)')  trim(file),trim(fieldname)
 
  525                call mpp_error(fatal,
'time_interp_ext, file/field '//trim(msg)//
' x dim doesnt match model')
 
  528             loaded_fields(num_fields)%axisname(2) = axisname(j)
 
  529             if (.not.
PRESENT(domain) .and. .not.
PRESENT(override)) 
then 
  534                loaded_fields(num_fields)%jsc=jscomp;loaded_fields(num_fields)%jec=jecomp
 
  535             elseif (
PRESENT(override)) 
then 
  537                if (
PRESENT(axis_sizes)) axis_sizes(2) = len
 
  539             if(use_comp_domain1) 
then 
  540                loaded_fields(num_fields)%siz(2) = ny
 
  542                loaded_fields(num_fields)%siz(2) = dysize
 
  544             if (len /= gysize) 
then 
  545                write(msg,
'(a,"/",a)')  trim(file),trim(fieldname)
 
  546                call mpp_error(fatal,
'time_interp_ext, file/field '//trim(msg)//
' y dim doesnt match model')
 
  549             loaded_fields(num_fields)%axisname(3) = axisname(j)
 
  550             loaded_fields(num_fields)%siz(3) = len
 
  552             loaded_fields(num_fields)%axisname(4) = axisname(j)
 
  553             loaded_fields(num_fields)%siz(4) = ntime
 
  554             loaded_fields(num_fields)%tdim   = j
 
  557       siz = loaded_fields(num_fields)%siz
 
  558       if(
PRESENT(axis_names)) axis_names = loaded_fields(num_fields)%axisname
 
  559       if (
PRESENT(axis_sizes) .and. .not.
PRESENT(override)) 
then 
  560          axis_sizes = loaded_fields(num_fields)%siz
 
  563       if (verb) 
write(outunit,
'(a,4i6)') 
'field x,y,z,t local size= ',siz
 
  564       if (verb) 
write(outunit,*) 
'field contains data in units = ',trim(loaded_fields(num_fields)%units)
 
  565       if (transpose_xy) 
call mpp_error(fatal,
'axis ordering not supported')
 
  566       if (num_io_buffers .le. 1) 
call mpp_error(fatal,
'time_interp_ext:num_io_buffers should be at least 2')
 
  567       nbuf = min(num_io_buffers,siz(4))
 
  569       loaded_fields(num_fields)%numwindows = numwindows
 
  570       allocate(loaded_fields(num_fields)%need_compute(nbuf, numwindows))
 
  571       loaded_fields(num_fields)%need_compute = .true.
 
  573       allocate(loaded_fields(num_fields)%domain_data(isdata:iedata,jsdata:jedata,siz(3),nbuf),&
 
  574                loaded_fields(num_fields)%mask(isdata:iedata,jsdata:jedata,siz(3),nbuf) )
 
  575       loaded_fields(num_fields)%mask = .false.
 
  576       loaded_fields(num_fields)%domain_data = 0.0_r8_kind
 
  577          slope=1.0_r8_kind;intercept=0.0_r8_kind
 
  583       loaded_fields(num_fields)%slope = slope
 
  584       loaded_fields(num_fields)%intercept = intercept
 
  585       allocate(loaded_fields(num_fields)%ibuf(nbuf))
 
  586       loaded_fields(num_fields)%ibuf = -1
 
  587       loaded_fields(num_fields)%nbuf =  0 
 
  588       if(
PRESENT(override)) 
then 
  589          loaded_fields(num_fields)%is_src = 1
 
  590          loaded_fields(num_fields)%ie_src = gxsize
 
  591          loaded_fields(num_fields)%js_src = 1
 
  592          loaded_fields(num_fields)%je_src = gysize
 
  593          allocate(loaded_fields(num_fields)%src_data(gxsize,gysize,siz(3),nbuf))
 
  595          loaded_fields(num_fields)%is_src = isdata
 
  596          loaded_fields(num_fields)%ie_src = iedata
 
  597          loaded_fields(num_fields)%js_src = jsdata
 
  598          loaded_fields(num_fields)%je_src = jedata
 
  599          allocate(loaded_fields(num_fields)%src_data(isdata:iedata,jsdata:jedata,siz(3),nbuf))
 
  602       allocate(loaded_fields(num_fields)%time(ntime))
 
  603       allocate(loaded_fields(num_fields)%period(ntime))
 
  604       allocate(loaded_fields(num_fields)%start_time(ntime))
 
  605       allocate(loaded_fields(num_fields)%end_time(ntime))
 
  608          loaded_fields(num_fields)%time(j)       = 
get_cal_time(tstamp(j),trim(timeunits),trim(calendar_type), &
 
  609               & permit_calendar_conversion)
 
  610          loaded_fields(num_fields)%start_time(j) = 
get_cal_time(tstart(j),trim(timeunits),trim(calendar_type), &
 
  611               & permit_calendar_conversion)
 
  612          loaded_fields(num_fields)%end_time(j)   = 
get_cal_time(  tend(j),trim(timeunits),trim(calendar_type), &
 
  613               & permit_calendar_conversion)
 
  616       if (loaded_fields(num_fields)%modulo_time) 
then 
  617          call set_time_modulo(loaded_fields(num_fields)%Time)
 
  618          call set_time_modulo(loaded_fields(num_fields)%start_time)
 
  619          call set_time_modulo(loaded_fields(num_fields)%end_time)
 
  622       if(
present(correct_leap_year_inconsistency)) 
then 
  623         loaded_fields(num_fields)%correct_leap_year_inconsistency = correct_leap_year_inconsistency
 
  625         loaded_fields(num_fields)%correct_leap_year_inconsistency = .false.
 
  628       if(have_modulo_time) 
then 
  630             loaded_fields(num_fields)%modulo_time_beg = 
set_time(timebeg)
 
  631             loaded_fields(num_fields)%modulo_time_end = 
set_time(timeend)
 
  633             loaded_fields(num_fields)%modulo_time_beg = 
set_date(timebeg)
 
  634             loaded_fields(num_fields)%modulo_time_end = 
set_date(timeend)
 
  636          loaded_fields(num_fields)%have_modulo_times = .true.
 
  638          loaded_fields(num_fields)%have_modulo_times = .false.
 
  641          call mpp_error(note, 
'time_interp_external_mod: file '//trim(file)//
'  has only one time level')
 
  644             loaded_fields(num_fields)%period(j) = loaded_fields(num_fields)%end_time(j) &
 
  645                                                 - loaded_fields(num_fields)%start_time(j)
 
  646             if (loaded_fields(num_fields)%period(j) > 
set_time(0,0)) 
then 
  647                call get_time(loaded_fields(num_fields)%period(j), sec, day)
 
  648                sec = sec/2+mod(day,2)*43200
 
  650                loaded_fields(num_fields)%time(j) = loaded_fields(num_fields)%start_time(j)+&
 
  653                if (j > 1 .and. j < ntime) 
then 
  654                   tdiff = loaded_fields(num_fields)%time(j+1) -  loaded_fields(num_fields)%time(j-1)
 
  656                   sec = sec/2+mod(day,2)*43200
 
  658                   loaded_fields(num_fields)%period(j) = 
set_time(sec,day)
 
  659                   sec = sec/2+mod(day,2)*43200
 
  661                   loaded_fields(num_fields)%start_time(j) = loaded_fields(num_fields)%time(j) - 
set_time(sec,day)
 
  662                   loaded_fields(num_fields)%end_time(j) = loaded_fields(num_fields)%time(j) + 
set_time(sec,day)
 
  663                elseif ( j == 1) 
then 
  664                   tdiff = loaded_fields(num_fields)%time(2) -  loaded_fields(num_fields)%time(1)
 
  666                   loaded_fields(num_fields)%period(j) = 
set_time(sec,day)
 
  667                   sec = sec/2+mod(day,2)*43200
 
  669                   loaded_fields(num_fields)%start_time(j) = loaded_fields(num_fields)%time(j) - 
set_time(sec,day)
 
  670                   loaded_fields(num_fields)%end_time(j) = loaded_fields(num_fields)%time(j) + 
set_time(sec,day)
 
  672                   tdiff = loaded_fields(num_fields)%time(ntime) -  loaded_fields(num_fields)%time(ntime-1)
 
  674                   loaded_fields(num_fields)%period(j) = 
set_time(sec,day)
 
  675                   sec = sec/2+mod(day,2)*43200
 
  677                   loaded_fields(num_fields)%start_time(j) = loaded_fields(num_fields)%time(j) - 
set_time(sec,day)
 
  678                   loaded_fields(num_fields)%end_time(j) = loaded_fields(num_fields)%time(j) + 
set_time(sec,day)
 
  685          if (loaded_fields(num_fields)%time(j) >= loaded_fields(num_fields)%time(j+1)) 
then 
  686             write(msg,
'(A,i20)') 
"times not monotonically increasing. Filename: " &
 
  687                  //trim(file)//
"  field:  "//trim(fieldname)//
" timeslice: ", j
 
  692       loaded_fields(num_fields)%modulo_time = 
get_axis_modulo(fileobj, timename)
 
  695          if (loaded_fields(num_fields)%modulo_time) 
write(outunit,*) 
'data are being treated as modulo in time' 
  697             write(outunit,*) 
'time index,  ', j
 
  698             call get_date(loaded_fields(num_fields)%start_time(j),yr,mon,day,hr,minu,sec)
 
  699             write(outunit,
'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
 
  700                  'start time: yyyy/mm/dd hh:mm:ss= ',yr,
'/',mon,
'/',day,hr,
':',minu,
':',sec
 
  701             call get_date(loaded_fields(num_fields)%time(j),yr,mon,day,hr,minu,sec)
 
  702             write(outunit,
'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
 
  703                  'mid time: yyyy/mm/dd hh:mm:ss= ',yr,
'/',mon,
'/',day,hr,
':',minu,
':',sec
 
  704             call get_date(loaded_fields(num_fields)%end_time(j),yr,mon,day,hr,minu,sec)
 
  705             write(outunit,
'(a,i4,a,i2,a,i2,1x,i2,a,i2,a,i2)') &
 
  706                  'end time: yyyy/mm/dd hh:mm:ss= ',yr,
'/',mon,
'/',day,hr,
':',minu,
':',sec
 
  710       deallocate(axisname, axislen)
 
  711       deallocate(tstamp, tstart, tend, tavg)
 
  718     function get_external_fileobj(filename, fileobj)
 
  719        character(len=*),             
intent(in) :: filename
 
  720        type(fmsnetcdffile_t), 
intent(out) :: fileobj
 
  721        logical                                  :: get_external_fileobj
 
  724        get_external_fileobj = .false.
 
  726           if(trim(opened_files(i)%filename) == trim(filename)) 
then 
  727             fileobj = opened_files(i)%fileobj
 
  728             get_external_fileobj = .true.
 
  733     end function get_external_fileobj
 
  736     subroutine set_time_modulo(Time)
 
  738       type(
time_type), 
intent(inout), 
dimension(:) :: time
 
  741       integer :: yr, mon, dy, hr, minu, sec
 
  743       ntime = 
size(time(:))
 
  746          call get_date(time(n), yr, mon, dy, hr, minu, sec)
 
  748          time(n) = 
set_date(yr, mon, dy, hr, minu, sec)
 
  752     end subroutine set_time_modulo
 
  757 subroutine load_record(field, rec, interp, is_in, ie_in, js_in, je_in, window_id_in)
 
  759   integer            ,     
intent(in)           :: rec    
 
  761   integer,                 
intent(in), 
optional :: is_in, ie_in, js_in, je_in
 
  762   integer,                 
intent(in), 
optional :: window_id_in
 
  766   integer :: isw,iew,jsw,jew 
 
  767   integer :: is_region, ie_region, js_region, je_region, i, j
 
  768   integer :: start(4), nread(4)
 
  769   logical :: need_compute
 
  771   real(r8_kind)    :: mask_in(size(field%src_data,1),size(field%src_data,2),size(field%src_data,3))
 
  772   real(r8_kind), 
allocatable :: mask_out(:,:,:)
 
  773   real(r4_kind), 
allocatable :: hi_tmp_data(:,:,:,:)
 
  774   real(r4_kind), 
allocatable :: hi_tmp_msk_out(:,:,:)
 
  777   if( 
PRESENT(window_id_in) ) window_id = window_id_in
 
  778   need_compute = .true.
 
  785      need_compute = .false.
 
  788      field%nbuf = field%nbuf + 1
 
  789      if(field%nbuf > 
size(field%domain_data,4).or.field%nbuf <= 0) field%nbuf = 1
 
  792      field%need_compute(ib,:) = .true.
 
  794      if (debug_this_module) 
write(outunit,*) 
'reading record without domain for field ',trim(field%name)
 
  796      start(1) = field%is_src; nread(1) = field%ie_src - field%is_src + 1
 
  797      start(2) = field%js_src; nread(2) = field%je_src - field%js_src + 1
 
  798      start(3) = 1;            nread(3) = 
size(field%src_data,3)
 
  799      start(field%tdim) = rec; nread(field%tdim) = 1
 
  800      call read_data(field%fileobj,field%name,field%src_data(:,:,:,ib:ib),corner=start,edge_lengths=nread)
 
  803   isw=field%isc;iew=field%iec
 
  804   jsw=field%jsc;jew=field%jec
 
  806   if( field%numwindows > 1) 
then 
  807      if( .NOT. 
PRESENT(is_in) .OR. .NOT. 
PRESENT(ie_in) .OR. .NOT. 
PRESENT(js_in) .OR. .NOT. 
PRESENT(je_in) ) 
then 
  809                   &  
'time_interp_external(load_record): is_in, ie_in, js_in, je_in must be present when numwindows>1')
 
  811      isw = isw + is_in - 1
 
  812      iew = isw + ie_in - is_in
 
  813      jsw = jsw + js_in - 1
 
  814      jew = jsw + je_in - js_in
 
  819   need_compute = field%need_compute(ib, window_id)
 
  820   if(need_compute) 
then 
  821      if(
PRESENT(interp)) 
then 
  822         is_region = field%is_region; ie_region = field%ie_region
 
  823         js_region = field%js_region; je_region = field%je_region
 
  824         mask_in = 0.0_r8_kind
 
  825         where (is_valid(field%src_data(:,:,:,ib), field%valid)) mask_in = 1.0_r8_kind
 
  826         if ( field%region_type .NE. no_region ) 
then 
  827            if( any(mask_in == 0.0_r8_kind) ) 
then 
  828              call mpp_error(fatal, 
"time_interp_external: mask_in should be all 1 when region_type is not NO_REGION")
 
  830            if(interp%interp_method == conserve) 
then 
  832                "conservative interpolation when region_type is not NO_REGION.  Make sure interp matches the region")
 
  834            if( field%region_type == outside_region) 
then 
  835               do j = js_region, je_region
 
  836                  do i = is_region, ie_region
 
  837                     mask_in(i,j,:) = 0.0_r8_kind
 
  841               do j = 1, 
size(mask_in,2)
 
  842                  do i = 1, 
size(mask_in,1)
 
  843                     if( j<js_region .OR. j>je_region .OR. i<is_region .OR. i>ie_region ) mask_in(i,j,:) = 0.0_r8_kind
 
  852         if (interp%horizInterpReals4_type%is_allocated) 
then 
  854             allocate(hi_tmp_msk_out(isw:iew,jsw:jew, 
SIZE(field%src_data,3)))
 
  855             allocate(hi_tmp_data(lbound(field%domain_data,1):ubound(field%domain_data,1), &
 
  856                                  lbound(field%domain_data,2):ubound(field%domain_data,2), &
 
  857                                  lbound(field%domain_data,3):ubound(field%domain_data,3), &
 
  858                                  lbound(field%domain_data,4):ubound(field%domain_data,4)))
 
  860             hi_tmp_data = real(field%domain_data, r4_kind)
 
  862             if(interp%interp_method == conserve) 
then 
  863               call horiz_interp(interp, real(field%src_data(:,:,:,ib),r4_kind), hi_tmp_data(isw:iew,jsw:jew,:,ib))
 
  864               hi_tmp_msk_out = 1.0_r4_kind
 
  866               call horiz_interp(interp, real(field%src_data(:,:,:,ib),r4_kind), hi_tmp_data(isw:iew,jsw:jew,:,ib), &
 
  867                 mask_in=real(mask_in,r4_kind), mask_out=hi_tmp_msk_out)
 
  870             field%domain_data = real(hi_tmp_data, r8_kind)
 
  871             field%mask(isw:iew,jsw:jew,:,ib) = hi_tmp_msk_out(isw:iew,jsw:jew,:) > 0.0_r4_kind
 
  873             if(
allocated(hi_tmp_data))     
deallocate(hi_tmp_data)
 
  874             if(
allocated(hi_tmp_msk_out))  
deallocate(hi_tmp_msk_out)
 
  876             allocate(mask_out(isw:iew,jsw:jew, 
size(field%src_data,3)))
 
  877             if(interp%interp_method == conserve) 
then 
  878               call horiz_interp(interp, field%src_data(:,:,:,ib), field%domain_data(isw:iew,jsw:jew,:,ib))
 
  879               mask_out = 1.0_r8_kind
 
  881               call horiz_interp(interp, field%src_data(:,:,:,ib),field%domain_data(isw:iew,jsw:jew,:,ib), &
 
  882                 mask_in=mask_in, mask_out=mask_out)
 
  884             field%mask(isw:iew,jsw:jew,:,ib) = mask_out(isw:iew,jsw:jew,:) > 0.0_r8_kind
 
  889         if ( field%region_type .NE. no_region ) 
then 
  890            call mpp_error(fatal, 
"time_interp_external: region_type should be NO_REGION when interp is not present")
 
  892         field%domain_data(isw:iew,jsw:jew,:,ib) = field%src_data(isw:iew,jsw:jew,:,ib)
 
  893         field%mask(isw:iew,jsw:jew,:,ib) = is_valid(field%domain_data(isw:iew,jsw:jew,:,ib),field%valid)
 
  896      where(field%mask(isw:iew,jsw:jew,:,ib)) field%domain_data(isw:iew,jsw:jew,:,ib) = &
 
  897           field%domain_data(isw:iew,jsw:jew,:,ib)*field%slope + field%intercept
 
  898      field%need_compute(ib, window_id) = .false.
 
  907   integer            ,     
intent(in)           :: rec    
 
  910   integer :: start(4), nread(4)
 
  918      field%nbuf = field%nbuf + 1
 
  919      if(field%nbuf > 
size(field%domain_data,4).or.field%nbuf <= 0) field%nbuf = 1
 
  923      if (debug_this_module) 
write(outunit,*) 
'reading record without domain for field ',trim(field%name)
 
  925      start(3) = 1;            nread(3) = 
size(field%src_data,3)
 
  926      start(field%tdim) = rec; nread(field%tdim) = 1
 
  927      call read_data(field%fileobj,field%name,field%src_data(:,:,:,ib),corner=start,edge_lengths=nread)
 
  928      if ( field%region_type .NE. no_region ) 
then 
  929         call mpp_error(fatal, 
"time_interp_external: region_type should be NO_REGION when field is scalar")
 
  931      field%domain_data(1,1,:,ib) = field%src_data(1,1,:,ib)
 
  932      field%mask(1,1,:,ib) = is_valid(field%domain_data(1,1,:,ib),field%valid)
 
  934      where(field%mask(1,1,:,ib)) field%domain_data(1,1,:,ib) = &
 
  935           field%domain_data(1,1,:,ib)*field%slope + field%intercept
 
  942    integer, 
intent(in) :: index
 
  943    integer, 
intent(in) :: is, ie, js, je
 
  946    if( is == loaded_fields(index)%is_src .AND. ie == loaded_fields(index)%ie_src .AND. &
 
  947        js == loaded_fields(index)%js_src .AND. ie == loaded_fields(index)%je_src ) 
return 
  949    if( .NOT. 
ASSOCIATED(loaded_fields(index)%src_data) ) 
call mpp_error(fatal, &
 
  950        "time_interp_external: field(index)%src_data is not associated")
 
  951    nk = 
size(loaded_fields(index)%src_data,3)
 
  952    nbuf = 
size(loaded_fields(index)%src_data,4)
 
  953    deallocate(loaded_fields(index)%src_data)
 
  954    allocate(loaded_fields(index)%src_data(is:ie,js:je,nk,nbuf))
 
  955    loaded_fields(index)%is_src = is
 
  956    loaded_fields(index)%ie_src = ie
 
  957    loaded_fields(index)%js_src = js
 
  958    loaded_fields(index)%je_src = je
 
  963 subroutine set_override_region(index, region_type, is_region, ie_region, js_region, je_region)
 
  964    integer, 
intent(in) :: index, region_type
 
  965    integer, 
intent(in) :: is_region, ie_region, js_region, je_region
 
  967    loaded_fields(index)%region_type = region_type
 
  968    loaded_fields(index)%is_region   = is_region
 
  969    loaded_fields(index)%ie_region   = ie_region
 
  970    loaded_fields(index)%js_region   = js_region
 
  971    loaded_fields(index)%je_region   = je_region
 
  975 end subroutine set_override_region
 
  979   integer, 
intent(in) :: n 
 
  984   if (
associated(opened_files)) 
then 
  985      if (n <= 
size(opened_files)) 
return  
  991      if(
Associated(ptr(i)%fileobj)) 
then 
  993         deallocate(ptr(i)%fileobj)
 
  997   if (
associated(opened_files))
then 
  998      ptr(1:
size(opened_files)) = opened_files(:)
 
  999      deallocate(opened_files)
 
 1007   integer, 
intent(in) :: n 
 
 1012   if (
associated(loaded_fields)) 
then 
 1013      if (n <= 
size(loaded_fields)) 
return  
 1018      ptr(i)%fileobj => null()
 
 1023      ptr(i)%domain = null_domain2d
 
 1024      if (
ASSOCIATED(ptr(i)%time))       
DEALLOCATE(ptr(i)%time, stat=ier)
 
 1025      if (
ASSOCIATED(ptr(i)%start_time)) 
DEALLOCATE(ptr(i)%start_time, stat=ier)
 
 1026      if (
ASSOCIATED(ptr(i)%end_time))   
DEALLOCATE(ptr(i)%end_time, stat=ier)
 
 1027      if (
ASSOCIATED(ptr(i)%period)) 
DEALLOCATE(ptr(i)%period, stat=ier)
 
 1028      ptr(i)%modulo_time=.false.
 
 1029      if (
ASSOCIATED(ptr(i)%domain_data)) 
DEALLOCATE(ptr(i)%domain_data, stat=ier)
 
 1030      if (
ASSOCIATED(ptr(i)%ibuf)) 
DEALLOCATE(ptr(i)%ibuf, stat=ier)
 
 1031      if (
ASSOCIATED(ptr(i)%src_data)) 
DEALLOCATE(ptr(i)%src_data, stat=ier)
 
 1033      ptr(i)%domain_present=.false.
 
 1034      ptr(i)%slope=1.0_r8_kind
 
 1035      ptr(i)%intercept=0.0_r8_kind
 
 1036      ptr(i)%isc=-1;ptr(i)%iec=-1
 
 1037      ptr(i)%jsc=-1;ptr(i)%jec=-1
 
 1039   if (
associated(loaded_fields)) 
then 
 1040      ptr(1:
size(loaded_fields)) = loaded_fields(:)
 
 1041      deallocate(loaded_fields)
 
 1051    integer, 
dimension(:) :: buf
 
 1061       if (buf(i) == indx) 
then 
 1077     if (index .lt. 1 .or. index .gt. num_fields) &
 
 1078         call mpp_error(fatal,
'invalid index in call to get_external_field_size')
 
 1094     if (index .lt. 1 .or. index .gt. num_fields) &
 
 1095         call mpp_error(fatal,
'invalid index in call to get_external_field_size')
 
 1103   integer        , 
intent(in)  :: index
 
 1108   if (index < 1.or.index > num_fields) &
 
 1109        call mpp_error(fatal,
'invalid index in call to get_time_axis')
 
 1111   n = min(
size(time),
size(loaded_fields(index)%time))
 
 1113   time(1:n) = loaded_fields(index)%time(1:n)
 
 1126         deallocate(loaded_fields(i)%time,loaded_fields(i)%start_time,loaded_fields(i)%end_time,&
 
 1127                    loaded_fields(i)%period,loaded_fields(i)%domain_data,loaded_fields(i)%mask,loaded_fields(i)%ibuf)
 
 1128         if (
ASSOCIATED(loaded_fields(i)%src_data)) 
deallocate(loaded_fields(i)%src_data)
 
 1129         loaded_fields(i)%domain = null_domain2d
 
 1130         loaded_fields(i)%nbuf = 0
 
 1131         loaded_fields(i)%slope = 0.0_r8_kind
 
 1132         loaded_fields(i)%intercept = 0.0_r8_kind
 
 1138         deallocate(opened_files(i)%fileobj)
 
 1141     deallocate(loaded_fields)
 
 1142     deallocate(opened_files)
 
 1146     module_initialized = .false.
 
 1150 #include "time_interp_external2_r4.fh" 
 1151 #include "time_interp_external2_r8.fh" 
 1153 #include "time_interp_external2_bridge_r4.fh" 
 1154 #include "time_interp_external2_bridge_r8.fh" 
 1156 end module time_interp_external2_mod
 
subroutine, public get_axis_cart(fileobj, axisname, cart)
Returns X,Y,Z or T cartesian attribute.
logical function, public get_axis_modulo(fileobj, axisname)
Checks if 'modulo' variable exists for a given axis.
logical function, public get_axis_modulo_times(fileobj, axisname, tbeg, tend)
Close a netcdf or domain file opened with open_file or open_virtual_file.
Opens a given netcdf or domain file.
Read data from a defined field in a file.
integer function, public check_nml_error(IOSTAT, NML_NAME)
Checks the iostat argument that is returned after reading a namelist and determines if the error code...
subroutine, public write_version_number(version, tag, unit)
Prints to the log file (or a specified unit) the version id string and tag name.
Added for mixed precision support. Updates force time_manager math to be done with kind=8 reals _wrap...
Subroutine for performing the horizontal interpolation between two grids.
Holds data pointers and metadata for horizontal interpolations, passed between the horiz_interp modul...
These routines retrieve the axis specifications associated with the compute domains....
These routines retrieve the axis specifications associated with the data domains. The domain is a der...
These routines retrieve the axis specifications associated with the global domains....
The domain2D type contains all the necessary information to define the global, compute and data domai...
integer function stdout()
This function returns the current standard fortran unit numbers for output.
integer function stdlog()
This function returns the current standard fortran unit numbers for log messages. Log messages,...
integer function mpp_npes()
Returns processor count for current pelist.
integer function mpp_pe()
Returns processor ID.
Perform parallel broadcasts.
subroutine, public time_interp_external_exit()
exit time_interp_external_mod. Close all open files and release storage
integer function, dimension(4), public get_external_field_size(index)
Returns size of field after call to init_external_field. Ordering is X/Y/Z/T. This call only makes se...
subroutine, public reset_src_data_region(index, is, ie, js, je)
Reallocates src_data for field from module level loaded_fields array.
subroutine load_record(field, rec, interp, is_in, ie_in, js_in, je_in, window_id_in)
load specified record from file Always loads in r8, casts down for horiz_interp if interp argument is...
subroutine, public time_interp_external_init()
Initialize the time_interp_external_mod module.
real(r8_kind) function, public get_external_field_missing(index)
return missing value
subroutine realloc_files(n)
reallocates array of fields, increasing its size
integer function, public init_external_field(file, fieldname, domain, desired_units, verbose, axis_names, axis_sizes, override, correct_leap_year_inconsistency, permit_calendar_conversion, use_comp_domain, ierr, nwindows, ignore_axis_atts, ongrid)
Initialize an external field. Buffer "num_io_buffers" (default=2) in memory to reduce memory allocati...
integer function, private find_buf_index(indx, buf)
simple linear search for given value in given list TODO should use better search if this list is bigg...
subroutine, public get_time_axis(index, time)
subroutine realloc_fields(n)
reallocates array of fields,increasing its size
subroutine load_record_0d(field, rec)
Given a initialized ext_fieldtype and record number, loads the given index into fieldsrc_data.
Provide data from external file interpolated to current model time. Data may be local to current proc...
Represents external fields.
Holds filename and file object.
Returns a weight and dates or indices for interpolating between two dates. The interface fraction_of_...
subroutine, public get_time(Time, seconds, days, ticks, err_msg)
Returns days and seconds ( < 86400 ) corresponding to a time. err_msg should be checked for any error...
subroutine, public get_date(time, year, month, day, hour, minute, second, tick, err_msg)
Gets the date for different calendar types. Given a time_interval, returns the corresponding date und...
integer function, public days_in_year(Time)
Returns the number of days in the calendar year corresponding to the date represented by time for the...
type(time_type) function, public increment_time(Time, seconds, days, ticks, err_msg, allow_neg_inc)
Increments a time by seconds and days.
integer function, public days_in_month(Time, err_msg)
Given a time, computes the corresponding date given the selected date time mapping algorithm.
integer function, public get_calendar_type()
Returns default calendar type for mapping from time to date.
Given an input date in year, month, days, etc., creates a time_type that represents this time interva...
Given some number of seconds and days, returns the corresponding time_type.
Type to represent amounts of time. Implemented as seconds and days to allow for larger intervals.