test_fms/diag_manager/check_time_pow.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  Checks the output file after running test_reduction_methods using the "time_pow" reduction method
  21 !! Pow reductions are run with a different dataset to simplify checking
  22 !! Each element in sent arrays is just the sum of its indices
  23 program check_time_pow
  24   use fms_mod,           only: fms_init, fms_end, string
  25   use fms2_io_mod,       only: FmsNetcdfFile_t, read_data, close_file, open_file
  26   use mpp_mod,           only: mpp_npes, mpp_error, FATAL, mpp_pe, input_nml_file, NOTE
  27   use platform_mod,      only: r4_kind, r8_kind
  28   use testing_utils,     only: allocate_buffer, test_normal, test_openmp, test_halos, no_mask, logical_mask, real_mask
  29
  30   implicit none
  31
  32   type(FmsNetcdfFile_t)              :: fileobj            !< FMS2 fileobj
  33   type(FmsNetcdfFile_t)              :: fileobj1           !< FMS2 fileobj for subregional file 1
  34   type(FmsNetcdfFile_t)              :: fileobj2           !< FMS2 fileobj for subregional file 2
  35   real(kind=r4_kind), allocatable    :: cdata_out(:,:,:,:) !< Data in the compute domain
  36   integer                            :: nx                 !< Number of points in the x direction
  37   integer                            :: ny                 !< Number of points in the y direction
  38   integer                            :: nz                 !< Number of points in the z direction
  39   integer                            :: nw                 !< Number of points in the 4th dimension
  40   integer                            :: ti                  !< For looping through time levels
  41   integer                            :: io_status          !< Io status after reading the namelist
  42   logical                            :: use_mask           !< .true. if using masks
  43   integer, parameter :: file_freq = 6 !< file frequency as set in diag_table.yaml
  44   integer, parameter :: pow_value = 2 !< pow value as set in reduction method (ie. pow2)
  45
  46   integer :: test_case = test_normal !< Indicates which test case to run
  47   integer :: mask_case = no_mask     !< Indicates which masking option to run
  48   integer, parameter :: kindl = KIND(0.0) !< compile-time default kind size
  49
  50   namelist / test_reduction_methods_nml / test_case, mask_case
  51
  52   call fms_init()
  53
  54   read (input_nml_file, test_reduction_methods_nml, iostat=io_status)
  55
  56   select case(mask_case)
  57   case (no_mask)
  58     use_mask = .false.
  59   case (logical_mask, real_mask)
  60     use_mask = .true.
  61   end select
  62   nx = 96
  63   ny = 96
  64   nz = 5
  65   nw = 2
  66
  67   if (.not. open_file(fileobj, "test_pow.nc", "read")) &
  68     call mpp_error(FATAL, "unable to open test_pow.nc")
  69
  70   if (.not. open_file(fileobj1, "test_pow_regional.nc.0004", "read")) &
  71     call mpp_error(FATAL, "unable to open test_pow_regional.nc.0004")
  72
  73   if (.not. open_file(fileobj2, "test_pow_regional.nc.0005", "read")) &
  74     call mpp_error(FATAL, "unable to open test_pow_regional.nc.0005")
  75
  76   cdata_out = allocate_buffer(1, nx, 1, ny, nz, nw)
  77
  78   do ti = 1, 8
  79     cdata_out = -999_r4_kind
  80     print *, "Checking answers for var0_pow - time_level:", string(ti)
  81     call read_data(fileobj, "var0_pow", cdata_out(1,1,1,1), unlim_dim_level=ti)
  82     call check_data_0d(cdata_out(1,1,1,1), ti)
  83
  84     cdata_out = -999_r4_kind
  85     print *, "Checking answers for var1_pow - time_level:", string(ti)
  86     call read_data(fileobj, "var1_pow", cdata_out(:,1,1,1), unlim_dim_level=ti)
  87     call check_data_1d(cdata_out(:,1,1,1), ti)
  88
  89     cdata_out = -999_r4_kind
  90     print *, "Checking answers for var2_pow - time_level:", string(ti)
  91     call read_data(fileobj, "var2_pow", cdata_out(:,:,1,1), unlim_dim_level=ti)
  92     call check_data_2d(cdata_out(:,:,1,1), ti)
  93
  94     cdata_out = -999_r4_kind
  95     print *, "Checking answers for var3_pow - time_level:", string(ti)
  96     call read_data(fileobj, "var3_pow", cdata_out(:,:,:,1), unlim_dim_level=ti)
  97     call check_data_3d(cdata_out(:,:,:,1), ti, .false.)
  98
  99     cdata_out = -999_r4_kind
 100     print *, "Checking answers for var3_Z - time_level:", string(ti)
 101     call read_data(fileobj, "var3_Z", cdata_out(:,:,1:2,1), unlim_dim_level=ti)
 102     call check_data_3d(cdata_out(:,:,1:2,1), ti, .true., nz_offset=1)
 103
 104     cdata_out = -999_r4_kind
 105     print *, "Checking answers for var3_pow in the first regional file- time_level:", string(ti)
 106     call read_data(fileobj1, "var3_pow", cdata_out(1:4,1:3,1:2,1), unlim_dim_level=ti)
 107     call check_data_3d(cdata_out(1:4,1:3,1:2,1), ti, .true., nx_offset=77, ny_offset=77, nz_offset=1)
 108
 109     cdata_out = -999_r4_kind
 110     print *, "Checking answers for var3_pow in the second regional file- time_level:", string(ti)
 111     call read_data(fileobj2, "var3_pow", cdata_out(1:4,1:1,1:2,1), unlim_dim_level=ti)
 112     call check_data_3d(cdata_out(1:4,1:1,1:2,1), ti, .true., nx_offset=77, ny_offset=80, nz_offset=1)
 113   enddo
 114
 115   call fms_end()
 116
 117 contains
 118
 119   ! sent data set to:
 120   !          buffer(ii-is+1+nhalo, j-js+1+nhalo, k, l) = i + j + k + l
 121   !                                                      + time_index/100
 122   ! sum of squares for 1..n can be calculated with:
 123   !  P(n) = (n^3 / 3) + (n^2 / 2) + (n/6)
 124   !> @brief Check that the 0d data read in is correct
 125   subroutine check_data_0d(buffer, time_level)
 126     real(kind=r4_kind), intent(inout) :: buffer        !< Buffer read from the table
 127     integer,            intent(in)    :: time_level    !< Time level read in
 128
 129     real(kind=r4_kind)                :: buffer_exp    !< Expected result
 130     integer :: i, step_pow = 0 !< pow of time step increments to use in generating reference data
 131
 132     ! only one index(1,1,1,1) = sums to 4
 133     buffer_exp = get_answer_from_index(4)
 134
 135     if (abs(buffer - buffer_exp) > 0.0) then
 136       print *, "time_level", time_level, "expected", buffer_exp, "read", buffer
 137       call mpp_error(FATAL, "Check_time_pow::check_data_0d:: Data is not correct")
 138     endif
 139   end subroutine check_data_0d
 140
 141   !> @brief Check that the 1d data read in is correct
 142   subroutine check_data_1d(buffer, time_level)
 143     real(kind=r4_kind), intent(in)    :: buffer(:)     !< Buffer read from the table
 144     integer,            intent(in)    :: time_level    !< Time level read in
 145     real(kind=r4_kind)                :: buffer_exp    !< Expected result
 146     integer :: step_sum !< pow of time step increments to use in generating reference data
 147     integer :: ii, i, j, k, l !< For looping
 148     integer :: n
 149
 150     ! 1d answer is
 151     !  (((i * 1000 + 11) * frequency) + (sum of time steps)) / frequency
 152     !  or
 153     !  => (i * 1000 + 11) + (sum of time_steps/frequency/100)
 154     do ii = 1, size(buffer, 1)
 155       buffer_exp = get_answer_from_index(ii + 3)
 156
 157       if (use_mask .and. ii .eq. 1) buffer_exp = -666_r4_kind
 158       if (abs(buffer(ii) - buffer_exp) > 0.0) then
 159         print *, "i:", ii, "read in:", buffer(ii), "expected:", buffer_exp, "time level:", time_level
 160         print *, "diff:", abs(buffer(ii) - buffer_exp)
 161         call mpp_error(FATAL, "Check_time_pow::check_data_1d:: Data is not exact")
 162       endif
 163     enddo
 164   end subroutine check_data_1d
 165
 166   !> @brief Check that the 2d data read in is correct
 167   subroutine check_data_2d(buffer, time_level)
 168     real(kind=r4_kind), intent(in)    :: buffer(:,:)   !< Buffer read from the table
 169     integer,            intent(in)    :: time_level    !< Time level read in
 170     real(kind=r4_kind)                :: buffer_exp    !< Expected result
 171
 172     integer :: ii,i, j, k, l !< For looping
 173     integer :: step_pow !< pow of time step increments to use in generating reference data
 174
 175     do ii = 1, size(buffer, 1)
 176       do j = 1, size(buffer, 2)
 177         buffer_exp = get_answer_from_index(ii + j + 2)
 178         if (use_mask .and. ii .eq. 1 .and. j .eq. 1) buffer_exp = -666_r4_kind
 179         if (abs(buffer(ii, j) - buffer_exp) > 0.0) then
 180           print *, "indices:", ii, j, "expected:", buffer_exp, "read in:",buffer(ii, j)
 181           call mpp_error(FATAL, "Check_time_pow::check_data_2d:: Data is not correct")
 182         endif
 183       enddo
 184     enddo
 185   end subroutine check_data_2d
 186
 187   !> @brief Check that the 3d data read in is correct
 188   subroutine check_data_3d(buffer, time_level, is_regional, nx_offset, ny_offset, nz_offset)
 189     real(kind=r4_kind), intent(in)    :: buffer(:,:,:) !< Buffer read from the table
 190     integer,            intent(in)    :: time_level    !< Time level read in
 191     logical,            intent(in)    :: is_regional   !< .True. if the variable is subregional
 192     real(kind=r4_kind)                :: buffer_exp    !< Expected result
 193     integer, optional,  intent(in)    :: nx_offset     !< Offset in the x direction
 194     integer, optional,  intent(in)    :: ny_offset     !< Offset in the y direction
 195     integer, optional,  intent(in)    :: nz_offset     !< Offset in the z direction
 196
 197     integer :: ii, i, j, k, l !< For looping
 198     integer :: nx_oset !< Offset in the x direction (local variable)
 199     integer :: ny_oset !< Offset in the y direction (local variable)
 200     integer :: nz_oset !< Offset in the z direction (local variable)
 201     integer :: step_pow!< pow of time step increments to use in generating reference data
 202
 203     step_pow = 0
 204     do i=(time_level-1)*file_freq+1, time_level*file_freq
 205       step_pow = step_pow + i
 206     enddo
 207
 208     nx_oset = 0
 209     if (present(nx_offset)) nx_oset = nx_offset
 210
 211     ny_oset = 0
 212     if (present(ny_offset)) ny_oset = ny_offset
 213
 214     nz_oset = 0
 215     if (present(nz_offset)) nz_oset = nz_offset
 216
 217     ! 3d answer is
 218     !  ((i * 1000 + j * 10 + k) * frequency) + (pow of time steps)
 219     do ii = 1, size(buffer, 1)
 220       do j = 1, size(buffer, 2)
 221         do k = 1, size(buffer, 3)
 222           buffer_exp = get_answer_from_index(ii + j + k + 1 + ny_oset + nx_oset + nz_oset)
 223           if (use_mask .and. ii .eq. 1 .and. j .eq. 1 .and. k .eq. 1 .and. .not. is_regional) buffer_exp = -666_r4_kind
 224           if (abs(buffer(ii, j, k) - buffer_exp) > 0.0) then
 225             print *, mpp_pe(),'indices:',ii, j, k, "read in:", buffer(ii, j, k), "expected:",buffer_exp
 226             call mpp_error(FATAL, "Check_time_pow::check_data_3d:: Data is not correct")
 227           endif
 228         enddo
 229       enddo
 230     enddo
 231   end subroutine check_data_3d
 232
 233   function get_answer_from_index(index_sum) &
 234     result(answ)
 235     integer, intent(in) :: index_sum !< sum of indices
 236     real(r4_kind) :: answ
 237     integer :: i
 238     answ = 0
 239     do i=1, file_freq
 240       answ = answ + real(index_sum, r4_kind) ** 2.0
 241     enddo
 242     answ = answ / file_freq
 243   end function
 244
 245 end program