test_fms/mosaic2/test_mosaic2.F90

   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 !> @brief  This programs tests calls to get_mosaic_ntiles, get_mosaic_ncontacts,
  21 !! get_mosaic_grid_sizes, get_mosaic_contact.  The subroutines are tested with
  22 !! made up C1 grids and exchange grids.  See write_files mod for grid details.
  23
  24 #include "write_files.inc" !> including write_files.mod because I don't know how to compile when write_files.mod is
  25                            !! in a separate file.
  26 program test_mosaic
  27
  28 use mosaic2_mod
  29 use grid2_mod
  30 use write_files
  31 use mpp_mod,       only : mpp_init, mpp_error, FATAL, mpp_pe, mpp_root_pe
  32 use fms2_io_mod,   only : open_file, close_file, FmsNetcdfFile_t, fms2_io_init, read_data
  33 use fms_mod,       only : fms_init, fms_end
  34 use constants_mod, only : DEG_TO_RAD
  35 use platform_mod,  only : r4_kind, r8_kind
  36
  37 implicit none
  38
  39 !> create mosaic and grid files
  40 !! In orderr to create the mosaic and grid files, fms2_io needs to be initialized first
  41 call fms2_io_init()
  42 call write_all()
  43 !< fms_init calls grid_init which reads in the grid_spec file
  44 !! In this case, the grid_version is VERSION_OCN_MOSAIC_FILE.
  45 call fms_init()
  46
  47 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_MOSAIC_GRID_SIZES'
  48 call test_get_mosaic_grid_sizes()
  49
  50 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_MOSAIC_CONTACT'
  51 call test_get_mosaic_contact()
  52
  53 !> does not work, results in negative areas for r4_kind.  Figure out why later
  54 !if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_GRID_GREAT_CIRCLE_AREA'
  55 !call test_get_grid_great_circle_area()
  56
  57 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST CALC_MOSAIC_GRID_AREA'
  58 call test_calc_mosaic_grid_area()
  59
  60 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_MOSAIC_XGRID'
  61 call test_get_mosaic_xgrid()
  62
  63 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST IS_INSIDE_POLYGON'
  64 call test_is_inside_polygon()
  65
  66 call fms_end()
  67
  68 contains
  69 !------------------------------------------------------!
  70 subroutine test_get_mosaic_grid_sizes
  71
  72   !> test get_mosaic_grid_sizes
  73
  74   integer              :: ntiles !< number of tiles
  75   integer              :: n      !< counter
  76   integer, allocatable :: nx_out(:), ny_out(:) !< number of grid points for each tile
  77   type(FmsNetcdfFile_t):: fileobj
  78
  79   !-- ocean --!
  80   if( .not. open_file(fileobj, 'INPUT/'//trim(ocn_mosaic_file), 'read') ) &
  81        call mpp_error(FATAL, 'test_mosaic: error in opening file '//'INPUT/'//trim(ocn_mosaic_file))
  82
  83   allocate( nx_out(ocn_ntiles), ny_out(ocn_ntiles) )
  84   !> get_mosaic_grid_sizes reads in the grid file
  85   call get_mosaic_grid_sizes(fileobj, nx_out, ny_out )
  86   do n=1, ocn_ntiles
  87      call check_answer(ocn_nx/2, nx_out(n), 'ocn TEST_GET_MOSAIC_GRID_SIZES')
  88      call check_answer(ocn_nY/2, ny_out(n), 'ocn TEST_GET_MOSAIC_GRID_SIZES')
  89   end do
  90   deallocate(nx_out, ny_out)
  91   call close_file(fileobj)
  92
  93   !-- atm --!
  94   if( .not. open_file(fileobj, 'INPUT/'//trim(c1_mosaic_file), 'read') ) &
  95        call mpp_error(FATAL, 'test_mosaic: error in opening file '//'INPUT/'//trim(c1_mosaic_file))
  96
  97   allocate( nx_out(c1_ntiles), ny_out(c1_ntiles) )
  98   call get_mosaic_grid_sizes(fileobj, nx_out, ny_out)
  99   ntiles = get_mosaic_ntiles(fileobj)
 100   do n=1, ntiles
 101      call check_answer(c1_nx/2, nx_out(n), 'atm TEST_GET_MOSAIC_GRID_SIZES')
 102      call check_answer(c1_nx/2, ny_out(n), 'atm TEST_GET_MOSAIC_GRID_SIZES')
 103   end do
 104   deallocate(nx_out, ny_out)
 105   call close_file(fileobj)
 106
 107 end subroutine test_get_mosaic_grid_sizes
 108 !------------------------------------------------------!
 109 subroutine test_get_mosaic_contact
 110
 111   !< @uriel.ramirez
 112
 113   integer              :: ntiles         !< Number of tiles
 114   integer              :: ncontacts      !< Number of contacts
 115   integer              :: n              !< For do loops
 116   integer, allocatable :: tile1(:)       !< tile number for first contact
 117   integer, allocatable :: tile2(:)       !< tile number of the second contact
 118   integer, allocatable :: nx(:), ny(:)   !< Number of x/y points for each tile
 119   integer, allocatable :: istart1(:), iend1(:), jstart1(:), jend1(:) !< Indexes of first contact point
 120   integer, allocatable :: istart2(:), iend2(:), jstart2(:), jend2(:) !< Indexes of second contact point
 121
 122   integer              :: answers(2, 8)  !< Expected results
 123
 124   type(FmsNetcdfFile_t):: ocn_fileobj
 125
 126   if( .not. open_file(ocn_fileobj, 'INPUT/'//trim(ocn_mosaic_file), 'read') ) &
 127        call mpp_error(FATAL, 'test_mosaic: error in opening file '//'INPUT/'//trim(ocn_mosaic_file))
 128
 129   answers(1,:) = (/1440, 1440, 1, 1080, 1, 1, 1, 1080 /)
 130   answers(2,:) = (/1, 720, 1080, 1080, 1440, 721, 1080, 1080 /)
 131
 132   ntiles = get_mosaic_ntiles(ocn_fileobj)
 133   ncontacts = get_mosaic_ncontacts(ocn_fileobj)
 134
 135   allocate(nx(ntiles), ny(ntiles))
 136   allocate(tile1(ncontacts), tile2(ncontacts) )
 137   allocate(istart1(ncontacts), iend1(ncontacts), jstart1(ncontacts), jend1(ncontacts) )
 138   allocate(istart2(ncontacts), iend2(ncontacts), jstart2(ncontacts), jend2(ncontacts) )
 139
 140   call get_mosaic_contact(ocn_fileobj, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2)
 141
 142   !< Compare with expected results:
 143   if (ntiles .ne. 1)    call mpp_error(FATAL, "ntiles is not equal to 1")
 144   if (ncontacts .ne. 2) call mpp_error(FATAL, "ncontacts is not the expected result")
 145   do n = 1, ncontacts
 146      if (istart1(n) .ne. answers(n,1)) call mpp_error(FATAL, "istart1 is not the expected result")
 147      if (iend1(n)   .ne. answers(n,2)) call mpp_error(FATAL, "iend1 is not the expected result")
 148
 149      if (jstart1(n) .ne. answers(n,3)) call mpp_error(FATAL, "jstart1 is not the expected result")
 150      if (jend1(n)   .ne. answers(n,4)) call mpp_error(FATAL, "jend1 is not the expected result")
 151
 152      if (istart2(n) .ne. answers(n,5)) call mpp_error(FATAL, "istart2 is not the expected result")
 153      if (iend2(n)   .ne. answers(n,6)) call mpp_error(FATAL, "iend2 is not the expected result")
 154
 155      if (jstart2(n) .ne. answers(n,7)) call mpp_error(FATAL, "jstart2 is not the expected result")
 156      if (jend2(n)   .ne. answers(n,8)) call mpp_error(FATAL, "jend2 is not the expected result")
 157   end do
 158
 159   deallocate(tile1, tile2, nx, ny)
 160   deallocate(istart1, iend1, jstart1, jend1)
 161   deallocate(istart2, iend2, jstart2, jend2)
 162
 163 end subroutine test_get_mosaic_contact
 164 !------------------------------------------------------!
 165 subroutine test_calc_mosaic_grid_area
 166
 167   !> This subroutine tests get_grid_area
 168
 169   implicit none
 170
 171   real(TEST_MOS_KIND_) :: x_rad(c1_nx, c1_ny), y_rad(c1_nx, c1_ny) !< x and y in radians
 172   real(TEST_MOS_KIND_) :: area_out(1,1) !< area to be computed
 173
 174   !> x_rad and y_rad can be set to be be the entire cell
 175   !! x_rad = x(1:3:2, 1:3:2) and y_rad = y(1:3:2, 1:3:2)
 176   !! The answer will then be 4.0*area(1,1)
 177   x_rad = real( real(x(1:2,1:2),r8_kind)*DEG_TO_RAD,TEST_MOS_KIND_) !< set coordinates
 178   y_rad = real( real(y(1:2,1:2),r8_kind)*DEG_TO_RAD,TEST_MOS_KIND_) !< set coordinates
 179
 180   call calc_mosaic_grid_area(x_rad, y_rad, area_out)
 181   call check_answer(area(1,1), area_out(1,1), 'TEST_CALC_MOSAIC_GRID_AREA')
 182
 183 end subroutine test_calc_mosaic_grid_area
 184 !------------------------------------------------------!
 185 subroutine test_get_grid_great_circle_area
 186
 187   !> This subroutine tests calc_mosaic_grid_great_circle_area
 188
 189   implicit none
 190
 191   real(TEST_MOS_KIND_) :: x_rad(c1_nx, c1_ny), y_rad(c1_nx, c1_ny) !< x and y in radians
 192   real(TEST_MOS_KIND_) :: area_out(1,1) !< area to be computed
 193
 194   !> x_rad and y_rad can be set to be be the entire cell
 195   !! x_rad = x(1:3:2, 1:3:2) and y_rad = y(1:3:2, 1:3:2)
 196   !! The answer will then be 4.0*area(1,1)
 197   x_rad = real( real(x(1:2,1:2),r8_kind)*DEG_TO_RAD,TEST_MOS_KIND_) !< set coordinates
 198   y_rad = real( real(y(1:2,1:2),r8_kind)*DEG_TO_RAD,TEST_MOS_KIND_) !< set coordinates
 199   call calc_mosaic_grid_great_circle_area(x_rad, y_rad, area_out)
 200   call check_answer(area(1,1), area_out(1,1), 'TEST_GET_GRID_GREAT_CIRCLE_AREA')
 201
 202 end subroutine test_get_grid_great_circle_area
 203 !------------------------------------------------------!
 204 subroutine test_get_mosaic_xgrid
 205
 206   !> Test get_mosaic_xgrid
 207
 208   implicit none
 209
 210   integer, dimension(ncells) :: i1, j1, i2, j2 !< indices of parent cells
 211   real(TEST_MOS_KIND_), dimension(ncells) :: area !< area to be returned
 212   real(r8_kind) :: garea, get_global_area !< global area
 213   integer :: i !< counter
 214
 215   type(FmsNetcdfFile_t) x_fileobj
 216
 217   garea = get_global_area()
 218
 219   if( .not. open_file(x_fileobj, 'INPUT/'//trim(exchange_file), 'read')) &
 220        call mpp_error(FATAL, 'test_mosaic: error in opening file '//'INPUT/'//trim(exchange_file))
 221
 222   call get_mosaic_xgrid(x_fileobj, i1, j1, i2, j2, area)
 223
 224   !> check answers
 225   do i=1, ncells
 226      call check_answer( real(real(xgrid_area(i),r8_kind)/garea,lkind), area(i),"TEST_GET_MOSAIC_XGRID area")
 227      call check_answer(tile1_cell(1,i), i1(i), "TEST_GET_MOSAIC_XGRID i1")
 228      call check_answer(tile2_cell(1,i), i2(i), "TEST_GET_MOSAIC_XGRID i2")
 229      call check_answer(tile1_cell(1,i), j1(i), "TEST_GET_MOSAIC_XGRID j1")
 230      call check_answer(tile2_cell(1,i), j2(i), "TEST_GET_MOSAIC_XGRID j2")
 231   end do
 232
 233   call close_file(x_fileobj)
 234
 235
 236 end subroutine test_get_mosaic_xgrid
 237 !------------------------------------------------------!
 238 subroutine test_is_inside_polygon
 239
 240   !> cheating a little.  starting with xyz coordinates (cause easier to understand)
 241
 242   implicit none
 243
 244   integer, parameter :: n=5
 245   integer :: i
 246   real(TEST_MOS_KIND_) :: lat1, lon1, x1, y1, z1, r
 247   real(TEST_MOS_KIND_), dimension(n) :: lon2, lat2, x2, y2, z2
 248   logical :: answer, is_inside
 249
 250   integer, parameter :: lkind=TEST_MOS_KIND_ !< local kind
 251
 252   !> polygon
 253   x2=0.0_lkind
 254   y2(1)=1.0_lkind ; y2(2)=1.0_lkind ; y2(3)=4.0_lkind ; y2(4)=4.0_lkind ; y2(5)=1.0_lkind
 255   z2(1)=2.0_lkind ; z2(2)=4.0_lkind ; z2(3)=4.0_lkind ; z2(4)=2.0_lkind ; z2(5)=2.0_lkind
 256   do i=1, n
 257      r = sqrt( x2(i)**2 + y2(i)**2 + z2(i)**2 )
 258      lon2(i)=atan2(y2(i), x2(i))
 259      lat2(i)=asin(z2(i)/r)
 260   end do
 261
 262   !> point outside of the polygon
 263   x1=2.0_lkind
 264   y1=5.0_lkind
 265   z1=4.2_lkind
 266   r = sqrt(x1**2+y1**2+z1**2)
 267   lon1=atan2(y1, x1)
 268   lat1=asin(z1/r)
 269
 270   answer=.false.
 271   is_inside=is_inside_polygon(lon1, lat1, lon2, lat2)
 272   call check_answer(answer,is_inside,' TEST_IS_INSIDE_POLYGON')
 273
 274   !> point inside the polygon
 275   x1=0.0_lkind
 276   y1=3.0_lkind
 277   z1=2.5_lkind
 278   r = sqrt(x1**2+y1**2+z1**2)
 279   lon1=atan2(y1, x1)
 280   lat1=asin(z1/r)
 281
 282   answer=.true.
 283   is_inside=is_inside_polygon(lon1, lat1, lon2, lat2)
 284   call check_answer(answer,is_inside,'TEST_IS_INSIDE_POLYGON')
 285
 286 end subroutine test_is_inside_polygon
 287 !------------------------------------------------------!
 288 subroutine check_answer(answer, myvalue, whoami)
 289
 290   implicit none
 291   class(*) :: answer
 292   class(*) :: myvalue
 293   character(*) :: whoami
 294
 295   select type(answer)
 296   type is ( logical )
 297      select type(myvalue)
 298      type is( logical )
 299         if( answer .neqv. myvalue ) then
 300            write(*,*) '*************************************'
 301            write(*,*) 'EXPECTED ', answer, 'but got ', myvalue
 302            call mpp_error( FATAL,'failed '//trim(whoami) )
 303         end if
 304      end select
 305   type is( real(r4_kind) )
 306      select type( myvalue)
 307         type is(real(r4_kind) )
 308            if( answer .ne. myvalue ) then
 309               write(*,*) '*************************************'
 310               write(*,*) 'EXPECTED ', answer, 'but got ', myvalue
 311               write(*,*) 'difference of', abs(answer-myvalue)
 312               call mpp_error( FATAL,'failed '//trim(whoami) )
 313            end if
 314         end select
 315   type is( real(r8_kind) )
 316      select type( myvalue)
 317         type is(real(r4_kind) )
 318            if( answer .ne. myvalue ) then
 319               write(*,*) '*************************************'
 320               write(*,*) 'EXPECTED ', answer, 'but got ', myvalue
 321               write(*,*) 'difference of', abs(answer-myvalue)
 322               call mpp_error( FATAL,'failed '//trim(whoami) )
 323            end if
 324         end select
 325   end select
 326
 327 end subroutine check_answer
 328 !------------------------------------------------------!
 329 end program test_mosaic