program edit_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
!
! edit grid topography.
!
!
! This program can edit the topography of input grid_spec file "orig_grid"
! according to the ascii input file "grid_edits". Then it will output the
! new grid_spec file "mod_grid". The program read file "grid_edits" line
! by line. Each line contains grid points position and new topography value
! of those grid points. The grid points position is specified by the grid
! index. You can specify a point or a region at each line. For example,
!
! 100, 60, 0
! will set the depth at point (100,60) to 0.
! 40:45, 30:34, 1000
! will set the depth at region ( index i from 54 to 50 and j from 30 to 34 ) to 1000.
!
!
!
use mpp_mod, only : mpp_error, FATAL, NOTE, mpp_pe, mpp_npes
use mpp_io_mod, only : mpp_open, mpp_close, mpp_read, mpp_write, mpp_write_meta
use mpp_io_mod, only : MPP_RDONLY, MPP_NETCDF, MPP_SINGLE, MPP_ASCII, MPP_OVERWR
use mpp_io_mod, only : axistype, fieldtype, atttype, mpp_get_atts, mpp_get_info
use mpp_io_mod, only : mpp_get_axes, mpp_get_fields, mpp_get_axis_data, mpp_copy_meta
use mpp_io_mod, only : mpp_get_att_char, mpp_get_att_name, mpp_get_att_real_scalar
use mpp_io_mod, only : mpp_get_att_type, mpp_get_att_real
use fms_mod, only : fms_init, fms_end, file_exist, close_file, stdout
use fms_mod, only : open_namelist_file, check_nml_error, write_version_number
use constants_mod, only : constants_init
use topog_mod, only : process_topo, show_deepest, set_topog_nml
implicit none
#include
!--- namelist interface
!
!
!
! original grid file
!
!
! output grid file after modification.
!
!
! input text file. Each line is in the format as
! "is:ie, js:je, depth", which means set the depth at region
! (is:ie, js:je) to value "depth". is and ie can be equal or
! different. js and je can be same or different.
!
!
! Control standard output. Default value is false.
!
!
character(len=128) :: orig_grid = 'orig_grid' ! original grid file
character(len=128) :: mod_grid = 'mod_grid ' ! modified grid file
character(len=128) :: grid_edits = 'grid_edits.txt' ! test file used to edit grid
logical :: debug = .FALSE. ! will fail at compilation if included
! in the namelist
namelist /edit_grid_nml/ orig_grid, mod_grid, grid_edits, debug
!--- version information variables -----------------------------------
character(len=128) :: version = '$Id: edit_grid.F90,v 19.0 2012/01/06 22:07:48 fms Exp $'
character(len=128) :: tagname = '$Name: tikal $'
!---------------------------------------------------------------------
logical :: tripolar_grid =.false. ! indicate the grid is tripolar grid or not.
logical :: cyclic_x =.false. ! true indicate cyclic in x-direction
logical :: cyclic_y =.false. ! true indicate cyclic in y-direction
logical :: full_cell =.false. ! do not generate partial bottom cells
logical :: fill_isolated_cells =.false. ! Do not allow non-advective tracer cells
logical :: dont_change_landmask =.false. ! Do not change land/sea mask when filling isolated cells
logical :: fill_shallow =.false. ! Make cells less than minimum depth land
logical :: deepen_shallow =.false. ! Make cells less than minimum depth equal to minimum depth
logical :: round_shallow =.false. ! Make cells land if depth is less than 1/2 mimumim depth,
! otherwise make ocean
logical :: adjust_topo =.false. ! adjust topography
logical :: fill_first_row =.false. ! make first row of ocean model all land points for ice model
integer :: kmt_min = 2 ! minimum number of vertical levels
real, dimension(:), allocatable :: zw ! vertical grid at T-cell bound
real, dimension(:,:), allocatable :: ht ! depth at T-cell center
real, dimension(:,:), allocatable :: hu ! depth at UV-cell center
real, dimension(:,:), allocatable :: geolon_t ! geographical longitude of T-cell center
real, dimension(:,:), allocatable :: geolat_t ! geographical latitude of T-cell center
real, dimension(:,:), allocatable :: dxte, dytn
real, dimension(:,:), allocatable :: wet ! land/sea mask
real, dimension(:,:), allocatable :: wet_c ! land/sea mask
integer, dimension(:,:), allocatable :: kmt ! number of vertical levels
integer, dimension(:,:), allocatable :: kmu ! number of vertical levels
character(len=64) :: topography ! type of topography
type(axistype), dimension(:), allocatable :: axes_in, axes_out
type(fieldtype), dimension(:), allocatable :: fields_in, fields_out
type(atttype), dimension(:), allocatable :: global_atts
integer :: ni, nj, nk ! grid size
integer :: i,j
integer :: unit_input ! corresponding to the file orig_grid.
logical :: hu_exist
!--- begin of the program
call fms_init
call constants_init
call edit_grid_init
!--- first read from orig_grid
call read_orig_grid
if(topography == 'idealized') then
call mpp_error(FATAL,'edit_grid: can not edit grid when topography == idealized')
endif
!--- read the grid edit file and edit ht
call reset_depth
!--- set topog_mod namelist option
call set_topog_nml(full_cell, fill_isolated_cells, dont_change_landmask, fill_shallow, &
deepen_shallow, round_shallow, adjust_topo, fill_first_row, kmt_min, debug )
! Compare "ht" to bottom of deepest model level.
if(debug) call show_deepest(zw, ht)
call process_topo(ht, kmt, zw, dxte, dytn, geolon_t, geolat_t, tripolar_grid, cyclic_x, cyclic_y)
if( hu_exist ) then
! MJH
do j=1,nj-1
do i=1,ni-1
hu(i,j) = min(ht(i,j),ht(i+1,j), ht(i,j+1),ht(i+1,j+1))
kmu(i,j) = min(kmt(i,j), kmt(i+1,j), kmt(i,j+1), kmt(i+1,j+1))
enddo
hu(ni,j) = hu(ni-1,j)
kmu(ni,j) = kmu(ni-1,j)
enddo
hu(:,nj) = hu(:,nj-1)
kmu(:,nj) = kmu(:,nj-1)
wet_c=0.0
where (hu .gt. 0.0) wet_c = 1.0
endif
!--- define land/sea mask according to depth
wet = 0.0
where(ht(1:ni,1:nj) .gt. 0.0 ) wet = 1.0
call write_mod_grid ()
call fms_end
contains
!#####################################################################
! --- read the namelist and write the version information to logfile. Also
! --- write the namelist to standard output
subroutine edit_grid_init
integer :: io_status, unit, ierr
!--- This program do not support parallel, so always run on 1 pe
if(mpp_npes() .gt. 1) call mpp_error(FATAL,'program edit_grid: set npes = 1 in the runscripts')
! provide for namelist over-ride of defaults
if(file_exist('input.nml')) then
unit = open_namelist_file()
read (unit,edit_grid_nml,IOSTAT=io_status)
write (stdout(),'(/)')
write (stdout(),edit_grid_nml)
ierr = check_nml_error(io_status, 'edit_grid_nml')
call close_file(unit)
else
call mpp_error(NOTE, 'edit_grid: file input.nml does not exist' )
endif
!--- write out version information ---------------------------------
call write_version_number(version, tagname)
end subroutine edit_grid_init
!#####################################################################
!--- read the grid information from the file orig_grid
subroutine read_orig_grid
integer :: ndim, nvar, natt, ntime, i, len
character(len=128) :: name
logical :: zw_found, gridlon_found, gridlat_found
logical :: ht_found, hu_found, dxte_found, kmu_found
logical :: dytn_found, xt_found, yt_found
real, dimension(:,:), allocatable :: tmp2d
if(.not. file_exist(trim(orig_grid))) &
call mpp_error(FATAL, 'file '//trim(orig_grid)//' does not exist')
call mpp_open(unit_input,trim(orig_grid),MPP_RDONLY,MPP_NETCDF,threading=MPP_SINGLE,&
fileset=MPP_SINGLE)
call mpp_get_info(unit_input,ndim,nvar,natt,ntime)
allocate (global_atts(natt))
call mpp_get_atts(unit_input,global_atts)
allocate(axes_in(ndim), axes_out(ndim) )
call mpp_get_axes(unit_input,axes_in)
axes_out = axes_in
allocate(fields_in(nvar), fields_out(nvar))
call mpp_get_fields(unit_input,fields_in)
fields_out = fields_in
do i=1,natt
select case (trim(mpp_get_att_name(global_atts(i))))
case ('x_boundary_type')
if (trim(mpp_get_att_char(global_atts(i))) == 'cyclic') then
cyclic_x = .true.
else
cyclic_x = .false.
endif
case ('y_boundary_type')
if (trim(mpp_get_att_char(global_atts(i))) == 'fold_north_edge') then
tripolar_grid = .true.
else if(trim(mpp_get_att_char(global_atts(i))) == 'cyclic') then
cyclic_y = .true.
endif
case ('topography')
topography = mpp_get_att_char(global_atts(i))
case ('full_cell')
full_cell = .true.
case ('fill_isolated_cells')
fill_isolated_cells = .true.
case ('dont_change_landmask')
dont_change_landmask = .true.
case ('fill_shallow')
fill_shallow = .true.
case ('deepen_shallow')
deepen_shallow = .true.
case ('round_shallow')
round_shallow = .true.
case ('kmt_min')
kmt_min = mpp_get_att_real_scalar(global_atts(i))
case ('adjust_topo')
adjust_topo = .true.
case ('fill_first_row')
fill_first_row = .true.
end select
enddo
do i=1,ndim
call mpp_get_atts(axes_in(i),name=name,len=len)
select case (trim(name))
case ('zb')
allocate(zw(len))
nk = len
call mpp_get_axis_data(axes_in(i),zw)
zw_found = .true.
case ('grid_x_T')
ni = len
gridlon_found = .true.
case ('grid_y_T')
nj = len
gridlat_found=.true.
end select
enddo
if (.not.zw_found ) call mpp_error(FATAL,'edit_grid: axis zb not found in the file '//trim(orig_grid) )
if (.not.gridlon_found ) call mpp_error(FATAL,'edit_grid: axis grid_x_T not found in the file '//trim(orig_grid) )
if (.not.gridlat_found ) call mpp_error(FATAL,'edit_grid: axis grid_y_T not found in the file '//trim(orig_grid) )
allocate( ht(0:ni+1,0:nj+1), hu(ni,nj),kmt(0:ni+1,0:nj+1), kmu(ni,nj), wet(ni,nj), wet_c(ni,nj) )
allocate( geolon_t(ni,nj), geolat_t(ni,nj), dxte(ni,nj), dytn(ni,nj) )
ht = 0.0
hu = 0.0
kmt = 0
kmu = 0
ht_found = .false.
hu_found = .false.
dxte_found = .false.
dytn_found = .false.
xt_found = .false.
yt_found = .false.
kmu_found = .false.
do i=1,nvar
call mpp_get_atts(fields_in(i),name=name)
select case (trim(name))
case('depth_t')
call mpp_read(unit_input,fields_in(i),ht(1:ni,1:nj) )
ht_found = .true.
case ('ds_01_21_E')
call mpp_read(unit_input,fields_in(i),dxte)
dxte_found = .true.
case ('ds_10_12_N')
call mpp_read(unit_input,fields_in(i),dytn)
dytn_found = .true.
case ('x_T')
call mpp_read(unit_input,fields_in(i),geolon_t(1:ni,1:nj))
xt_found = .true.
case ('y_T')
call mpp_read(unit_input,fields_in(i),geolat_t(1:ni,1:nj))
yt_found = .true.
case ('depth_c')
call mpp_read(unit_input,fields_in(i),hu(1:ni,1:nj))
hu_found = .true.
case ('num_levels_c')
allocate(tmp2d(ni,nj))
call mpp_read(unit_input,fields_in(i),tmp2d(1:ni,1:nj))
kmu(1:ni,1:nj) = tmp2d(1:ni,1:nj)
deallocate(tmp2d)
kmu_found = .true.
end select
enddo
if (.not.ht_found) call mpp_error(FATAL,'edit_grid: depth_t not found in the file '//trim(orig_grid) )
if (.not.dxte_found) call mpp_error(FATAL,'edit_grid: ds_01_21_E not found in the file '//trim(orig_grid) )
if (.not.dytn_found) call mpp_error(FATAL,'edit_grid: ds_10_12_N not found in the file '//trim(orig_grid) )
if (.not.xt_found) call mpp_error(FATAL,'edit_grid: x_T not found in the file '//trim(orig_grid) )
if (.not.yt_found) call mpp_error(FATAL,'edit_grid: y_T not found in the file '//trim(orig_grid) )
if ( hu_found .neqv. kmu_found ) &
call mpp_error(FATAL,'edit_grid: depth_c and num_levels_c should co-exist in the file '//trim(orig_grid) )
hu_exist = hu_found
if( hu_exist ) then
write(stdout(),*)"NOTE from edit_grid: field depth_c and num_levels exist in the file "//trim(orig_grid)//".", &
"Please be aware that the output topography should only be used in radiating open boundary condition model"
endif
end subroutine read_orig_grid
!#####################################################################
!--- write out the grid information after modification according to
! grid_edits to file mod_grid.
subroutine write_mod_grid
integer :: unit, i, n
real, dimension(:,:) , allocatable :: tmp2d
real, dimension(:,:,:) , allocatable :: tmp3d
integer, dimension(4) :: siz_in = 1
character(len=128) :: name
call mpp_open(unit,trim(mod_grid),MPP_OVERWR,MPP_NETCDF,threading=MPP_SINGLE,&
fileset=MPP_SINGLE)
! ----------------------------------------------------------------------
! write global atts
! ----------------------------------------------------------------------
do i=1,size(global_atts(:))
if (mpp_get_att_type(global_atts(i)) == NF_FLOAT) &
call mpp_write_meta(unit,mpp_get_att_name(global_atts(i)), &
rval=mpp_get_att_real(global_atts(i)))
if (mpp_get_att_type(global_atts(i)) == NF_CHAR) then
if (trim(mpp_get_att_name(global_atts(i))) == 'filename') then
call mpp_write_meta(unit,mpp_get_att_name(global_atts(i)),cval=trim(mod_grid))
else
call mpp_write_meta(unit,mpp_get_att_name(global_atts(i)), cval=trim(mpp_get_att_char(global_atts(i))))
end if
endif
enddo
!-----------------------------------------------------------------------
! write axis metadata
! ----------------------------------------------------------------------
do i=1,size(axes_out(:))
call mpp_copy_meta(unit,axes_out(i))
enddo
! ----------------------------------------------------------------------
! write variable metadata
! ----------------------------------------------------------------------
do i=1,size(fields_out(:))
call mpp_copy_meta(unit,fields_out(i))
enddo
! ----------------------------------------------------------------------
! write axis data
! ----------------------------------------------------------------------
do i=1,size(axes_out(:))
call mpp_write(unit,axes_out(i))
enddo
do i=1,size(fields_out(:))
call mpp_get_atts(fields_out(i),name=name,siz=siz_in)
select case (trim(name))
case ('depth_t')
call mpp_write(unit,fields_out(i),ht(1:ni,1:nj))
case ('depth_c')
call mpp_write(unit,fields_out(i),hu(1:ni,1:nj))
case ('num_levels')
allocate(tmp2d(ni,nj))
tmp2d = kmt(1:ni,1:nj)
call mpp_write(unit,fields_out(i),tmp2d )
deallocate(tmp2d)
case ('num_levels_c')
allocate(tmp2d(ni,nj))
tmp2d = kmu(1:ni,1:nj)
call mpp_write(unit,fields_out(i),tmp2d )
deallocate(tmp2d)
case ('wet')
call mpp_write(unit,fields_out(i),wet)
case ('wet_c')
call mpp_write(unit,fields_out(i),wet_c)
case default
if(siz_in(3) > 1 ) then
allocate(tmp3d(siz_in(1),siz_in(2), siz_in(3)) )
call mpp_read(unit_input,fields_in(i),tmp3d)
call mpp_write(unit,fields_out(i),tmp3d)
deallocate(tmp3d)
else
allocate(tmp2d(siz_in(1),siz_in(2)) )
call mpp_read(unit_input,fields_in(i),tmp2d)
call mpp_write(unit,fields_out(i),tmp2d)
deallocate(tmp2d)
endif
end select
enddo
call mpp_close(unit)
call mpp_close(unit_input)
end subroutine write_mod_grid
!#####################################################################
!--- reset depth according to file grid_edits
subroutine reset_depth
integer :: unit, i, j, is, ie, js, je
character(len=128) :: txt
logical :: flag
real :: ht_new
if(.not. file_exist(trim(grid_edits))) &
call mpp_error(FATAL, 'file '//trim(grid_edits)//' does not exist')
call mpp_open(unit,trim(grid_edits),MPP_RDONLY,MPP_ASCII,threading=MPP_SINGLE,fileset=MPP_SINGLE)
do
read(unit,'(a)',end=99,err=99) txt
call parse_edits(txt,is,ie,js,je,ht_new,flag)
if (flag) then
do j=js,je
do i=is,ie
if (i < 0 .or. i > ni .or. j < 0 .or. j > nj) call mpp_error(FATAL,'indices exceed grid bounds')
if (ht_new < 0.0 .or. ht_new > zw(nk)) call mpp_error(FATAL,'depth exceeds allowable depth')
write(stdout(),*) 'Resetting depth at location (i,j) : ',i,j,' = ',ht_new
ht(i,j) = ht_new
enddo
enddo
endif
enddo
99 call mpp_close(unit)
end subroutine reset_depth
!#####################################################################
!--- read each line of the file grid_edits
subroutine parse_edits(txt,is,ie,js,je,depth,flag)
character(len=*), intent(in) :: txt
integer, intent(inout) :: is,ie,js,je
real, intent(inout) :: depth
logical, intent(inout) :: flag
integer :: i1,i2, i3
character(len=128) :: txt2
flag = .true.
i1 = scan(txt,',')
txt2 = txt(1:i1-1)
if (i1 <= 0) goto 90
i2 = scan(txt2,':')
if (i2 <= 0) then
read(txt2,*,err=90) is
ie = is
else
read(txt2(1:i2-1),*,err=90) is
read(txt2(i2+1:),*,err=90) ie
endif
i2 = scan(txt(i1+1:),',')
txt2 = txt(i1+1:i1+i2-1)
if (i2 <= 0) goto 90
i3 = scan(txt2,':')
if (i3 <= 0) then
read(txt2,*,err=90) js
je = js
else
read(txt2(1:i3-1),*,err=90) js
read(txt2(i3+1:),*,err=90) je
endif
read(txt(i1+i2+1:),*,err=90) depth
return
90 continue
flag = .false.
return
end subroutine parse_edits
!#####################################################################
end program edit_grid