!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! 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 General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program compare_grid
!-----------------------------------------------------------------------
! GNU General Public License
!
! This program is free software; you can redistribute it and/or modify it and
! are expected to follow the terms of the GNU General Public License
! as published by the Free Software Foundation; either version 2 of
! the License, or (at your option) any later version.
!
! MOM 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.
! For the full text of the GNU General Public License,
! write to: Free Software Foundation, Inc.,
! 675 Mass Ave, Cambridge, MA 02139, USA.
! or see: http://www.gnu.org/licenses/gpl.html
!-----------------------------------------------------------------------
! Z. Liang
! S. M. Griffies
!
! compare depth and land/sea mask of two grid_spec file.
!
!
! This program reads in two grid descriptor files (generated via ocean_grid_generator)
! and creates a text file output listing line-by-line differences between the
! two files. Output file format is the same as the grid_edits file used by
! edit_grid.F90. These two files should have same grid size.
!
! Originally developed by Jeffery B. Greenblatt on 12/11/2001 at Princeton University
!
use mpp_mod, only : mpp_error, FATAL, NOTE, mpp_pe, mpp_init, mpp_exit
use mpp_io_mod, only : mpp_open, mpp_read, mpp_close, fieldtype
use mpp_io_mod, only : mpp_get_info, mpp_get_atts, mpp_get_fields
use mpp_io_mod, only : MPP_RDONLY, MPP_NETCDF, MPP_SINGLE, MPP_OVERWR, MPP_ASCII
use fms_mod, only : fms_init, fms_end, file_exist, stdout, close_file
use fms_mod, only : open_namelist_file, check_nml_error, write_version_number
implicit none
!--- namelist information --------------------------------------------
!
!
!
! First grid files to be compared with grid_file_2.
!
!
! Second grid files to be compared with grid_file_1.
!
!
! output text file. Each line is in the format as
! "i, j, depth_new, #was depth_old ". depth_new is
! the depth at point (i,j) of grid_file_2 and depth_old
! is the depth at point (i,j) of grid_file_1.
!
!
! output text file. Each line is in the format as
! "i, j, wet_new, #was wet_old ". wet_new is
! the land/sea mask at point (i,j) of grid_file_2 and wet_old
! is the land/sea mask at point (i,j) of grid_file_1.
!
!
character(len=128) :: grid_file_1 = 'grid_file_1'
character(len=128) :: grid_file_2 = 'grid_file_2'
character(len=128) :: grid_edits = 'grid_edits.txt'
character(len=128) :: mask_diff = 'mask_diff.txt'
namelist /compare_grid_nml/ grid_file_1, grid_file_2, grid_edits, mask_diff
!--- version information variables -----------------------------------
character(len=128) :: version = '$Id: compare_grid.f90,v 11.0 2004/09/28 20:07:16 fms Exp $'
character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $'
!--- compare_grid_type
type compare_grid_type
real, dimension(:,:), pointer :: ht => NULL()
integer, dimension(:,:), pointer :: wet => NULL()
real, dimension(:,:), pointer :: geolon_t => NULL()
real, dimension(:,:), pointer :: geolat_t => NULL()
integer :: ni, nj
end type compare_grid_type
!--- other variables
type(compare_grid_type) :: grid_1
type(compare_grid_type) :: grid_2
call fms_init ()
call compare_grid_init ()
! --- Read grid files ------------------------------------------------
call read_grid(grid_file_1, grid_1)
call read_grid(grid_file_2, grid_2)
! Ensure we are comparing same-size arrays
call grid_check()
! Now compare selected arrays and write out differences
call grid_compare()
call fms_end()
contains
!#####################################################################
! --- read the namelist and write the version information to logfile. Also
! --- write the namelist to standard output
subroutine compare_grid_init
integer :: io_status, unit, ierr
! --- provide for namelist over-ride of defaults ---------------------
if(file_exist('input.nml')) then
unit = open_namelist_file()
read (unit,compare_grid_nml,IOSTAT=io_status)
write (stdout(),'(/)')
write (stdout(),compare_grid_nml)
ierr = check_nml_error(io_status, 'compare_grid_nml')
call close_file(unit)
else
call mpp_error(NOTE, 'file input.nml does not exist' )
endif
!--- write out version information ---------------------------------
call write_version_number(version,tagname)
end subroutine compare_grid_init
!#####################################################################
!--- read the grid information from file.
subroutine read_grid(file, grid)
character(len=*), intent(in) :: file
type(compare_grid_type), intent(inout) :: grid
!--- Local variables
integer :: i, j, unit, len, ndim, nvar, natt, ntime
integer :: siz_in(3)
logical :: ht_found, wet_found, xt_found, yt_found
character(len=128) :: name
type(fieldtype), dimension(:), allocatable :: fields
real, dimension(:,:), allocatable :: tmp
if(.not. file_exist(trim(file)) ) &
call mpp_error(FATAL, 'compare_grid: file '//trim(file)//' does not exist')
call mpp_open(unit,trim(file),MPP_RDONLY,MPP_NETCDF,threading=MPP_SINGLE,&
fileset=MPP_SINGLE)
call mpp_get_info(unit,ndim,nvar,natt,ntime)
allocate(fields(nvar))
call mpp_get_fields(unit,fields)
ht_found=.false.
wet_found=.false.
do i=1,nvar
call mpp_get_atts(fields(i),name=name,siz=siz_in)
select case (trim(name))
case('depth_t')
grid%ni = siz_in(1); grid%nj = siz_in(2)
allocate(grid%ht(grid%ni,grid%nj))
call mpp_read(unit,fields(i),grid%ht)
ht_found = .true.
case ('wet')
allocate(grid%wet(siz_in(1), siz_in(2)),tmp(siz_in(1), siz_in(2)))
call mpp_read(unit,fields(i),tmp)
grid%wet = tmp
deallocate(tmp)
wet_found = .true.
case('x_T')
allocate(grid%geolon_t(siz_in(1), siz_in(2)))
call mpp_read(unit,fields(i),grid%geolon_t)
xt_found = .true.
case('y_T')
allocate(grid%geolat_t(siz_in(1), siz_in(2)))
call mpp_read(unit,fields(i),grid%geolat_t)
yt_found = .true.
end select
enddo
if (.not. ht_found) call mpp_error(FATAL,'compare_grid: depth_t not found in the file '//trim(file) )
if (.not. wet_found) call mpp_error(FATAL,'compare_grid: wet not found in the file '//trim(file) )
if (.not. xt_found) call mpp_error(FATAL,'compare_grid: x_T not found in the file '//trim(file) )
if (.not. yt_found) call mpp_error(FATAL,'compare_grid: x_T not found in the file '//trim(file) )
call mpp_close(unit)
end subroutine read_grid
!#####################################################################
! compare grid size and geographic grid.
subroutine grid_check
integer :: i, j
if(grid_1%ni .ne. grid_2%ni .or. grid_1%nj .ne. grid_2%nj ) then
call mpp_error(FATAL,'compare_grid: grid sizes of '//trim(grid_file_1)// &
' and'//trim(grid_file_2)//' do not match')
endif
!--- make sure it is the same grid
do j = 1, grid_1%nj
do i = 1, grid_1%ni
if(grid_1%geolon_t(i,j) .ne. grid_2%geolon_t(i,j) .or. grid_1%geolat_t(i,j) .ne. grid_2%geolat_t(i,j)) then
write(stdout(),*) 'grid_1%geolon_t(',i,',',j,') = ', grid_1%geolon_t(i,j), &
'grid_2%geolon_t(',i,',',j,') = ', grid_2%geolon_t(i,j)
write(stdout(),*) 'grid_1%geolat_t(',i,',',j,') = ', grid_1%geolat_t(i,j), &
'grid_2%geolat_t(',i,',',j,') = ', grid_2%geolat_t(i,j)
call mpp_error(FATAL, 'file '//trim(grid_file_1)//' and file '//trim(grid_file_2)// &
' do not have the same geographical grid' )
endif
enddo
enddo
end subroutine grid_check
!#####################################################################
!--- compare the two grid file and write out the difference to file grid_edits
subroutine grid_compare
logical :: flag
integer :: i, j, unit
!--- check the depth_t
call mpp_open(unit,trim(grid_edits),MPP_OVERWR,MPP_ASCII,threading=&
MPP_SINGLE,fileset=MPP_SINGLE, nohdrs = .TRUE.)
flag = .false.
write(unit, *) '# i j depth_new'
write(unit, *) '#'
do i = 1, grid_1%ni
do j = 1, grid_2%nj
if (grid_1%ht(i,j) .ne. grid_2%ht(i,j)) then
write(unit, '(i6,a,i6,a,f24.15,a,f24.15)') i,',',j,',',grid_2%ht(i,j),' # was ',grid_1%ht(i,j)
flag = .true.
endif
enddo
enddo
if (.not. flag) write(unit,*) '# NO CHANGES in ht'
call mpp_close(unit)
!--- check land_sea mask
call mpp_open(unit,trim(mask_diff),MPP_OVERWR,MPP_ASCII,threading=&
MPP_SINGLE,fileset=MPP_SINGLE)
write(unit, *) '# i j mask_new'
write(unit, *) '#'
flag = .false.
do i = 1, grid_1%ni
do j = 1, grid_1%nj
if (grid_1%wet(i,j) .ne. grid_2%wet(i,j)) then
write(unit, '(i6,a,i6,a,i6,a,i6)') i,',',j,',',grid_2%wet(i,j),' # was ', grid_1%wet(i,j)
flag = .true.
endif
enddo
enddo
if (.not. flag) write(unit,*) '# NO CHANGES in wet'
call mpp_close(unit)
end subroutine grid_compare
!#####################################################################
end program compare_grid