test_fms/mosaic2/test_grid2.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 !> @brief  This programs tests calls to get_mosaic_ntiles, get_mosaic_ncontacts,
  20 !! get_mosaic_grid_sizes, get_mosaic_contact.  All subroutines here are tested
  21 !! with C1 tiles where tiles 1-6 are identical.  The tile points are made up with
  22 !! values that result in simple answers.  See write_files module 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 mpp_mod,      only : mpp_init, mpp_error, FATAL, mpp_npes, mpp_pe, mpp_root_pe
  29 use mpp_domains_mod, only: domain2D, domainUG, mpp_define_domains, mpp_get_compute_domain, mpp_define_unstruct_domain
  30 use fms2_io_mod,  only : open_file, close_file, FmsNetcdfFile_t, fms2_io_init
  31 use fms2_io_mod,  only : register_axis, register_field, write_data, read_data
  32 use fms_mod,      only : fms_init, fms_end
  33 use platform_mod, only : r4_kind, r8_kind
  34 use grid2_mod
  35 use write_files
  36
  37 implicit none
  38
  39 !> write out netcdf files
  40 !! write_all sets up the grids
  41 call fms2_io_init()
  42 call write_all()
  43 call fms_init()
  44
  45 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_CELL_VERTICIES'
  46 call test_get_cell_vertices
  47
  48 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_CELL_CENTERS'
  49 call test_get_cell_centers
  50
  51 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_GRID_CELL_AREA_SG'
  52 call test_get_grid_cell_area_sg
  53
  54 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST_GET_GRID_CELL_AREA_UG'
  55 call test_get_grid_cell_area_ug
  56
  57 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_GRID_COMP_AREA_SG'
  58 call test_get_grid_comp_area_sg
  59
  60 if(mpp_pe() .eq. mpp_root_pe()) write(*,*) 'TEST GET_GRID_COMP_AREA_UG'
  61 call test_get_grid_comp_area_ug
  62
  63 contains
  64   !------------------------------------------!
  65   subroutine test_get_cell_vertices
  66
  67     !> This subroutine tests get_cell_verticees. This
  68     !! subroutine only tests for vertices in tile 1.
  69
  70     implicit none
  71
  72     real(TEST_MOS_KIND_) :: lonb_2d(c1_nx,c1_ny)       !< returned values for lon 2d
  73     real(TEST_MOS_KIND_) :: latb_2d(c1_nx,c1_ny)       !< returned values for lat 2d
  74     real(TEST_MOS_KIND_) :: answer_lon_2d(c1_nx,c1_ny) !< answers for lon 2d
  75     real(TEST_MOS_KIND_) :: answer_lat_2d(c1_nx,c1_ny) !< answers for lat 2d
  76
  77     integer :: i,j
  78
  79     !> answers
  80     answer_lon_2d=x(1:c1_nxp:2, 1:c1_nxp:2)
  81     answer_lat_2d=y(1:c1_nxp:2, 1:c1_nxp:2)
  82
  83     call get_grid_cell_vertices('ATM',1,lonb_2d,latb_2d)
  84     !> check
  85     do j=1, c1_ny
  86        do i=1, c1_nx
  87           call check_answer(answer_lon_2d(i,j), lonb_2d(i,j), 'TEST_GET_CELL_VERTICIES_2D lon')
  88           call check_answer(answer_lat_2d(i,j), latb_2d(i,j), 'TEST_GET_CELL_VERTICIES_2D lat')
  89        end do
  90     end do
  91
  92   end subroutine test_get_cell_vertices
  93   !------------------------------------------!
  94   subroutine test_get_cell_centers
  95
  96     !> This subroutine tests get_cell_centers.
  97     !! There is only one cell center point in a C1 tile.
  98
  99     implicit none
 100
 101     integer, parameter :: nx = c1_nx/2 !< number of center points
 102     integer, parameter :: ny = c1_ny/2 !< number of center points
 103
 104     real(TEST_MOS_KIND_) :: glon_2d(nx,ny) !< results from grid_cell_centers
 105     real(TEST_MOS_KIND_) :: glat_2d(nx,ny) !< results from grid_cell_centers
 106     real(TEST_MOS_KIND_) :: answer_glon_2d(nx,ny) !< answers for glon
 107     real(TEST_MOS_KIND_) :: answer_glat_2d(nx,ny) !< answers for glat
 108
 109     integer :: i, j
 110
 111     !--- 2d ---!
 112     answer_glon_2d=x(2:c1_nx:2, 2:c1_nx:2)
 113     answer_glat_2d=y(2:c1_nx:2, 2:c1_nx:2)
 114
 115     call get_grid_cell_centers('ATM', 1, glon_2d, glat_2d)
 116     do i=1, nx
 117        do j=1, ny
 118           call check_answer(answer_glon_2d(j,i), glon_2d(j,i), 'TEST_GET_CELL_CENTERS_2D lon')
 119           call check_answer(answer_glat_2d(j,i), glat_2d(j,i), 'TEST_GET_CELL_CENTERS_2D lat')
 120        end do
 121     end do
 122
 123   end subroutine test_get_cell_centers
 124   !------------------------------------------!
 125   subroutine test_get_grid_cell_area_sg
 126
 127     !> This subroutine tests get_grid_cell_area_SG
 128     !! first without the domain input argument and second
 129     !! with the domain input argument
 130
 131     implicit none
 132
 133     type(domain2D) :: SG_domain
 134     real(TEST_MOS_KIND_) :: area_out2(1,1)
 135     real(TEST_MOS_KIND_) :: answer
 136
 137     answer = real(2.0_r8_kind*PI*RADIUS*RADIUS,lkind)
 138
 139     !> total of 1 domain with 1 (center) point in the domain
 140     call mpp_define_domains((/1,1,1,1/), (/1,1/), SG_domain)
 141
 142     !> The area computed by get_grid_cell_area is for the entire cell
 143     !! The array area, set in write_files.F90, is the area for 1/4th of the cell
 144
 145     !> Test withtout SG_domain
 146     call get_grid_cell_area('ATM',2, area_out2)
 147     call check_answer(answer, area_out2(1,1), 'TEST_GRID_CELL_AREA_SG')
 148
 149     !> Test with SG_domain
 150     call get_grid_cell_area('ATM',2, area_out2, SG_domain)
 151     call check_answer(answer, area_out2(1,1), 'TEST_GRID_CELL_AREA_SG with SG_domain')
 152
 153   end subroutine test_get_grid_cell_area_sg
 154   !------------------------------------------!
 155   subroutine test_get_grid_cell_area_ug
 156
 157     !> This subroutine tests get_grid_cell_area_ug
 158
 159     implicit none
 160     type(domain2D) :: SG_domain
 161     type(domainUG) :: UG_domain !< UG_domain is the same as SG_domain
 162     real(TEST_MOS_KIND_) :: area_out1(1)
 163     real(TEST_MOS_KIND_) :: answer
 164     integer :: i
 165     integer :: npts_tile(1),grid_nlevel(1), ndivs, grid_index(1)
 166
 167     npts_tile=1
 168     grid_nlevel=1
 169     ndivs=1
 170     grid_index=1
 171
 172     answer = real( 4.0_r8_kind * area(1,1), TEST_MOS_KIND_)
 173
 174     !> The unstructured grid is the same as the structured grid; there's only one center point in the tile.
 175     call mpp_define_domains((/1,1,1,1/), (/1,1/), SG_domain)
 176     call mpp_define_unstruct_domain(UG_domain, SG_domain,npts_tile,grid_nlevel,&
 177                                     mpp_npes(),ndivs,grid_index,name='immadeup')
 178
 179     !> The area computed by get_grid_cell_area is for the entire cell
 180     !! The array area, set in write_files.F90, is the area for 1/4th of the cell
 181     call get_grid_cell_area('ATM',1, area_out1, SG_domain, UG_domain)
 182     call check_answer(answer, area_out1(1), 'TEST_GRID_CELL_AREA_UG')
 183
 184   end subroutine test_get_grid_cell_area_ug
 185   !------------------------------------------!
 186   subroutine test_get_grid_comp_area_sg
 187
 188     !> This subroutine tests get_grid_comp_area_sg
 189     !! first without the domain input argument and second
 190     !! with the domain input argument
 191
 192     implicit none
 193     type(domain2D) :: SG_domain
 194     real(TEST_MOS_KIND_) :: area_out2(1,1)
 195     real(TEST_MOS_KIND_) :: answer
 196
 197     answer = real( 4.0_r8_kind * area(1,1), TEST_MOS_KIND_)
 198
 199     call mpp_define_domains((/1,1,1,1/), (/1,1/), SG_domain)
 200
 201     !> The area computed by get_grid_cell_area is for the entire cell
 202     !! The array area, set in write_files.F90, is the area for 1/4th of the cell
 203     !! Test without SG_domain
 204     call get_grid_comp_area('ATM', 1, area_out2)
 205     call check_answer(answer, area_out2(1,1), 'TEST_GRID_COMP_AREA_SG')
 206
 207     !> The area computed by get_grid_cell_area is for the entire cell
 208     !! The array area, set in write_files.F90, is the area for 1/4th of the cell
 209     !! Test with SG_domain
 210     call get_grid_comp_area('ATM', 1, area_out2, SG_domain)
 211     call check_answer(answer, area_out2(1,1), 'TEST_GRID_COMP_AREA_SG with SG_domain')
 212
 213   end subroutine test_get_grid_comp_area_sg
 214   !------------------------------------------!
 215   subroutine test_get_grid_comp_area_ug
 216
 217     !> This subroutine tests get_grid_comp_area_ug
 218
 219     implicit none
 220     type(domain2D) :: SG_domain
 221     type(domainUG) :: UG_domain !< UG_domain is the same as SG_domain
 222     integer :: npts_tile(1), ntiles_grid(1), grid_index(1)
 223     real(TEST_MOS_KIND_) :: answer
 224     real(TEST_MOS_KIND_) :: area_out1(1)
 225
 226     npts_tile=1
 227     ntiles_grid=1
 228     grid_index(1)=1
 229     answer = real( 4.0_r8_kind * area(1,1), TEST_MOS_KIND_)
 230
 231     !> the unstructured grid is the same as the structured grid
 232     call mpp_define_domains((/1,1,1,1/), (/1,1/), SG_domain)
 233     call mpp_define_unstruct_domain(UG_domain, SG_domain, npts_tile, ntiles_grid,mpp_npes(),1,grid_index)
 234
 235     !> The area computed by get_grid_cell_area is for the entire cell
 236     !! The array area, set in write_files.F90, is the area for 1/4th of the cell
 237     call get_grid_comp_area('ATM',3,area_out1,SG_domain, UG_domain)
 238     call check_answer(answer, area_out1(1), 'TEST_GRID_CELL_AREA_UG')
 239
 240   end subroutine test_get_grid_comp_area_ug
 241   !------------------------------------------!
 242   subroutine check_answer(answer, myvalue, whoami)
 243
 244     implicit none
 245     real(TEST_MOS_KIND_) :: answer
 246     real(TEST_MOS_KIND_) :: myvalue
 247     character(*) :: whoami
 248
 249     if( answer .ne. myvalue ) then
 250        write(*,*) '*************************************'
 251        write(*,*) 'EXPECTED ', answer, 'but got ', myvalue
 252        write(*,*) 'difference of', abs(answer-myvalue)
 253        call mpp_error( FATAL,'failed '//trim(whoami) )
 254     end if
 255
 256   end subroutine check_answer
 257 !------------------------------------------------------!
 258 end program test_mosaic