#ifdef TEST_HORIZ_INTERP !z1l: currently only test bilinear interpolation. program test use mpp_mod, only : mpp_error, mpp_pe, mpp_npes, mpp_root_pe use mpp_mod, only : FATAL, stdout, stdlog 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 use mpp_io_mod, only : mpp_get_axes, mpp_get_axis_data, mpp_get_atts use mpp_io_mod, only : MPP_RDONLY, MPP_NETCDF, MPP_MULTI, MPP_SINGLE, MPP_OVERWR use mpp_io_mod, only : mpp_get_att_name, mpp_get_att_char, mpp_get_att_type, mpp_get_att_real use mpp_io_mod, only : mpp_write_meta, axistype, fieldtype, mpp_write, mpp_close use mpp_domains_mod, only : mpp_update_domains, mpp_define_domains, domain1d use mpp_domains_mod, only : domain2d, mpp_define_layout, mpp_get_compute_domain use mpp_domains_mod, only : mpp_get_domain_components use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_end, horiz_interp_type use horiz_interp_mod, only : horiz_interp_init use axis_utils_mod, only : get_axis_cart use fms_io_mod, only : read_data, write_data use fms_io_mod, only : field_size, fms_io_exit use fms_mod, only : fms_init, fms_end, open_namelist_file, close_file, file_exist, field_exist use fms_mod, only : check_nml_error, write_version_number, lowercase use constants_mod, only : constants_init, PI use horiz_interp_mod, only : horiz_interp_new, horiz_interp, horiz_interp_end, horiz_interp_type implicit none character(len=256) :: src_file = "" character(len=256) :: dst_grid = "INPUT/grid_spec.nc" character(len=256) :: field_name = "" character(len=256) :: dst_file = "output.nc" logical :: new_missing_handle = .false. character(len=256) :: interp_method = "bilinear" real, allocatable, dimension(:) :: x_src, y_src real, allocatable, dimension(:,:) :: x_dst, y_dst type(axistype), allocatable, dimension(:) :: axes type(fieldtype), allocatable, dimension(:) :: fields type(domain2d) :: Domain integer :: unit, ierr, io, src_unit, src_field_index, nk integer :: nx_src, ny_src, nx_dst, ny_dst, is, ie, js, je real :: missing_value integer :: n, ntimes namelist /test_horiz_interp_nml/ src_file, field_name, dst_file, dst_grid, new_missing_handle, & interp_method call fms_init call horiz_interp_init call constants_init #ifdef INTERNAL_FILE_NML read (input_nml_file, test_horiz_interp_nml, iostat=io) ierr = check_nml_error(io, 'test_horiz_interp_nml') #else if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr=1 do while (ierr /= 0) read(unit, nml=test_horiz_interp_nml, iostat=io, end=10) ierr = check_nml_error(io, 'test_horiz_interp_nml') enddo 10 call close_file (unit) endif #endif if( .not. file_exist(src_file) ) call mpp_error(FATAL, & "test_horiz_interp: src_file = "//trim(src_file)//" does not exist") if( .not. field_exist(src_file, field_name) ) call mpp_error(FATAL, & "test_horiz_interp: field_name = "//trim(field_name)//" does not exist in file "//trim(src_file) ) ! reading the grid information and missing value from src_file call mpp_open(src_unit, trim(src_file), & action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_MULTI, fileset=MPP_SINGLE) call read_src_file() ! reading the dst_grid file call read_dst_grid() !--- currently only test for first time level. The following will read the input data, !--- do the remapping and write out data call process_data() call mpp_close(src_unit) call fms_io_exit call fms_end() contains ! read the input data, do the remapping and write out data subroutine process_data() real, allocatable, dimension(:,:,:) :: src_data, dst_data type(horiz_interp_type) :: Interp type(axistype) :: xaxis, yaxis, zaxis, taxis type(domain1D) :: xdom, ydom type(fieldtype) :: field 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_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( unit, xaxis ) call mpp_write( unit, yaxis ) call mpp_write( unit, zaxis ) D2R = PI/180. 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, & interp_method = trim(interp_method), grid_at_center = .true. ) do n = 1, ntimes call mpp_read(src_unit,fields(src_field_index),src_data, tindex=n) do k = 1, nk call horiz_interp(interp, src_data(:,:,k), dst_data(:,:,k), & missing_value=missing_value, new_missing_handle=new_missing_handle) enddo call mpp_write(unit, field, domain, dst_data, n*1.0) enddo call mpp_close(unit) deallocate(src_data, dst_data) end subroutine process_data subroutine read_dst_grid integer :: layout(2), start(4), nread(4) character(len=256) :: tile_file integer :: i, j, siz(4) real, allocatable :: tmp(:,:) if (field_exist(dst_grid, "gridfiles" )) then call read_data(dst_grid, "gridfiles", tile_file) tile_file = 'INPUT/'//trim(tile_file) else call mpp_error(FATAL, "test_horiz_interp: field gridfiles does not exist in file "//trim(dst_grid) ) endif call field_size(tile_file, "x", siz) nx_dst = (siz(1)-1)/2 ny_dst = (siz(2)-1)/2 if(layout(1)*layout(2) .NE. mpp_npes() ) then call mpp_define_layout( (/1,nx_dst,1,ny_dst/), mpp_npes(), layout ) end if call mpp_define_domains( (/1,nx_dst,1,ny_dst/), layout, Domain, name='test_data_override') call mpp_get_compute_domain(Domain, is, ie, js, je) allocate(tmp(2*is-1:2*ie+1,2*js-1:2*je+1)) allocate(x_dst(is:ie,js:je), y_dst(is:ie,js:je)) start = 1; nread = 1 start(1) = 2*is-1; nread(1) = 2*(ie-is+1)+1 start(2) = 2*js-1; nread(2) = 2*(je-js+1)+1 call read_data(tile_file, "x", tmp, start, nread, domain) do j = js, je do i = is, ie x_dst(i,j) = tmp(2*i,2*j) enddo enddo call read_data(tile_file, "y", tmp, start, nread, domain) do j = js, je do i = is, ie y_dst(i,j) = tmp(2*i,2*j) enddo enddo deallocate(tmp) end subroutine read_dst_grid subroutine read_src_file integer :: ndim, nvar, natt, n integer :: nt, i, j, k, jj, len1 character(len=1) :: cart character(len=32) :: name call mpp_get_info(src_unit, ndim, nvar, natt, ntimes) allocate(fields(nvar)) call mpp_get_fields(src_unit, fields) src_field_index = 0 do i=1,nvar call mpp_get_atts(fields(i),name=name) if (lowercase(trim(field_name)) == lowercase(trim(name))) then src_field_index = i endif enddo if(src_field_index == 0) call mpp_error(FATAL, 'test_horiz_interp: field '& //trim(field_name)//' is not in the file '//trim(src_file) ) !--- get the src grid call mpp_get_atts(fields(src_field_index),ndim=ndim) allocate(axes(ndim)) call mpp_get_atts(fields(src_field_index),axes=axes) nx_src=0; ny_src=0; nk=1 do j=1,ndim call mpp_get_atts(axes(j),len=len1) call get_axis_cart(axes(j),cart) select case (cart) case ('X') nx_src = len1 allocate(x_src(nx_src)) call mpp_get_axis_data(axes(j),x_src) case('Y') ny_src = len1 allocate(y_src(ny_src)) call mpp_get_axis_data(axes(j),y_src) case('Z') 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" ') !--- get the missing value call mpp_get_atts(fields(src_field_index),missing=missing_value) end subroutine read_src_file end program test #else module null_test_horiz_interp end module #endif