!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! 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 !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module atmos_grid_mod !----------------------------------------------------------------------- ! 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 ! ! ! atmos_grid_mod Generate horizontal grid ( either bgrid or spectral grid ) ! for atmosphere and land model. ! ! ! There are four subgrids, labeled T (for tracer), C (corner of T), N (north of T) ! and E (east of T). The following schematic describes the grid cell notation. ! ! !
  !                          Ni,j
  !               +----------+-----------+Ci,j
  !               |                      |     
  !               |                      |
  !               |                      |
  !               +          +Ti,j       +Ei,j
  !               |                      |
  !               |                      |
  !               +----------+-----------+
  !
  !                          
  !
! ! The grid_spec file would contains all of the following information. ! !
  !   x_T, y_T           = Geographic location of T-cell center
  !   x_vert_T, y_vert_T = Geographic location of T-cell vertices(each cell has 4 vertices)
  !   x_C, y_C           = Geographic location of C-cell center
  !
  !
! !
use mpp_mod, only : mpp_pe, mpp_root_pe, mpp_chksum, mpp_error, FATAL, NOTE use mpp_io_mod, only : mpp_write_meta, mpp_write, axistype, fieldtype use fms_mod, only : write_version_number, file_exist, stdout, field_size, string use fms_mod, only : open_namelist_file, close_file, check_nml_error, read_data use transforms_mod, only : transforms_init, get_grid_boundaries, get_deg_lon, get_deg_lat use constants_mod, only : PI, RADIUS implicit none private !------ namelist interface --------------------------------------------- ! ! ! type of grid. Its value can be "bgrid" or "spectral". Its default value is "bgrid". ! ! ! control standard output. ! ! ! number of longitude points. ! ! ! number of latitude points. ! ! ! the fourier wavenumbers in the truncation are set equal to fourier_inc*m, where m = 0, ! num_fourier therefore, the total number of fourier modes is num_fourier +1. ! This namelist option is only for grid_type = 'spectral' ! ! ! the wavenumber increment (see num_fourier above). This namelist option is only ! for grid_type = 'spectral' ! ! ! the maximum meridional wavenumber Retained meridional wavewnumbers are n = 0, ! num_spherical The total spherical wavenumber is L = n+m. This namelist option is only ! for grid_type = 'spectral' ! ! ! range of the longitude. Default value is : lon_begin=0, lon_end=360. If you want to ! generate regional grid, you may need to set these namelists. Otherwise use default value. ! This namelist option is only for grid_type = 'bgrid' ! ! ! range of the latitude. Default value is : lat_begin=-90, lat_end=90. If you want to ! generate regional grid, you may need to set these namelists. Otherwise use default value. ! This namelist option is only for grid_type = 'bgrid' ! ! character(len=24) :: grid_type = 'bgrid' logical :: debug = .false. integer :: num_lon = 0 integer :: num_lat = 0 integer :: num_fourier = 0 integer :: num_spherical = 0 integer :: fourier_inc = 1 real :: lon_begin = 0.0 real :: lon_end = 360.0 real :: lat_begin = -90.0 real :: lat_end = 90.0 namelist /atmos_grid_nml/ debug, grid_type, num_lon, num_lat, num_fourier, num_spherical, & fourier_inc, lon_begin, lon_end, lat_begin, lat_end ! type atmos_grid_type real, dimension(:,:), pointer :: x_T => NULL() ! geographical longitude of T-cell center real, dimension(:,:), pointer :: y_T => NULL() ! geographical latitude of T-cell center real, dimension(:,:), pointer :: x_C => NULL() ! geographical longitude of C-cell center real, dimension(:,:), pointer :: y_C => NULL() ! geographical latitude of C-cell center real, dimension(:,:,:), pointer :: x_vert_T => NULL() ! geographical longitude of T-cell vertices real, dimension(:,:,:), pointer :: y_vert_T => NULL() ! geographical latitude of T-cell vertices end type atmos_grid_type ! !----------------------------------------------------------------------- !--------private data--------------------------------------------------- logical :: module_is_initialized = .false. !--- axis and field data type for output file ----------------------- type(axistype),save :: axis_xt, axis_yt, axis_xc, axis_yc, axis_v type(fieldtype),save :: fld_x_T, fld_y_T, fld_x_vert_T, fld_y_vert_T, fld_x_C, fld_y_C !---------version information------------------------------------------- character(len=128) :: version = '$Id: atmos_grid.f90,v 10.0 2003/10/24 22:01:43 fms Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' !---------public interface---------------------------------------------- public :: generate_atmos_grid, atmos_grid_init, atmos_grid_end, atmos_grid_type public :: write_atmos_grid_meta, write_atmos_grid_data contains !####################################################################### ! ! ! Initialization routine. ! ! ! Read namelist, write out version and namelist informaiton. ! ! ! subroutine atmos_grid_init !--- local variables ------------------------------------------------- integer :: unit, ierr, io !---- read namelist -------------------------------------------------- if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr=1 do while (ierr /= 0) read (unit, nml=atmos_grid_nml, iostat=io, end=10) ierr = check_nml_error(io,'atmos_grid_nml') ! also initializes nml error codes enddo 10 call close_file (unit) else call mpp_error(FATAL,'grid_generator: file input.nml does not exist') endif !--- write version info and namelist to logfile ---------------------- call write_version_number(version,tagname) write (stdout(), nml=atmos_grid_nml) module_is_initialized = .true. !--- check namelist value if(trim(grid_type) .ne. 'bgrid' .and. trim(grid_type) .ne. 'spectral') & call mpp_error(FATAL,'atmos_grid_mod: '//trim(grid_type)//' is not a valid option of nml "grid_type" ') if(num_lon .le. 0 .or. num_lat .le. 0) then call mpp_error(FATAL, 'atmos_grid_mod: nml num_lon= '//trim(string(num_lon))//' and num_lat= ' & //trim(string(num_lat))//' should be positive number') endif if(trim(grid_type) == 'spectral') then if(num_fourier == 0 .or. num_spherical == 0) & call mpp_error(FATAL, 'atmos_grid_mod: nml "num_fourier"= '//trim(string(num_fourier))// & ' and "num_spherical"= '// trim(string(num_spherical))//' shoulbe be positive number') else ! grid_type == bgrid if(lon_begin .ge. lon_end ) & call mpp_error(FATAL, 'atmos_grid_mod: nml lon_begin= '//trim(string(lon_begin))// & ' should be less than nml "lon_end"= '//trim(string(lon_end)) ) if(lat_begin .ge. lat_end ) & call mpp_error(FATAL, 'atmos_grid_mod: nml lat_begin= '//trim(string(lat_begin))// & ' should be less than nml "lat_end"= '//trim(string(lat_end)) ) endif end subroutine atmos_grid_init !####################################################################### ! ! ! Generate horizontal grid. ! ! ! Calculate geographic locations of T and C-cell center. ! ! ! ! A derived-type variable that contains horizontal grid information. ! ! subroutine generate_atmos_grid (Hgrid) type(atmos_grid_type), intent(inout) :: Hgrid select case(trim(grid_type)) case('bgrid') call generate_bgrid (Hgrid) case('spectral') call generate_spectral_grid (Hgrid) end select end subroutine generate_atmos_grid !####################################################################### !--- generate spectral horizontal grid subroutine generate_spectral_grid (Hgrid) type(atmos_grid_type), intent(inout) :: Hgrid real :: dlon, dlat, R2D integer :: i, j, pe, root_pe integer :: chksum_x_t, chksum_y_t, chksum_x_c, chksum_y_c, chksum_x_vert_t, chksum_y_vert_t real, allocatable, dimension(:) :: lonb, latb, deg_lon, deg_lat R2D = 180./PI pe = mpp_pe() root_pe = mpp_root_pe() call transforms_init(RADIUS, num_lat, num_lon, num_fourier, fourier_inc, num_spherical) allocate(lonb(num_lon+1), latb(num_lat+1)) allocate(deg_lon(num_lon), deg_lat(num_lat)) call get_deg_lon(deg_lon) call get_deg_lat(deg_lat) call get_grid_boundaries(lonb, latb, global=.true.) allocate(hgrid%x_C(num_lon,num_lat), hgrid%y_C(num_lon,num_lat), & hgrid%x_T(num_lon,num_lat), hgrid%y_T(num_lon,num_lat), & hgrid%x_vert_T(num_lon,num_lat,4), hgrid%y_vert_T(num_lon,num_lat,4)) do i=1,num_lon hgrid%x_C(i,:) = lonb(i+1)*R2D hgrid%x_T(i,:) = deg_lon(i) hgrid%x_vert_T(i,:,1) = lonb(i)*R2D hgrid%x_vert_T(i,:,2) = lonb(i+1)*R2D hgrid%x_vert_T(i,:,3) = lonb(i+1)*R2D hgrid%x_vert_T(i,:,4) = lonb(i)*R2D enddo do j=1,num_lat hgrid%y_C(:,j) = latb(j+1)*R2D hgrid%y_T(:,j) = deg_lat(j) hgrid%y_vert_T(:,j,1) = latb(j)*R2D hgrid%y_vert_T(:,j,2) = latb(j)*R2D hgrid%y_vert_T(:,j,3) = latb(j+1)*R2D hgrid%y_vert_T(:,j,4) = latb(j+1)*R2D enddo ! Print all grid coordinates write (stdout(),'(//,36x,a,/)') 'H G R I D G E N E R A T I O N' write (stdout(),'(/)') write (stdout(),9105) num_lat write (stdout(),9001) (hgrid%y_T(1,j),j=1,num_lat) write (stdout(),9106) num_lat write (stdout(),9001) (hgrid%y_C(1,j),j=1,num_lat) write (stdout(),9107) num_lon write (stdout(),9001) (hgrid%x_T(i,1),i=1,num_lon) write (stdout(),9108) num_lon write (stdout(),9001) (hgrid%x_T(i,1),i=1,num_lon) write (stdout(),'(/)') if(debug) then chksum_x_t = mpp_chksum(Hgrid%x_T,(/pe/)) chksum_y_t = mpp_chksum(Hgrid%y_T,(/pe/)) chksum_x_c = mpp_chksum(Hgrid%x_C,(/pe/)) chksum_y_c = mpp_chksum(Hgrid%y_C,(/pe/)) chksum_x_vert_t = mpp_chksum(Hgrid%x_vert_T,(/pe/)) chksum_y_vert_t = mpp_chksum(Hgrid%y_vert_T,(/pe/)) write (stdout(),'(/(a,I20/))')'spectral_hgrid chksum: chksum_x_T = ', chksum_x_t write (stdout(),'(/(a,I20/))')'spectral_hgrid chksum: chksum_y_T = ', chksum_y_t write (stdout(),'(/(a,I20/))')'spectral_hgrid chksum: chksum_x_C = ', chksum_x_c write (stdout(),'(/(a,I20/))')'spectral_hgrid chksum: chksum_y_c = ', chksum_y_c write (stdout(),'(/(a,I20/))')'spectral_hgrid chksum: chksum_x_vert_T = ', chksum_x_vert_t write (stdout(),'(/(a,I20/))')'spectral_hgrid chksum: chksum_y_vert_T = ', chksum_y_vert_t endif deallocate(lonb, latb, deg_lon, deg_lat) return 9001 format (1x,10f10.4) 9105 format (/,' Latitude of T points (deg): yt(j) j=1,',i4) 9106 format (/,' Latitude of C points (deg): yu(j) j=1,',i4) 9107 format (/,' Longitude of T points (deg): xt(i) i=1,',i4) 9108 format (/,' Longitude of C points (deg): xu(i) i=1,',i4) end subroutine generate_spectral_grid !##################################################################### !--- generate horizontal bgrid subroutine generate_bgrid (hgrid) type(atmos_grid_type), intent(inout) :: hgrid real :: dlon, dlat integer :: i, j, pe integer :: chksum_x_t, chksum_y_t, chksum_x_c, chksum_y_c, chksum_x_vert_t, chksum_y_vert_t pe = mpp_pe() allocate(hgrid%x_C(num_lon,num_lat), hgrid%y_C(num_lon,num_lat), & hgrid%x_T(num_lon,num_lat), hgrid%y_T(num_lon,num_lat), & hgrid%x_vert_T(num_lon,num_lat,4), hgrid%y_vert_T(num_lon,num_lat,4)) dlon = (lon_end-lon_begin)/num_lon dlat = (lat_end-lat_begin)/num_lat do i=1,num_lon hgrid%x_C(i,:) = lon_begin + dlon*real(i) hgrid%x_T(i,:) = lon_begin + dlon*(i-.5) hgrid%x_vert_T(i,:,1) = lon_begin + dlon*(i-1.0) hgrid%x_vert_T(i,:,2) = lon_begin + dlon*real(i) hgrid%x_vert_T(i,:,3) = lon_begin + dlon*real(i) hgrid%x_vert_T(i,:,4) = lon_begin + dlon*(i-1.0) enddo do j=1,num_lat hgrid%y_C(:,j) = lat_begin + dlat*real(j) hgrid%y_T(:,j) = lat_begin + dlat*(j-.5) hgrid%y_vert_T(:,j,1) = lat_begin + dlat*real(j-1.0) hgrid%y_vert_T(:,j,2) = lat_begin + dlat*real(j-1.0) hgrid%y_vert_T(:,j,3) = lat_begin + dlat*real(j) hgrid%y_vert_T(:,j,4) = lat_begin + dlat*real(j) enddo ! Print all grid coordinates write (stdout(),'(//,36x,a,/)') 'H G R I D G E N E R A T I O N' write (stdout(),'(/)') write (stdout(),9105) num_lat write (stdout(),9001) (hgrid%y_T(1,j),j=1,num_lat) write (stdout(),9106) num_lat write (stdout(),9001) (hgrid%y_C(1,j),j=1,num_lat) write (stdout(),9107) num_lon write (stdout(),9001) (hgrid%x_T(i,1),i=1,num_lon) write (stdout(),9108) num_lon write (stdout(),9001) (hgrid%x_C(i,1),i=1,num_lon) write (stdout(),'(/)') if(debug ) then chksum_x_t = mpp_chksum(Hgrid%x_T,(/pe/)) chksum_y_t = mpp_chksum(Hgrid%y_T,(/pe/)) chksum_x_c = mpp_chksum(Hgrid%x_C,(/pe/)) chksum_y_c = mpp_chksum(Hgrid%y_C,(/pe/)) chksum_x_vert_t = mpp_chksum(Hgrid%x_vert_T,(/pe/)) chksum_y_vert_t = mpp_chksum(Hgrid%y_vert_T,(/pe/)) write (stdout(),'(/(a,I20/))')'bgrid_hgrid chksum: chksum_x_T = ', chksum_x_t write (stdout(),'(/(a,I20/))')'bgrid_hgrid chksum: chksum_y_T = ', chksum_y_t write (stdout(),'(/(a,I20/))')'bgrid_hgrid chksum: chksum_x_C = ', chksum_x_c write (stdout(),'(/(a,I20/))')'bgrid_hgrid chksum: chksum_y_c = ', chksum_y_c write (stdout(),'(/(a,I20/))')'bgrid_hgrid chksum: chksum_x_vert_T = ', chksum_x_vert_t write (stdout(),'(/(a,I20/))')'bgrid_hgrid chksum: chksum_y_vert_T = ', chksum_y_vert_t endif return 9001 format (1x,10f10.4) 9105 format (/,' Latitude of T points (deg): yt(j) j=1,',i4) 9106 format (/,' Latitude of C points (deg): yu(j) j=1,',i4) 9107 format (/,' Longitude of T points (deg): xt(i) i=1,',i4) 9108 format (/,' Longitude of C points (deg): xu(i) i=1,',i4) end subroutine generate_bgrid !####################################################################### ! ! ! write the Hgrid data to netcdf file ! ! ! ! The unit corresponding the output netcdf file. Always is returned by mpp_open. ! ! ! A derived-type variable that contains horizontal grid information. ! ! subroutine write_atmos_grid_data (unit,Hgrid) integer, intent(in) :: unit type(atmos_grid_type), intent(in) :: Hgrid !--- write out axis data call mpp_write(unit, axis_xt) call mpp_write(unit, axis_yt) call mpp_write(unit, axis_xc) call mpp_write(unit, axis_yc) call mpp_write(unit, axis_v) call mpp_write(unit, fld_x_T, hgrid%x_T) call mpp_write(unit, fld_y_T, hgrid%y_T) call mpp_write(unit, fld_x_C, hgrid%x_C) call mpp_write(unit, fld_y_C, hgrid%y_C) call mpp_write(unit, fld_x_vert_T, hgrid%x_vert_T) call mpp_write(unit, fld_y_vert_T, hgrid%y_vert_T) return end subroutine write_atmos_grid_data !####################################################################### ! ! ! Write out horizontal grid meta data. ! ! ! ! The unit corresponding the output netcdf file. Always is returned by mpp_open. ! ! ! A derived-type variable that contains horizontal grid information. ! ! subroutine write_atmos_grid_meta(unit, Hgrid ) integer, intent(in) :: unit type(atmos_grid_type), intent(in) :: Hgrid !--- write out axis meta ---------------------------------------- call mpp_write_meta(unit, axis_xt,'grid_x_T','degree_east','Nominal Longitude of T-cell center', & cartesian ='X', data = hgrid%x_T(:,1) ) call mpp_write_meta(unit, axis_yt,'grid_y_T','degree_north','Nominal Latitude of T-cell center', & cartesian ='Y', data = hgrid%y_T(1,:) ) call mpp_write_meta(unit, axis_xc,'grid_x_C','degree_east','Nominal Longitude of C-cell center', & cartesian ='X', data = hgrid%x_C(:,1) ) call mpp_write_meta(unit, axis_yc,'grid_y_C','degree_north','Nominal Latitude of C-cell center', & cartesian ='Y', data = hgrid%y_C(1,:) ) call mpp_write_meta(unit, axis_v, 'vertex', 'none ', 'Vertex position from southwest couterclockwise', & data = (/1.,2.,3.,4./) ) call mpp_write_meta(unit, fld_x_T, (/axis_xt, axis_yt/), 'x_T', 'degree_east', & 'Geographic longitude of T_cell centers', pack=1) call mpp_write_meta(unit, fld_y_T, (/axis_xt, axis_yt/), 'y_T', 'degree_north', & 'Geographic latitude of T_cell centers', pack=1) call mpp_write_meta(unit, fld_x_C, (/axis_xc, axis_yc/), 'x_C', 'degree_east', & 'Geographic longitude of C_cell centers', pack=1) call mpp_write_meta(unit, fld_y_C, (/axis_xc, axis_yc/), 'y_C', 'degree_north', & 'Geographic latitude of C_cell centers', pack=1) call mpp_write_meta(unit, fld_x_vert_T, (/axis_xt, axis_yt, axis_v/), 'x_vert_T', & 'degree_east', 'Geographic longitude of T_cell vertices begin southwest counterclockwise', pack=1) call mpp_write_meta(unit, fld_y_vert_T, (/axis_xt, axis_yt, axis_v/), 'y_vert_T', & 'degree_north', 'Geographic latitude of T_cell vertices begin southwest counterclockwise', pack=1) return end subroutine write_atmos_grid_meta !####################################################################### ! ! ! Destruction routine. ! ! ! Deallocates memory used by "atmos_grid_type" variables. ! ! ! ! A derived-type variable that contains horizontal grid information. ! ! subroutine atmos_grid_end(Hgrid) type(atmos_grid_type), intent(inout) :: Hgrid !--- release memory of module variables ------------------------------ deallocate(Hgrid%x_C, Hgrid%y_C, Hgrid%x_T, Hgrid%y_T, Hgrid%x_vert_T, Hgrid%y_vert_T) module_is_initialized = .false. return end subroutine atmos_grid_end !####################################################################### end module atmos_grid_mod