test_fms/diag_manager/check_time_diurnal.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 !! TODO more complicated cases and data
  20 !> @brief  Checks the output file after running test_reduction_methods using the "time_diurnal" reduction method
  21 program check_time_diurnal
  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, NOTE
  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   real(kind=r4_kind), allocatable    :: cdata_out_5d(:,:,:,:,:) !< Data in the compute domain
  35   integer                            :: nx                 !< Number of points in the x direction
  36   integer                            :: ny                 !< Number of points in the y direction
  37   integer                            :: nz                 !< Number of points in the z direction
  38   integer                            :: nw                 !< Number of points in the w direction
  39   integer                            :: nd                 !< Number of points in the diurnal axis
  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
  45   integer :: test_case = test_normal !< Indicates which test case to run
  46   integer :: mask_case = no_mask     !< Indicates which masking option to run
  47   integer, parameter :: kindl = KIND(0.0) !< compile-time default kind size
  48   integer :: nmonths !< number of months the test ran for
  49   namelist / test_diag_diurnal_nml / test_case, mask_case
  50
  51   call fms_init()
  52
  53   read (input_nml_file, test_diag_diurnal_nml, iostat=io_status)
  54
  55   select case(mask_case)
  56   case (no_mask)
  57     use_mask = .false.
  58   case (logical_mask, real_mask)
  59     use_mask = .true.
  60   end select
  61   nx = 96
  62   ny = 96
  63   nz = 5
  64   nw = 2
  65   nmonths = 3
  66   nd = 3 !< diurnal sample size
  67
  68   if (.not. open_file(fileobj, "test_diurnal.nc", "read")) &
  69     call mpp_error(FATAL, "unable to open test_diurnal.nc")
  70
  71   if (.not. open_file(fileobj1, "test_diurnal_regional.nc.0004", "read")) &
  72     call mpp_error(FATAL, "unable to open test_diurnal_regional.nc.0004")
  73
  74   if (.not. open_file(fileobj2, "test_diurnal_regional.nc.0005", "read")) &
  75     call mpp_error(FATAL, "unable to open test_diurnal_regional.nc.0005")
  76
  77   !cdata_out = allocate_buffer(1, nx, 1, ny, nz, nd)
  78   allocate(cdata_out(nx, ny, nz, nd))
  79   allocate(cdata_out_5d(nx, ny, nz, nw, nd))
  80
  81   do ti = 1, nmonths
  82     cdata_out = -999_r4_kind
  83     print *, "Checking answers for var1 - time_level:", string(ti)
  84     call read_data(fileobj, "var1", cdata_out(:,1:nd,1,1), unlim_dim_level=ti)
  85     call check_data_1d(cdata_out(:,1:nd,1,1), ti, sample_size=nd)
  86
  87     cdata_out = -999_r4_kind
  88     print *, "Checking answers for var2 - time_level:", string(ti)
  89     call read_data(fileobj, "var2", cdata_out(:,:,1:nd,1), unlim_dim_level=ti)
  90     call check_data_2d(cdata_out(:,:,1:nd,1), ti, sample_size=nd)
  91
  92     cdata_out = -999_r4_kind
  93     print *, "Checking answers for var3 - time_level:", string(ti)
  94     call read_data(fileobj, "var3", cdata_out(:,:,:,:), unlim_dim_level=ti)
  95     call check_data_3d(cdata_out(:,:,:,:), ti, .false., sample_size=nd)
  96
  97     cdata_out = -999_r4_kind
  98     print *, "Checking answers for var4_diurnal - time_level:", string(ti)
  99     call read_data(fileobj, "var4", cdata_out_5d, unlim_dim_level=ti)
 100     call check_data_4d(cdata_out_5d(:,:,:,:,:), ti, .false., sample_size=nd)
 101
 102     cdata_out = -999_r4_kind
 103     print *, "Checking answers for var3_diurnal in the first regional file- time_level:", string(ti)
 104     call read_data(fileobj1, "var3_diurnal", cdata_out(1:4,1:3,1:2,1:1), unlim_dim_level=ti)
 105     call check_data_3d(cdata_out(1:4,1:3,1:2,1:1), ti, .true., sample_size=nd, nx_offset=77, ny_offset=77, nz_offset=1)
 106
 107     cdata_out = -999_r4_kind
 108     print *, "Checking answers for var3_diurnal in the second regional file- time_level:", string(ti)
 109     call read_data(fileobj2, "var3_diurnal", cdata_out(1:4,1:1,1:2,1:1), unlim_dim_level=ti)
 110     call check_data_3d(cdata_out(1:4,1:1,1:2,1:1), ti, .true., sample_size=nd, nx_offset=77, ny_offset=80, nz_offset=1)
 111   enddo
 112
 113   call fms_end()
 114
 115 contains
 116
 117
 118   !> @brief Check that the 2d data read in is correct
 119   subroutine check_data_2d(buffer, time_level, sample_size)
 120     real(kind=r4_kind), intent(in)    :: buffer(:,:,:)   !< Buffer read from the table (2d + the diurnal axis)
 121     integer,            intent(in)    :: time_level    !< Time level read in
 122     integer,            intent(in)    :: sample_size   !< diurnal sample size of variable to check
 123     real(kind=r4_kind)                :: buffer_exp    !< Expected result
 124
 125     integer :: ii,i, j, k, l, d!< For looping
 126     integer :: step_avg !< avg of time step increments to use in generating reference data
 127     integer :: d_index
 128     real(r8_kind) :: hrly_sums(sample_size)
 129
 130     ! sum of hours in diurnal section
 131     hrly_sums = 0
 132     do i=1, 23
 133       d_index = i / (24/sample_size) + 1
 134       hrly_sums(d_index) = hrly_sums(d_index) + i
 135     enddo
 136     hrly_sums = hrly_sums / (24/sample_size)
 137
 138     ! 2d answer is the
 139     do ii = 1, size(buffer, 1)
 140       do j = 1, size(buffer, 2)
 141         do d = 1, sample_size
 142           buffer_exp = hrly_sums(d)
 143           if (use_mask .and. ii .eq. 1 .and. j .eq. 1) buffer_exp = -666_r4_kind
 144           if (abs(buffer(ii, j, d) - buffer_exp) > 0.0) then
 145             print *, "indices:", ii, j, d, "expected:", buffer_exp, "read in:",buffer(ii, j, d)
 146             call mpp_error(FATAL, "Check_time_diurnal::check_data_2d:: Data is not correct")
 147           endif
 148         enddo
 149       enddo
 150     enddo
 151   end subroutine check_data_2d
 152
 153   !> @brief Check that the 3d data read in is correct
 154   subroutine check_data_3d(buffer, time_level, is_regional, sample_size, nx_offset, ny_offset, nz_offset)
 155     real(kind=r4_kind), intent(in)    :: buffer(:,:,:,:) !< Buffer read from the table
 156     integer,            intent(in)    :: time_level    !< Time level read in
 157     logical,            intent(in)    :: is_regional   !< .True. if the variable is subregional
 158     integer, intent(in)               :: sample_size   !< diurnal sample size
 159     real(kind=r4_kind)                :: buffer_exp    !< Expected result
 160     integer, optional,  intent(in)    :: nx_offset     !< Offset in the x direction
 161     integer, optional,  intent(in)    :: ny_offset     !< Offset in the y direction
 162     integer, optional,  intent(in)    :: nz_offset     !< Offset in the z direction
 163
 164     integer :: ii, i, j, k, l, d!< For looping
 165     integer :: nx_oset !< Offset in the x direction (local variable)
 166     integer :: ny_oset !< Offset in the y direction (local variable)
 167     integer :: nz_oset !< Offset in the z direction (local variable)
 168     integer :: step_avg!< avg of time step increments to use in generating reference data
 169     real(r8_kind) :: hrly_sums(24/sample_size) !< can i even do this (yes)
 170     integer :: d_index !< diurnal index
 171
 172     ! data is just the hour it was sent at
 173     ! sum of hours in each diurnal section
 174     hrly_sums = 0
 175     do i=1, 23
 176       d_index = i / (24/sample_size) + 1
 177       hrly_sums(d_index) = hrly_sums(d_index) + i
 178     enddo
 179     hrly_sums = hrly_sums / (24/sample_size)
 180
 181     nx_oset = 0
 182     if (present(nx_offset)) nx_oset = nx_offset
 183
 184     ny_oset = 0
 185     if (present(ny_offset)) ny_oset = ny_offset
 186
 187     nz_oset = 0
 188     if (present(nz_offset)) nz_oset = nz_offset
 189
 190     ! 3d answer is
 191     !
 192     do ii = 1, size(buffer, 1)
 193       do j = 1, size(buffer, 2)
 194         do k = 1, size(buffer, 3)
 195           do d=1, size(buffer, 4)
 196             buffer_exp = hrly_sums(d)
 197             if (use_mask .and. ii .eq. 1 .and. j .eq. 1 .and. k .eq. 1 .and. .not. is_regional) &
 198               buffer_exp = -666_r4_kind
 199             if (abs(buffer(ii, j, k, d) - buffer_exp) > 0.0) then
 200               print *, mpp_pe(),'indices:',ii, j, k, d, "read in:", buffer(ii, j, k, d), "expected:",buffer_exp
 201               call mpp_error(FATAL, "Check_time_diurnal::check_data_3d:: Data is not correct")
 202             endif
 203           enddo
 204         enddo
 205       enddo
 206     enddo
 207   end subroutine check_data_3d
 208
 209   !> @brief Check that the 1d data read in is correct
 210   subroutine check_data_1d(buffer, time_level, sample_size)
 211     real(kind=r4_kind), intent(in)    :: buffer(:,:)   !< Buffer read from the table
 212     integer,            intent(in)    :: time_level    !< Time level read in
 213     integer,            intent(in)    :: sample_size   !< diurnal sample size of variable to check
 214     real(kind=r4_kind)                :: buffer_exp    !< Expected result
 215     integer :: ii,i, j, k, l, d!< For looping
 216     integer :: step_avg !< avg of time step increments to use in generating reference data
 217     integer :: d_index
 218     real(r8_kind) :: hrly_sums(sample_size)
 219
 220     ! sum of hours in diurnal section
 221     hrly_sums = 0
 222     do i=1, 23
 223       d_index = i / (24/sample_size) + 1
 224       hrly_sums(d_index) = hrly_sums(d_index) + i
 225     enddo
 226     hrly_sums = hrly_sums / (24/sample_size)
 227
 228     do ii = 1, size(buffer, 1)
 229       do d = 1, sample_size
 230         buffer_exp = hrly_sums(d)
 231         if (use_mask .and. ii .eq. 1) buffer_exp = -666_r4_kind
 232         if (abs(buffer(ii,d) - buffer_exp) > 0.0) then
 233           print *, "indices:", ii, d, "expected:", buffer_exp, "read in:",buffer(ii,d)
 234           call mpp_error(FATAL, "Check_time_diurnal::check_data_1d:: Data is not correct")
 235         endif
 236       enddo
 237     enddo
 238   end subroutine check_data_1d
 239
 240   !> @brief Check that the 4d data read in is correct
 241   subroutine check_data_4d(buffer, time_level, is_regional, sample_size, nx_offset, ny_offset, nz_offset)
 242     real(kind=r4_kind), intent(in)    :: buffer(:,:,:,:,:) !< Buffer read from the table
 243     integer,            intent(in)    :: time_level    !< Time level read in
 244     logical,            intent(in)    :: is_regional   !< .True. if the variable is subregional
 245     integer, intent(in)               :: sample_size   !< diurnal sample size
 246     real(kind=r4_kind)                :: buffer_exp    !< Expected result
 247     integer, optional,  intent(in)    :: nx_offset     !< Offset in the x direction
 248     integer, optional,  intent(in)    :: ny_offset     !< Offset in the y direction
 249     integer, optional,  intent(in)    :: nz_offset     !< Offset in the z direction
 250
 251     integer :: ii, i, j, k, l, d, w!< For looping
 252     integer :: nx_oset !< Offset in the x direction (local variable)
 253     integer :: ny_oset !< Offset in the y direction (local variable)
 254     integer :: nz_oset !< Offset in the z direction (local variable)
 255     integer :: step_avg!< avg of time step increments to use in generating reference data
 256     real(r8_kind) :: hrly_sums(24/sample_size) !< calculated hourly sums for each diurnal section
 257     integer :: d_index !< diurnal index
 258
 259     ! data is just the hour it was sent at
 260     ! sum of hours in each diurnal section
 261     hrly_sums = 0
 262     do i=1, 23
 263       d_index = i / (24/sample_size) + 1
 264       hrly_sums(d_index) = hrly_sums(d_index) + i
 265     enddo
 266     hrly_sums = hrly_sums / (24/sample_size)
 267
 268     nx_oset = 0
 269     if (present(nx_offset)) nx_oset = nx_offset
 270
 271     ny_oset = 0
 272     if (present(ny_offset)) ny_oset = ny_offset
 273
 274     nz_oset = 0
 275     if (present(nz_offset)) nz_oset = nz_offset
 276
 277     do ii = 1, size(buffer, 1)
 278       do j = 1, size(buffer, 2)
 279         do k = 1, size(buffer, 3)
 280           do w = 1, size(buffer, 4)
 281             do d = 1, sample_size
 282               buffer_exp = hrly_sums(d)
 283               if (use_mask .and. ii .eq. 1 .and. j .eq. 1 .and. k .eq. 1 .and. &
 284                  .not. is_regional) then
 285                 buffer_exp = -666_r4_kind
 286               endif
 287               if (abs(buffer(ii, j, k, w, d) - buffer_exp) > 0.0) then
 288                 print *, mpp_pe(),'indices:',ii, j, k, w, d, "read in:", buffer(ii, j, k, w, d), "expected:",buffer_exp
 289                 call mpp_error(FATAL, "Check_time_diurnal::check_data_4d:: Data is not correct")
 290               endif
 291             enddo
 292           enddo
 293         enddo
 294       enddo
 295     enddo
 296   end subroutine check_data_4d
 297 end program