!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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 create_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 !----------------------------------------------------------------------- ! ! Zhi Liang ! ! ! This program works together with program make_xgrids, runoff_regrid to ! remap runoff data from a spherical grid onto any grid. This program can ! generate a spherical grid netcdf file that can be used to generate exchange ! grid grid_spec.nc with the destination grid file ( ocean grid file). ! ! ! ! This program expects to get grid data from a netcdf file, which is specfied ! by the namelist variable "src_data". The name of the runoff field is specified ! by the namelist variable "src_fld_name". The output file is a netcdf file with ! name "atmos_grid.nc". ! ! implicit none #include "netcdf.inc" !--- namelist interface ---------------------------------------------- ! ! ! Name of input file containing runoff data. ! ! ! Name of runoff field. ! ! character(len=128) :: src_fld_name = 'runoff' ! name of the runoff field character(len=128) :: src_data = 'src_data.nc' ! name of the runoff source data file namelist /create_grid_nml/ src_fld_name, src_data !--- version information --------------------------------------------- character(len=128) :: version= '$ID$' character(len=128) :: tagname= '$Name: mom4p1_pubrel_dec2009_nnz $' !--- other variables ------------------------------------------------- real, allocatable :: lon(:), lat(:) integer :: ni, nj integer :: stdout = 6 !--- initialize routine call create_grid_init() !--- read source data file to get longitude and latitude ------------- call read_src_data() !--- generate the simple latlon atmosphere grid !--- the source data should be a simple latlon grid call generate_grid() contains !####################################################################### !--- initialization routine : read namelist subroutine create_grid_init integer :: unit=7, iostat logical :: opened !--- read namelist do inquire( unit, opened = opened ) if( .not. opened ) exit unit = unit + 1 if( unit.EQ.100 )call error_handler( 'Error: unable to locate unit number.' ) enddo open (unit, file='input.nml', iostat = iostat ) if(iostat .ne. 0) call error_handler('error in open input.nml file') read (unit, create_grid_nml, iostat=iostat) if(iostat .ne. 0) call error_handler('error in while reading create_grid_nml ') close(unit) !--- write out namelist and version information write( stdout,'(/a)' )'create_grid '//trim(version)//trim(tagname) write( stdout, create_grid_nml) end subroutine create_grid_init !####################################################################### !--- read the longitude and latitude from src_data subroutine read_src_data integer :: dims(4) integer :: rcode, ncid, nvar, id_runoff, i character(len=64) :: lon_name, lat_name, name rcode = nf_open(trim(src_data), NF_NOWRITE, ncid) call error_handler('error in opening file '//trim(src_data), rcode) rcode = nf_inq_nvars(ncid, nvar) call error_handler('error in inquring nvars', rcode) rcode = nf_inq_varid(ncid, trim(src_fld_name), id_runoff ) call error_handler('error in inquring variable id', rcode) rcode = nf_inq_vardimid(ncid, id_runoff, dims) call error_handler('error in inquring dimension id', rcode) rcode = nf_inq_dimlen(ncid, dims(1), ni) call error_handler('error in inquring dimension length ni', rcode) rcode = nf_inq_dimlen(ncid, dims(2), nj) call error_handler('error in inquring dimension length nj', rcode) allocate(lon(ni), lat(nj) ) rcode = nf_inq_dimname(ncid, dims(1), lon_name) call error_handler('error in inquring lon_name', rcode) rcode = nf_inq_dimname(ncid, dims(2), lat_name) call error_handler('error in inquring lat_name', rcode) do i = 1, nvar rcode = nf_inq_varname(ncid,i,name) if(trim(name) == trim(lon_name)) then rcode = nf_get_var_double(ncid,i,lon) else if(trim(name) == trim(lat_name)) then rcode = nf_get_var_double(ncid,i,lat) endif enddo end subroutine read_src_data !##################################################################### !--- generate the spherical grid subroutine generate_grid integer :: i, j, rcode, ncid, start(4), nwrite(4), dims(3), id_vertex integer :: id_x_t, id_y_t, id_x_vert_t, id_y_vert_t, id_grid_x_t, id_grid_y_t real, allocatable :: x_t(:,:), y_t(:,:), x_vert_t(:,:,:), y_vert_t(:,:,:) real, allocatable :: lonb(:), latb(:) allocate(x_t(ni, nj), x_vert_t(ni, nj,4) ) allocate(y_t(ni, nj), y_vert_t(ni, nj,4) ) allocate(lonb(ni+1), latb(nj+1)) do i = 1, ni x_t(i,:) = lon(i) enddo do j = 1, nj y_t(:,j) = lat(j) enddo lonb(1) = (3.0*lon(1) - lon(2))*0.5 lonb(ni+1) = (3.0*lon(ni) - lon(ni-1))*0.5 do i = 2, ni lonb(i) = 0.5*(lon(i-1) + lon(i)) enddo latb(1) = (3.0*lat(1) - lat(2)) * 0.5 latb(nj+1) = (3.0*lat(nj) - lat(nj-1)) *0.5 do j = 2, nj latb(j) = 0.5*(lat(j-1) + lat(j)) enddo do i = 1, ni x_vert_t(i,:,1) = lonb(i) x_vert_t(i,:,2) = lonb(i+1) x_vert_t(i,:,3) = lonb(i+1) x_vert_t(i,:,4) = lonb(i) enddo do j = 1, nj y_vert_t(:,j,1) = latb(j) y_vert_t(:,j,2) = latb(j) y_vert_t(:,j,3) = latb(j+1) y_vert_t(:,j,4) = latb(j+1) enddo rcode = nf_create('atmos_grid.nc', NF_WRITE, ncid) call error_handler('error in create file atmos_grid.nc', rcode) rcode = nf_def_dim(ncid, 'grid_x_T', ni, dims(1)) call error_handler('error in defining x dimension', rcode) rcode = nf_def_dim(ncid, 'grid_y_T', nj, dims(2)) call error_handler('error in defining y dimension', rcode) rcode = nf_def_dim(ncid, 'vertex', 4, dims(3)) call error_handler('error in defining vertex dimension', rcode) rcode = nf_def_var(ncid, 'grid_x_T', NF_FLOAT, 1, dims(1), id_grid_x_t) call error_handler('error in defining axis grid_x_T', rcode) rcode = nf_def_var(ncid, 'grid_y_T', NF_FLOAT, 1, dims(2), id_grid_y_t) call error_handler('error in defining axis grid_y_T', rcode) rcode = nf_def_var(ncid, 'vertex', NF_FLOAT, 1, dims(3), id_vertex) call error_handler('error in defining axis vertex', rcode) rcode = nf_def_var(ncid, 'x_T', NF_DOUBLE, 2, dims(1:2), id_x_t) call error_handler('error in defining field x_T', rcode) rcode = nf_def_var(ncid, 'y_T', NF_DOUBLE, 2, dims(1:2), id_y_t) call error_handler('error in defining field y_T', rcode) rcode = nf_def_var(ncid, 'x_vert_T', NF_DOUBLE, 3, dims(1:3), id_x_vert_t) call error_handler('error in defining field x_vert_T', rcode) rcode = nf_def_var(ncid, 'y_vert_T', NF_DOUBLE, 3, dims(1:3), id_y_vert_t) call error_handler('error in defining field y_vert_T', rcode) rcode = nf_enddef(ncid) start = 1; nwrite = 1 nwrite(1) = ni rcode = nf_put_vara_double(ncid, id_grid_x_t, start, nwrite, lon) call error_handler('error in putting the data for axis grid_x_T', rcode) nwrite(1) = nj rcode = nf_put_vara_double(ncid, id_grid_y_t, start, nwrite, lat) call error_handler('error in putting the data for axis grid_y_T', rcode) nwrite(1) = ni; nwrite(2) = nj rcode = nf_put_vara_double(ncid, id_x_t, start, nwrite, x_t) call error_handler('error in putting the data for field x_T', rcode) rcode = nf_put_vara_double(ncid, id_y_t, start, nwrite, y_t) call error_handler('error in putting the data for field y_T', rcode) nwrite(3) = 4 rcode = nf_put_vara_double(ncid, id_x_vert_t, start, nwrite, x_vert_t) call error_handler('error in putting the data for field x_vert_T', rcode) rcode = nf_put_vara_double(ncid, id_y_vert_t, start, nwrite, y_vert_t) call error_handler('error in put the data for field y_vert_T', rcode) rcode = nf_close(ncid) deallocate(x_t, y_t, x_vert_t, y_vert_t, lonb, latb, lon, lat) end subroutine generate_grid !##################################################################### ! error handling routine. subroutine error_handler(mesg, status) character(len=*), intent(in) :: mesg integer, optional, intent(in) :: status character(len=128) :: msg if(present(status)) then if(status == 0) return write(msg,*) mesg, ', status code is ', status else write(msg,*) mesg endif write(stdout,*) msg call ABORT() end subroutine error_handler !####################################################################### end program create_grid