diff --git a/horiz_interp/horiz_interp_bilinear.F90 b/horiz_interp/horiz_interp_bilinear.F90
index 0365bc6b85..b745461a53 100644
--- a/horiz_interp/horiz_interp_bilinear.F90
+++ b/horiz_interp/horiz_interp_bilinear.F90
@@ -555,7 +555,7 @@ subroutine find_neighbor( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo )
jlat(1) = js
else
npts = 0
- !--- bottom and top boundary
+ !--- bottom boundary
jstart = max(js-step,1)
jend = min(js+step,nlat_in)
@@ -576,9 +576,6 @@ subroutine find_neighbor( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo )
npts = npts + 1
ilon(npts) = i
jlat(npts) = jstart
- npts = npts + 1
- ilon(npts) = i
- jlat(npts) = jend
enddo
!--- right and left boundary -----------------------------------------------
@@ -593,7 +590,7 @@ subroutine find_neighbor( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo )
endif
do l = -step, step
j = js+l
- if( j < 1 .or. j > nlat_in) cycle
+ if( j < 1 .or. j > nlat_in .or. j==jstart .or. j==jend) cycle
npts = npts+1
ilon(npts) = istart
jlat(npts) = j
@@ -601,6 +598,29 @@ subroutine find_neighbor( Interp, lon_in, lat_in, lon_out, lat_out, src_modulo )
ilon(npts) = iend
jlat(npts) = j
end do
+
+ !--- top boundary
+
+ do l = -step, step
+ i = is+l
+ if(src_modulo)then
+ if( i < 1) then
+ i = i + nlon_in
+ else if (i > nlon_in) then
+ i = i - nlon_in
+ endif
+ if( i < 1 .or. i > nlon_in) call mpp_error(FATAL, &
+ 'horiz_interp_bilinear_mod: max_step is too big, decrease max_step' )
+ else
+ if( i < 1 .or. i > nlon_in) cycle
+ endif
+
+ npts = npts + 1
+ ilon(npts) = i
+ jlat(npts) = jend
+ enddo
+
+
end if
!--- find the surrouding points
diff --git a/horiz_interp/test_horiz_interp.F90 b/horiz_interp/test_horiz_interp.F90
index 9efa898ed4..cc88dbe480 100644
--- a/horiz_interp/test_horiz_interp.F90
+++ b/horiz_interp/test_horiz_interp.F90
@@ -4,7 +4,7 @@
program test
use mpp_mod, only : mpp_error, mpp_pe, mpp_npes, mpp_root_pe
- use mpp_mod, only : FATAL, stdout, stdlog
+ use mpp_mod, only : FATAL, stdout, stdlog, mpp_chksum
use mpp_io_mod, only : mpp_open, mpp_close, mpp_read
use mpp_io_mod, only : axistype, fieldtype
use mpp_io_mod, only : mpp_get_info, mpp_get_fields, mpp_get_times
@@ -33,8 +33,12 @@ program test
character(len=256) :: dst_file = "output.nc"
logical :: new_missing_handle = .false.
character(len=256) :: interp_method = "bilinear"
+ logical :: use_2d_version = .false.
+ logical :: write_remap_index = .false.
+ integer :: layout(2) = (/1,1/)
real, allocatable, dimension(:) :: x_src, y_src
+ real, allocatable, dimension(:,:) :: x_src_2d, y_src_2d
real, allocatable, dimension(:,:) :: x_dst, y_dst
type(axistype), allocatable, dimension(:) :: axes
type(fieldtype), allocatable, dimension(:) :: fields
@@ -45,7 +49,7 @@ program test
integer :: n, ntimes
namelist /test_horiz_interp_nml/ src_file, field_name, dst_file, dst_grid, new_missing_handle, &
- interp_method
+ interp_method, use_2d_version, layout, write_remap_index
call fms_init
call horiz_interp_init
@@ -100,18 +104,26 @@ subroutine process_data()
type(horiz_interp_type) :: Interp
type(axistype) :: xaxis, yaxis, zaxis, taxis
type(domain1D) :: xdom, ydom
- type(fieldtype) :: field
+ type(fieldtype) :: field, field_istart,field_iend,field_jstart,field_jend
real :: D2R
integer :: k, i
call mpp_get_domain_components( domain, xdom, ydom )
!--- set up meta data
+! call mpp_open(unit,dst_file,action=MPP_OVERWR,form=MPP_NETCDF,threading=MPP_MULTI, fileset=MPP_MULTI)
call mpp_open(unit,dst_file,action=MPP_OVERWR,form=MPP_NETCDF,threading=MPP_MULTI, fileset=MPP_MULTI)
call mpp_write_meta( unit, xaxis, 'lon', 'none', 'X index', 'X', domain=xdom, data=(/(i-1.,i=1,nx_dst)/) )
call mpp_write_meta( unit, yaxis, 'lat', 'none', 'Y index', 'Y', domain=ydom, data=(/(i-1.,i=1,ny_dst)/) )
call mpp_write_meta( unit, zaxis, 'level', 'none', 'Z index', 'Z', data=(/(i-1.,i=1,nk)/) )
call mpp_write_meta( unit, taxis, 'time', 'none', 'T index', 'T' )
- call mpp_write_meta( unit, field, (/xaxis, yaxis, zaxis, taxis/), field_name, 'none', 'none')
+ call mpp_write_meta( unit, field, (/xaxis, yaxis, zaxis, taxis/), field_name, 'none', 'none', missing=missing_value)
+ if(write_remap_index) then
+ call mpp_write_meta( unit, field_istart, (/xaxis, yaxis/), "istart", 'none', 'none')
+ call mpp_write_meta( unit, field_iend, (/xaxis, yaxis/), "iend", 'none', 'none')
+ call mpp_write_meta( unit, field_jstart, (/xaxis, yaxis/), "jstart", 'none', 'none')
+ call mpp_write_meta( unit, field_jend, (/xaxis, yaxis/), "jend", 'none', 'none')
+ endif
+
call mpp_write( unit, xaxis )
call mpp_write( unit, yaxis )
call mpp_write( unit, zaxis )
@@ -122,9 +134,28 @@ subroutine process_data()
allocate(src_data(nx_src,ny_src,nk))
allocate(dst_data(is:ie,js:je,nk))
- call horiz_interp_new(Interp, x_src*D2R, y_src*D2R, x_dst*D2R, y_dst*D2R, &
+ if(trim(interp_method) == "bilinear" .and. use_2d_version) then
+ write(stdout(),*) "use 2-D version of bilinear interpolation"
+ call horiz_interp_new(Interp, x_src_2d*D2R, y_src_2d*D2R, x_dst*D2R, y_dst*D2R, &
+ interp_method = trim(interp_method) )
+ else
+ write(stdout(),*) "use 1-D version of interpolation"
+ call horiz_interp_new(Interp, x_src*D2R, y_src*D2R, x_dst*D2R, y_dst*D2R, &
interp_method = trim(interp_method), grid_at_center = .true. )
-
+ endif
+
+ if(write_remap_index) then
+ dst_data(:,:,1) = interp%i_lon(:,:,1)
+ call mpp_write(unit, field_istart, domain, dst_data(:,:,1))
+ dst_data(:,:,1) = interp%i_lon(:,:,2)
+ call mpp_write(unit, field_iend, domain, dst_data(:,:,1))
+ dst_data(:,:,1) = interp%j_lat(:,:,1)
+ call mpp_write(unit, field_jstart, domain, dst_data(:,:,1))
+ dst_data(:,:,1) = interp%j_lat(:,:,2)
+ call mpp_write(unit, field_jend, domain, dst_data(:,:,1))
+ endif
+
+ dst_data = 0
do n = 1, ntimes
call mpp_read(src_unit,fields(src_field_index),src_data, tindex=n)
@@ -133,6 +164,7 @@ subroutine process_data()
missing_value=missing_value, new_missing_handle=new_missing_handle)
enddo
call mpp_write(unit, field, domain, dst_data, n*1.0)
+ write(stdout(),*) "chksum at time = ", n, " = ", mpp_chksum(dst_data)
enddo
call mpp_close(unit)
@@ -142,7 +174,7 @@ end subroutine process_data
subroutine read_dst_grid
- integer :: layout(2), start(4), nread(4)
+ integer :: start(4), nread(4)
character(len=256) :: tile_file
integer :: i, j, siz(4)
real, allocatable :: tmp(:,:)
@@ -231,11 +263,21 @@ subroutine read_src_file
nk = len1
end select
enddo
+
if(nx_src==0) call mpp_error(FATAL,'test_horiz_interp: file ' &
//trim(src_file)//' does not contain axis with cartesian attributes = "X" ')
if(ny_src==0) call mpp_error(FATAL,'test_horiz_interp: file '&
//trim(src_file)//' does not contain axis with cartesian attributes = "Y" ')
+ allocate(x_src_2d(nx_src,ny_src), y_src_2d(nx_src,ny_src))
+ do i = 1, nx_src
+ x_src_2d(i,:) = x_src(i)
+ enddo
+ do i = 1, ny_src
+ y_src_2d(:,i) = y_src(i)
+ enddo
+
+
!--- get the missing value
call mpp_get_atts(fields(src_field_index),missing=missing_value)
diff --git a/mpp/include/mpp_io_connect.inc b/mpp/include/mpp_io_connect.inc
index ab1ab0940b..008e46c039 100644
--- a/mpp/include/mpp_io_connect.inc
+++ b/mpp/include/mpp_io_connect.inc
@@ -513,7 +513,7 @@
if( verbose )print '(a,i6,i16,i4)', 'MPP_OPEN: append to existing netCDF file: pe, ncid, time_axis_id=',&
pe, mpp_file(unit)%ncid, mpp_file(unit)%id
mpp_file(unit)%format=form_flag ! need this for mpp_read
-! call mpp_read_meta(unit)
+ call mpp_read_meta(unit, read_time=.FALSE.)
else if( action_flag.EQ.MPP_RDONLY )then
inquire(file=trim(mpp_file(unit)%name),EXIST=exists)
if (.NOT.exists) call mpp_error(FATAL,'MPP_OPEN:'&
@@ -522,7 +522,7 @@
if( verbose )print '(a,i6,i16,i4)', 'MPP_OPEN: opening existing netCDF file: pe, ncid, time_axis_id=',&
pe, mpp_file(unit)%ncid, mpp_file(unit)%id
mpp_file(unit)%format=form_flag ! need this for mpp_read
- call mpp_read_meta(unit)
+ call mpp_read_meta(unit, read_time=.TRUE.)
end if
mpp_file(unit)%opened = .TRUE.
#elif use_LARGEFILE
@@ -553,7 +553,7 @@
if( verbose )print '(a,i6,i16,i4)', 'MPP_OPEN: append to existing netCDF file: pe, ncid, time_axis_id=',&
pe, mpp_file(unit)%ncid, mpp_file(unit)%id
mpp_file(unit)%format=form_flag ! need this for mpp_read
-! call mpp_read_meta(unit)
+ call mpp_read_meta(unit, read_time=.FALSE.)
else if( action_flag.EQ.MPP_RDONLY )then
inquire(file=trim(mpp_file(unit)%name),EXIST=exists)
if (.NOT.exists) call mpp_error(FATAL,'MPP_OPEN:'&
@@ -562,7 +562,7 @@
if( verbose )print '(a,i6,i16,i4)', 'MPP_OPEN: opening existing netCDF file: pe, ncid, time_axis_id=',&
pe, mpp_file(unit)%ncid, mpp_file(unit)%id
mpp_file(unit)%format=form_flag ! need this for mpp_read
- call mpp_read_meta(unit)
+ call mpp_read_meta(unit, read_time=.TRUE.)
end if
mpp_file(unit)%opened = .TRUE.
#else
@@ -593,7 +593,7 @@
if( verbose )print '(a,i6,i16,i4)', 'MPP_OPEN: append to existing netCDF file: pe, ncid, time_axis_id=',&
pe, mpp_file(unit)%ncid, mpp_file(unit)%id
mpp_file(unit)%format=form_flag ! need this for mpp_read
-! call mpp_read_meta(unit)
+ call mpp_read_meta(unit, read_time=.FALSE.)
else if( action_flag.EQ.MPP_RDONLY )then
inquire(file=trim(mpp_file(unit)%name),EXIST=exists)
if (.NOT.exists) call mpp_error(FATAL,'MPP_OPEN:'&
@@ -602,7 +602,7 @@
if( verbose )print '(a,i6,i16,i4)', 'MPP_OPEN: opening existing netCDF file: pe, ncid, time_axis_id=',&
pe, mpp_file(unit)%ncid, mpp_file(unit)%id
mpp_file(unit)%format=form_flag ! need this for mpp_read
- call mpp_read_meta(unit)
+ call mpp_read_meta(unit, read_time=.TRUE.)
end if
mpp_file(unit)%opened = .TRUE.
diff --git a/mpp/include/mpp_io_read.inc b/mpp/include/mpp_io_read.inc
index 4ef6181db6..09179a602b 100644
--- a/mpp/include/mpp_io_read.inc
+++ b/mpp/include/mpp_io_read.inc
@@ -210,7 +210,7 @@
start = 1
do i = 1,size(field%axes(:))
axsiz(i) = field%size(i)
- if( field%axes(i)%did.EQ.field%time_axis_index .AND. field%time_axis_index >0 )start(i) = tlevel
+ if( i .EQ. field%time_axis_index )start(i) = tlevel
end do
if( PRESENT(domain) )then
@@ -448,7 +448,7 @@
! mpp_read_meta must be called prior to mpp_read.
!
!
- subroutine mpp_read_meta(unit)
+ subroutine mpp_read_meta(unit, read_time)
!
! read file attributes including dimension and variable attributes
! and store in filetype structure. All of the file information
@@ -458,13 +458,14 @@
! every PE is eligible to call mpp_read_meta
!
integer, intent(in) :: unit
-
+ logical, intent(in), optional :: read_time ! read_time is set to false when file action is appending or writing.
+ ! This is temporary fix for MOM6 reopen_file issue.
integer :: ncid,ndim,nvar_total,natt,recdim,nv,nvar,len
integer :: error, i, j, istat, check_exist
integer :: type, nvdims, nvatts, dimid
integer, allocatable, dimension(:) :: dimids
character(len=128) :: name, attname, unlimname, attval, bounds_name
- logical :: isdim, found_bounds
+ logical :: isdim, found_bounds, get_time_info
integer(LONG_KIND) :: checksumf
character(len=64) :: checksum_char
integer :: num_checksumf, last, is, k
@@ -474,6 +475,9 @@
real(FLOAT_KIND), allocatable :: rvals(:)
real(DOUBLE_KIND), allocatable :: r8vals(:)
+ get_time_info = .TRUE.
+ if(present(read_time)) get_time_info = read_time
+
#ifdef use_netCDF
if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
@@ -717,7 +721,7 @@
case default
call mpp_error( FATAL, 'Invalid data type for dimension' )
end select
- else
+ else if(get_time_info) then
len = mpp_file(unit)%time_level
if( len > 0 ) then
allocate(mpp_file(unit)%time_values(len), STAT=istat)
@@ -1060,7 +1064,8 @@
do j=1,nvdims
if(dimids(j).eq.mpp_file(unit)%recdimid .and. mpp_file(unit)%time_level/=-1)then
- mpp_file(unit)%Var(nv)%time_axis_index = dimids(j)
+ mpp_file(unit)%Var(nv)%time_axis_index = j !dimids(j). z1l: Should be j.
+ !This will cause problem when appending to existed file.
mpp_file(unit)%Var(nv)%size(j)=1 ! dimid length set to 1 here for consistency w/ mpp_write
else
mpp_file(unit)%Var(nv)%size(j)=mpp_file(unit)%Axis(dimids(j))%len
diff --git a/mpp/test_mpp_io.F90 b/mpp/test_mpp_io.F90
index 7774f2e97e..8ec4024d69 100644
--- a/mpp/test_mpp_io.F90
+++ b/mpp/test_mpp_io.F90
@@ -14,7 +14,7 @@ program test
use mpp_io_mod, only : MPP_RDONLY, mpp_open, MPP_OVERWR, MPP_ASCII, MPP_SINGLE
use mpp_io_mod, only : MPP_NETCDF, MPP_MULTI, mpp_get_atts, mpp_write, mpp_close
use mpp_io_mod, only : mpp_get_info, mpp_get_axes, mpp_get_fields, mpp_get_times
- use mpp_io_mod, only : mpp_read, mpp_io_exit
+ use mpp_io_mod, only : mpp_read, mpp_io_exit, MPP_APPEND
#ifdef INTERNAL_FILE_NML
USE mpp_mod, ONLY: input_nml_file
@@ -103,6 +103,8 @@ program test
pack_size = size(transfer(doubledata, realarray))
if( pack_size .NE. 1 .AND. pack_size .NE. 2) call mpp_error(FATAL,'test_mpp_io: pack_size should be 1 or 2')
+ call test_netcdf_io_append()
+
if(ntiles_x == 1 .and. ntiles_y == 1 .and. io_layout(1) == 1 .and. io_layout(2) == 1) then
call test_netcdf_io('Simple')
call test_netcdf_io('Symmetry_T_cell')
@@ -360,6 +362,52 @@ subroutine test_netcdf_io(type)
end subroutine test_netcdf_io
+ subroutine test_netcdf_io_append
+ integer :: ndim, nvar, natt, ntime
+ integer :: is, ie, js, je, isd, ied, jsd, jed, ism, iem, jsm, jem
+ integer :: position, msize(2), ioff, joff, nxg, nyg
+ logical :: symmetry
+ real :: data
+ type(fieldtype), allocatable :: vars(:)
+
+!netCDF distributed write
+ if( pe.EQ.mpp_root_pe() )print *, 'netCDF single thread write'
+ call mpp_open( unit, "timestats.nc", action=MPP_OVERWR, &
+ form=MPP_NETCDF, threading=MPP_SINGLE, fileset=MPP_SINGLE )
+ call mpp_write_meta( unit, t, 'T', 'sec', 'Time', 'T' )
+ call mpp_write_meta( unit, f, (/t/), 'Data', 'metres', 'Random data', pack=pack_size )
+ do i = 0,nt-1
+ time = i*1.
+ data = i*3.0
+ call mpp_write( unit, f, data, time )
+ end do
+ call mpp_close(unit)
+
+!--- append
+ if( pe.EQ.mpp_root_pe() )print *, 'netCDF single thread append'
+ call mpp_open( unit, "timestats.nc", action=MPP_APPEND, &
+ form=MPP_NETCDF, threading=MPP_SINGLE, fileset=MPP_SINGLE )
+ allocate(vars(1))
+ if(pe.EQ.mpp_root_pe() ) then
+ call mpp_get_info(unit, ndim, nvar, natt, ntime)
+
+ if (nvar /= 1) then
+ call mpp_error(FATAL, "test_netcdf_io_append: nvar should be 1")
+ endif
+ call mpp_get_fields(unit,vars(1:nvar))
+ endif
+
+ do i = nt,2*nt-1
+ time = i*1.
+ data = i*3.0
+ call mpp_write( unit, vars(1), data, time )
+ end do
+ call mpp_close(unit)
+ deallocate(vars)
+
+
+ end subroutine test_netcdf_io_append
+
subroutine test_netcdf_io_mosaic(type, layout, ntiles_x, ntiles_y, io_layout)
character(len=*), intent(in) :: type