Merge 'user/z1l/horiz_interp' into dev/master

* Fix a reproducing ability in the 2-D bilinear interpolation. This
  might change answer for some MOM6 test.
* Fix a crash issue for a CM4 experiment

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,14 +590,37 @@ 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
                          npts = npts+1
                          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             

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)
 

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.
 

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 @@
 !     <TT>mpp_read_meta</TT> must be called prior to <TT>mpp_read</TT>.
 !   </NOTE>
 ! </SUBROUTINE>
-    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

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