!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! 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 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 17.0 2009/07/21 03:22:35 fms Exp $'
character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $'
!---------------------------------------------------------------------
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.
!--- 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)
! 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
!--- 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
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
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)
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.hu_found) call mpp_error(FATAL,'edit_grid: depth_c 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) )
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