FMS  2025.03
Flexible Modeling System
get_grid_version.inc
1 !***********************************************************************
2 !* GNU Lesser General Public License
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* FMS is free software: you can redistribute it and/or modify it under
7 !* the terms of the GNU Lesser General Public License as published by
8 !* the Free Software Foundation, either version 3 of the License, or (at
9 !* your option) any later version.
10 !*
11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 !* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 !* for more details.
15 !*
16 !* You should have received a copy of the GNU Lesser General Public
17 !* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18 !***********************************************************************
19 
20 !> Get global lon and lat of three model (target) grids, with a given file name
21 subroutine get_grid_version_1_(grid_file, mod_name, domain, lon, lat, min_lon, max_lon)
22  integer, parameter :: lkind = fms_get_grid_version_kind_
23 
24  character(len=*), intent(in) :: grid_file !< name of grid file
25  character(len=*), intent(in) :: mod_name !< module name
26  type(domain2d), intent(in) :: domain !< 2D domain
27  real(lkind), intent(out), dimension(:,:), allocatable :: lon, lat
28  real(lkind), intent(out) :: min_lon, max_lon
29 
30  integer :: i, j, siz(4)
31  integer :: nlon, nlat !< size of global lon and lat
32  real(lkind), dimension(:,:,:), allocatable :: lon_vert, lat_vert !< of OCN grid vertices
33  real(lkind), dimension(:), allocatable :: glon, glat !< lon and lat of 1-D grid of atm/lnd
34  logical :: is_new_grid
35  integer :: is, ie, js, je
36  integer :: isg, ieg, jsg, jeg
37  integer :: isc, iec, jsc, jec
38  character(len=3) :: xname, yname
39  integer :: start(2), nread(2)
40  type(FmsNetcdfDomainFile_t) :: fileobj
41  integer :: ndims !< Number of dimensions
42 
43  if(.not. open_file(fileobj, grid_file, 'read', domain )) then
44  call mpp_error(fatal, 'data_override_mod(get_grid_version_1): Error in opening file '//trim(grid_file))
45  endif
46 
47  call mpp_get_global_domain(domain, isg, ieg, jsg, jeg)
48  call mpp_get_compute_domain(domain, isc, iec, jsc, jec)
49 
50  allocate(lon(isc:iec, jsc:jec))
51  allocate(lat(isc:iec, jsc:jec))
52 
53  select case(mod_name)
54  case('ocn', 'ice')
55  is_new_grid = .false.
56  if(variable_exists(fileobj, 'x_T')) then
57  is_new_grid = .true.
58  else if(variable_exists(fileobj, 'geolon_t')) then
59  is_new_grid = .false.
60  else
61  call mpp_error(fatal,'data_override: both x_T and geolon_t is not in the grid file '//trim(grid_file) )
62  endif
63 
64  if(is_new_grid) then
65  ndims = get_variable_num_dimensions(fileobj, 'x_T')
66  call get_variable_size(fileobj, 'x_T', siz(1:ndims))
67  nlon = siz(1); nlat = siz(2)
68  call check_grid_sizes(trim(mod_name)//'_domain ', domain, nlon, nlat)
69  allocate(lon_vert(isc:iec,jsc:jec,4), lat_vert(isc:iec,jsc:jec,4) )
70 
71  call read_data(fileobj, 'x_vert_T', lon_vert)
72  call read_data(fileobj, 'y_vert_T', lat_vert)
73 
74 !2 Global lon and lat of ocean grid cell centers are determined from adjacent vertices
75  lon(:,:) = (lon_vert(:,:,1) + lon_vert(:,:,2) + lon_vert(:,:,3) + lon_vert(:,:,4))*0.25_lkind
76  lat(:,:) = (lat_vert(:,:,1) + lat_vert(:,:,2) + lat_vert(:,:,3) + lat_vert(:,:,4))*0.25_lkind
77  else
78 
79  ndims = get_variable_num_dimensions(fileobj, 'geolon_vert_t')
80  call get_variable_size(fileobj, 'geolon_vert_t', siz(1:ndims))
81  nlon = siz(1) - 1; nlat = siz(2) - 1;
82  call check_grid_sizes(trim(mod_name)//'_domain ', domain, nlon, nlat)
83 
84  start(1) = isc; nread(1) = iec-isc+2
85  start(2) = jsc; nread(2) = jec-jsc+2
86 
87  allocate(lon_vert(isc:iec+1,jsc:jec+1,1))
88  allocate(lat_vert(isc:iec+1,jsc:jec+1,1))
89 
90  call read_data(fileobj, 'geolon_vert_t', lon_vert(:,:,1), corner=start, edge_lengths=nread)
91  call read_data(fileobj, 'geolat_vert_t', lat_vert(:,:,1), corner=start, edge_lengths=nread)
92 
93  do j = jsc, jec
94  do i = isc, iec
95  lon(i,j) = (lon_vert(i,j,1) + lon_vert(i+1,j,1) + &
96  lon_vert(i+1,j+1,1) + lon_vert(i,j+1,1))*0.25_lkind
97  lat(i,j) = (lat_vert(i,j,1) + lat_vert(i+1,j,1) + &
98  lat_vert(i+1,j+1,1) + lat_vert(i,j+1,1))*0.25_lkind
99  enddo
100  enddo
101  endif
102  deallocate(lon_vert)
103  deallocate(lat_vert)
104  case('atm', 'lnd')
105  if(trim(mod_name) == 'atm') then
106  xname = 'xta'; yname = 'yta'
107  else
108  xname = 'xtl'; yname = 'ytl'
109  endif
110  ndims = get_variable_num_dimensions(fileobj, xname)
111  call get_variable_size(fileobj, xname, siz(1:ndims))
112  nlon = siz(1); allocate(glon(nlon))
113  call read_data(fileobj, xname, glon)
114 
115  ndims = get_variable_num_dimensions(fileobj, xname)
116  call get_variable_size(fileobj, yname, siz(1:ndims))
117  nlat = siz(1); allocate(glat(nlat))
118  call read_data(fileobj, yname, glat)
119  call check_grid_sizes(trim(mod_name)//'_domain ', domain, nlon, nlat)
120 
121  is = isc - isg + 1; ie = iec - isg + 1
122  js = jsc - jsg + 1; je = jec - jsg + 1
123  do j = js, jec
124  do i = is, ie
125  lon(i,j) = glon(i)
126  lat(i,j) = glat(j)
127  enddo
128  enddo
129  deallocate(glon)
130  deallocate(glat)
131  case default
132  call mpp_error(fatal, "data_override_mod: mod_name should be 'atm', 'ocn', 'ice' or 'lnd' ")
133  end select
134 
135  call close_file(fileobj)
136 
137  ! convert from degree to radian
138  lon = lon * real(deg_to_rad, lkind)
139  lat = lat* real(deg_to_rad, lkind)
140  min_lon = minval(lon)
141  max_lon = maxval(lon)
142  call mpp_min(min_lon)
143  call mpp_max(max_lon)
144 end subroutine get_grid_version_1_
145 
146 !> Get global lon and lat of three model (target) grids from mosaic.nc.
147 !! Currently we assume the refinement ratio is 2 and there is one tile on each pe.
148 subroutine get_grid_version_2_(fileobj, mod_name, domain, lon, lat, min_lon, max_lon, use_center_grid_points)
149  integer, parameter :: lkind = fms_get_grid_version_kind_
150 
151  type(FmsNetcdfFile_t), intent(in) :: fileobj !< file object for grid file
152  character(len=*), intent(in) :: mod_name !< module name
153  type(domain2d), intent(in) :: domain !< 2D domain
154  real(lkind), intent(out), dimension(:,:), allocatable :: lon, lat
155  real(lkind), intent(out) :: min_lon, max_lon
156  logical, optional, intent(in) :: use_center_grid_points !< Flag indicating whether or not to use the
157  !! centroid values of the supergrid from the
158  !! grid file as opposed to calcuating it by
159  !! taking the average of the four corner points.
160  !! This is only relevant to OCN and ICE grids.
161 
162  integer :: i, j, siz(2)
163  integer :: nlon, nlat ! size of global grid
164  integer :: nlon_super, nlat_super ! size of global supergrid.
165  integer :: isc, iec, jsc, jec
166  integer :: isc2, iec2, jsc2, jec2
167  character(len=FMS_PATH_LEN) :: solo_mosaic_file, grid_file
168  real(lkind), allocatable :: tmpx(:,:), tmpy(:,:)
169  logical :: open_solo_mosaic
170  type(FmsNetcdfFile_t) :: mosaicfileobj, tilefileobj
171  integer :: start(2), nread(2)
172  logical :: use_center_grid_points_local
173 
174  use_center_grid_points_local = .false.
175  if (present(use_center_grid_points)) use_center_grid_points_local = use_center_grid_points
176 
177  if(trim(mod_name) .NE. 'atm' .AND. trim(mod_name) .NE. 'ocn' .AND. &
178  trim(mod_name) .NE. 'ice' .AND. trim(mod_name) .NE. 'lnd' ) call mpp_error(fatal, &
179  "data_override_mod: mod_name should be 'atm', 'ocn', 'ice' or 'lnd' ")
180 
181  call mpp_get_compute_domain(domain, isc, iec, jsc, jec)
182 
183  allocate(lon(isc:iec, jsc:jec))
184  allocate(lat(isc:iec, jsc:jec))
185 
186  ! get the grid file to read
187 
188  if(variable_exists(fileobj, trim(mod_name)//'_mosaic_file' )) then
189  call read_data(fileobj, trim(mod_name)//'_mosaic_file', solo_mosaic_file)
190 
191  solo_mosaic_file = 'INPUT/'//trim(solo_mosaic_file)
192  if(.not. open_file(mosaicfileobj, solo_mosaic_file, 'read')) then
193  call mpp_error(fatal, 'data_override_mod(get_grid_version_2): Error in opening solo mosaic file '// &
194  & trim(solo_mosaic_file))
195  endif
196  open_solo_mosaic=.true.
197  else
198  mosaicfileobj = fileobj
199  open_solo_mosaic = .false.
200  end if
201 
202  call get_mosaic_tile_grid(grid_file, mosaicfileobj, domain)
203 
204  if(.not. open_file(tilefileobj, grid_file, 'read')) then
205  call mpp_error(fatal, 'data_override_mod(get_grid_version_2): Error in opening tile file '//trim(grid_file))
206  endif
207 
208  call get_variable_size(tilefileobj, 'area', siz)
209  nlon_super = siz(1); nlat_super = siz(2)
210  if( mod(nlon_super,2) .NE. 0) call mpp_error(fatal, &
211  'data_override_mod: '//trim(mod_name)//' supergrid longitude size can not be divided by 2')
212  if( mod(nlat_super,2) .NE. 0) call mpp_error(fatal, &
213  'data_override_mod: '//trim(mod_name)//' supergrid latitude size can not be divided by 2')
214  nlon = nlon_super/2;
215  nlat = nlat_super/2;
216  call check_grid_sizes(trim(mod_name)//'_domain ', domain, nlon, nlat)
217  isc2 = 2*isc-1; iec2 = 2*iec+1
218  jsc2 = 2*jsc-1; jec2 = 2*jec+1
219 
220  start(1) = isc2; nread(1) = iec2-isc2+1
221  start(2) = jsc2; nread(2) = jec2-jsc2+1
222 
223  allocate(tmpx(isc2:iec2, jsc2:jec2), tmpy(isc2:iec2, jsc2:jec2) )
224 
225  call read_data( tilefileobj, 'x', tmpx, corner=start,edge_lengths=nread)
226  call read_data( tilefileobj, 'y', tmpy, corner=start,edge_lengths=nread)
227 
228  ! copy data onto model grid
229  if(trim(mod_name) == 'atm' .OR. trim(mod_name) == 'lnd' .OR. use_center_grid_points_local) then
230  do j = jsc, jec
231  do i = isc, iec
232  lon(i,j) = tmpx(i*2,j*2)
233  lat(i,j) = tmpy(i*2,j*2)
234  end do
235  end do
236  else
237  do j = jsc, jec
238  do i = isc, iec
239  lon(i,j) = (tmpx(i*2-1,j*2-1)+tmpx(i*2+1,j*2-1)+tmpx(i*2+1,j*2+1)+tmpx(i*2-1,j*2+1))*0.25_lkind
240  lat(i,j) = (tmpy(i*2-1,j*2-1)+tmpy(i*2+1,j*2-1)+tmpy(i*2+1,j*2+1)+tmpy(i*2-1,j*2+1))*0.25_lkind
241  end do
242  end do
243  endif
244 
245  ! convert to radian
246  lon = lon * real(deg_to_rad, lkind)
247  lat = lat * real(deg_to_rad, lkind)
248 
249  deallocate(tmpx, tmpy)
250  min_lon = minval(lon)
251  max_lon = maxval(lon)
252  call mpp_min(min_lon)
253  call mpp_max(max_lon)
254 
255  call close_file(tilefileobj)
256  if(open_solo_mosaic) call close_file(mosaicfileobj)
257 end subroutine get_grid_version_2_