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