FMS  2025.04
Flexible Modeling System
mosaic2.inc
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 
19 !> @file
20 
21 !> @brief Get exchange grid information from mosaic xgrid file.
22 !> Example usage:
23 !!
24 !! call get_mosaic_xgrid(fileobj, nxgrid, i1, j1, i2, j2, area)
25 !!
26  subroutine get_mosaic_xgrid_(fileobj, i1, j1, i2, j2, area, ibegin, iend)
27  type(FmsNetcdfFile_t), intent(in) :: fileobj !> The file that contains exchange grid information.
28  integer, intent(inout) :: i1(:), j1(:), i2(:), j2(:) !> i and j indices for grids 1 and 2
29  real(kind=fms_mos_kind_), intent(inout) :: area(:) !> area of the exchange grid. The area is scaled to
30  !! represent unit earth area
31  integer, optional, intent(in) :: ibegin, iend
32 
33  integer :: start(4), nread(4), istart
34  real(kind=fms_mos_kind_), dimension(2, size(i1(:))) :: tile1_cell, tile2_cell
35  integer :: nxgrid, n
36  real(kind=r8_kind) :: garea
37  real(kind=r8_kind) :: get_global_area
38 
39  garea = get_global_area() !< get_global_area returns a r8_kind
40 
41  ! When start and nread present, make sure nread(1) is the same as the size of the data
42  if(present(ibegin) .and. present(iend)) then
43  istart = ibegin
44  nxgrid = iend - ibegin + 1
45  if(nxgrid .NE. size(i1(:))) call mpp_error(fatal, .NE."get_mosaic_xgrid: nxgrid size(i1(:))")
46  if(nxgrid .NE. size(j1(:))) call mpp_error(fatal, .NE."get_mosaic_xgrid: nxgrid size(j1(:))")
47  if(nxgrid .NE. size(i2(:))) call mpp_error(fatal, .NE."get_mosaic_xgrid: nxgrid size(i2(:))")
48  if(nxgrid .NE. size(j2(:))) call mpp_error(fatal, .NE."get_mosaic_xgrid: nxgrid size(j2(:))")
49  if(nxgrid .NE. size(area(:))) call mpp_error(fatal, .NE."get_mosaic_xgrid: nxgrid size(area(:))")
50  else
51  istart = 1
52  nxgrid = size(i1(:))
53  endif
54 
55  start = 1; nread = 1
56  start(1) = istart; nread(1) = nxgrid
57 
58  call read_data(fileobj, 'xgrid_area', area, corner=start, edge_lengths=nread)
59 
60  start = 1; nread = 1
61  nread(1) = 2
62  start(2) = istart; nread(2) = nxgrid
63 
64  call read_data(fileobj, 'tile1_cell', tile1_cell, corner=start, edge_lengths=nread)
65  call read_data(fileobj, 'tile2_cell', tile2_cell, corner=start, edge_lengths=nread)
66 
67  do n = 1, nxgrid
68  i1(n) = int(tile1_cell(1,n))
69  j1(n) = int(tile1_cell(2,n))
70  i2(n) = int(tile2_cell(1,n))
71  j2(n) = int(tile2_cell(2,n))
72  area(n) = real( real(area(n),r8_kind)/garea, fms_mos_kind_ )
73  end do
74 
75  return
76 
77  end subroutine get_mosaic_xgrid_
78  !###############################################################################
79  !> @brief Calculate grid cell area.
80  !> Calculate the grid cell area. The purpose of this routine is to make
81  !! sure the consistency between model grid area and exchange grid area.
82  !> @param lon geographical longitude of grid cell vertices.
83  !> @param lat geographical latitude of grid cell vertices.
84  !> @param[inout] area grid cell area.
85  !> <br>Example usage:
86  !! call calc_mosaic_grid_area(lon, lat, area)
87  subroutine calc_mosaic_grid_area_(lon, lat, area)
88  real(kind=fms_mos_kind_), dimension(:,:), intent(in) :: lon
89  real(kind=fms_mos_kind_), dimension(:,:), intent(in) :: lat
90  real(kind=fms_mos_kind_), dimension(:,:), intent(out) :: area
91  integer :: nlon, nlat
92 
93  real(r8_kind) :: area_r8(size(area,1),size(area,2))
94 
95  nlon = size(area,1)
96  nlat = size(area,2)
97  ! make sure size of lon, lat and area are consitency
98  if( size(lon,1) .NE. nlon+1 .OR. size(lat,1) .NE. nlon+1 ) &
99  call mpp_error(fatal, "mosaic_mod: size(lon,1) and size(lat,1) should equal to size(area,1)+1")
100  if( size(lon,2) .NE. nlat+1 .OR. size(lat,2) .NE. nlat+1 ) &
101  call mpp_error(fatal, "mosaic_mod: size(lon,2) and size(lat,2) should equal to size(area,2)+1")
102 
103  ! get_grid_area only accepts double precision data
104  call get_grid_area( nlon, nlat, real(lon,r8_kind), real(lat,r8_kind), area_r8)
105 
106  area=real(area_r8,fms_mos_kind_)
107 
108  end subroutine calc_mosaic_grid_area_
109  !###############################################################################
110  !> @brief Calculate grid cell area using great cirlce algorithm
111  !> Calculate the grid cell area. The purpose of this routine is to make
112  !! sure the consistency between model grid area and exchange grid area.
113  !> @param lon geographical longitude of grid cell vertices.
114  !> @param lat geographical latitude of grid cell vertices.
115  !> @param[inout] area grid cell area.
116  !> <br>Example usage:
117  !! call calc_mosaic_grid_great_circle_area(lon, lat, area)
118  subroutine calc_mosaic_grid_great_circle_area_(lon, lat, area)
119  real(kind=fms_mos_kind_), dimension(:,:), intent(in) :: lon
120  real(kind=fms_mos_kind_), dimension(:,:), intent(in) :: lat
121  real(kind=fms_mos_kind_), dimension(:,:), intent(inout) :: area
122  integer :: nlon, nlat
123 
124  real(r8_kind) :: area_r8(size(area,1),size(area,2))
125 
126  nlon = size(area,1)
127  nlat = size(area,2)
128  ! make sure size of lon, lat and area are consitency
129  if( size(lon,1) .NE. nlon+1 .OR. size(lat,1) .NE. nlon+1 ) &
130  call mpp_error(fatal, "mosaic_mod: size(lon,1) and size(lat,1) should equal to size(area,1)+1")
131  if( size(lon,2) .NE. nlat+1 .OR. size(lat,2) .NE. nlat+1 ) &
132  call mpp_error(fatal, "mosaic_mod: size(lon,2) and size(lat,2) should equal to size(area,2)+1")
133 
134  ! get_grid_great_circle_area only accepts r8_kind arguments
135  call get_grid_great_circle_area( nlon, nlat, real(lon,r8_kind), real(lat,r8_kind), area_r8)
136 
137  area=real(area_r8, fms_mos_kind_)
138 
139  end subroutine calc_mosaic_grid_great_circle_area_
140  !#####################################################################
141  !> This function check if a point (lon1,lat1) is inside a polygon (lon2(:), lat2(:))
142  !! lon1, lat1, lon2, lat2 are in radians.
143  function is_inside_polygon_(lon1, lat1, lon2, lat2 )
144  real(kind=fms_mos_kind_), intent(in) :: lon1, lat1
145  real(kind=fms_mos_kind_), intent(in) :: lon2(:), lat2(:)
146  logical :: IS_INSIDE_POLYGON_
147  integer :: npts, isinside
148  integer :: inside_a_polygon
149 
150  npts = size(lon2(:))
151 
152  !> inside_a_polygon function only accepts r8_kind real variables
153 
154  isinside = inside_a_polygon(real(lon1,r8_kind), real(lat1,r8_kind), npts, real(lon2,r8_kind), real(lat2,r8_kind))
155  if(isinside == 1) then
156  is_inside_polygon_ = .true.
157  else
158  is_inside_polygon_ = .false.
159  endif
160 
161  return
162 
163  end function is_inside_polygon_
164 !> @}