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: tikal $' !---------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