test_fms/diag_manager/test_diag_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 send more complicated data than just the current hour
  20
  21 !> @brief  Program to test the diurnal reduction
  22 !! Similar to test_reduction_methods, but uses the variables and reductions
  23 !! from the test_diag_manager_time diurnal test (#25)
  24 program test_diag_diurnal
  25   use fms_mod,           only: fms_init, fms_end
  26   use testing_utils,     only: allocate_buffer, test_normal, test_openmp, test_halos, no_mask, logical_mask, real_mask
  27   use platform_mod,      only: r8_kind
  28   use block_control_mod, only: block_control_type, define_blocks
  29   use mpp_mod,           only: mpp_sync, FATAL, mpp_error, mpp_npes, mpp_pe, mpp_root_pe, mpp_broadcast, input_nml_file
  30   use time_manager_mod,  only: time_type, set_calendar_type, set_date, JULIAN, set_time, OPERATOR(+), days_in_month, &
  31                                get_time
  32   use diag_manager_mod,  only: diag_manager_init, diag_manager_end, diag_axis_init, register_diag_field, &
  33                                diag_send_complete, diag_manager_set_time_end, send_data
  34   use mpp_domains_mod,   only: domain2d, mpp_define_domains, mpp_define_io_domain, mpp_get_compute_domain, &
  35                                mpp_get_data_domain
  36
  37   implicit none
  38
  39   integer                            :: nx              !< Number of points in the x direction
  40   integer                            :: ny              !< Number of points in the y direction
  41   integer                            :: nz              !< Number of points in the z direction
  42   integer                            :: nw              !< Number of points in the 4th dimension
  43   integer                            :: layout(2)       !< Layout
  44   integer                            :: io_layout(2)    !< Io layout
  45   type(domain2d)                     :: Domain          !< 2D domain
  46   integer                            :: isc, isd        !< Starting x compute, data domain index
  47   integer                            :: iec, ied        !< Ending x compute, data domain index
  48   integer                            :: jsc, jsd        !< Starting y compute, data domaine index
  49   integer                            :: jec, jed        !< Ending y compute, data domain index
  50   integer                            :: nhalox          !< Number of halos in x
  51   integer                            :: nhaloy          !< Number of halos in y
  52   real(kind=r8_kind), allocatable    :: cdata(:,:,:,:)  !< Data in the compute domain
  53   real(kind=r8_kind), allocatable    :: ddata(:,:,:,:)  !< Data in the data domain
  54   real(kind=r8_kind), allocatable    :: crmask(:,:,:,:) !< Mask in the compute domain
  55   real(kind=r8_kind), allocatable    :: drmask(:,:,:,:) !< Mask in the data domain
  56   logical,            allocatable    :: clmask(:,:,:,:) !< Logical mask in the compute domain
  57   logical,            allocatable    :: dlmask(:,:,:,:) !< Logical mask in the data domain
  58   type(time_type)                    :: Time            !< Time of the simulation
  59   type(time_type)                    :: Time_step       !< Time of the simulation
  60   integer                            :: nmonths         !< number of months to run for (submits ntimes per month)
  61   integer                            :: ndays           !< number of days in the month
  62   integer                            :: nhours          !< number of hours in a day - 1
  63   integer                            :: id_x            !< axis id for the x dimension
  64   integer                            :: id_y            !< axis id for the y dimension
  65   integer                            :: id_z            !< axis id for the z dimension
  66   integer                            :: id_w            !< axis id for the w dimension
  67   integer                            :: id_var0         !< diag_field id for 0d var
  68   integer                            :: id_var1         !< diag_field id for 1d var
  69   integer                            :: id_var2         !< diag_field id for 2d var
  70   integer                            :: id_var3         !< diag_field id for 3d var
  71   integer                            :: id_var4         !< diag_field id for 4d var
  72   integer                            :: io_status       !< Status after reading the namelist
  73   type(block_control_type)           :: my_block        !< Returns instantiated @ref block_control_type
  74   logical                            :: message         !< Flag for outputting debug message
  75   integer                            :: isd1            !< Starting x data domain index (1-based)
  76   integer                            :: ied1            !< Ending x data domain index (1-based)
  77   integer                            :: jsd1            !< Starting y data domain index (1-based)
  78   integer                            :: jed1            !< Ending y data domain index (1-based)
  79   integer                            :: isw             !< Starting index for each thread in the x direction
  80   integer                            :: iew             !< Ending index for each thread in the x direction
  81   integer                            :: jsw             !< Starting index for each thread in the y direction
  82   integer                            :: jew             !< Ending index for each thread in the y direction
  83   integer                            :: is1             !< Starting index for each thread in the x direction (1-based)
  84   integer                            :: ie1             !< Ending index for each thread in the x direction (1-based)
  85   integer                            :: js1             !< Starting index for each thread in the y direction (1-based)
  86   integer                            :: je1             !< Ending index for each thread in the y direction (1-based)
  87   integer                            :: iblock          !< For looping through the blocks
  88   integer                            :: i               !< For do loops
  89   logical                            :: used            !< Dummy argument to send_data
  90   real(kind=r8_kind)                 :: missing_value   !< Missing value to use
  91   integer :: days_out, seconds_out
  92   integer :: m, h, d !< to iterate through months, hours, and days
  93
  94   !< Configuration parameters
  95   integer :: test_case = test_normal !< Indicates which test case to run
  96   integer :: mask_case = no_mask     !< Indicates which masking option to run
  97
  98   namelist / test_diag_diurnal_nml / test_case, mask_case
  99
 100   call fms_init
 101   call set_calendar_type(JULIAN)
 102   call diag_manager_init
 103
 104   read (input_nml_file, test_diag_diurnal_nml, iostat=io_status)
 105   if (io_status > 0) call mpp_error(FATAL,'=>test_modern_diag: Error reading input.nml')
 106
 107   nx = 96
 108   ny = 96
 109   nz = 5
 110   nw = 2
 111   layout = (/1, mpp_npes()/)
 112   io_layout = (/1, 1/)
 113   nhalox = 2
 114   nhaloy = 2
 115   nmonths = 3
 116   nhours = 23 !< Number of hours in a day - 1
 117
 118   !< Create a lat/lon domain
 119   call mpp_define_domains( (/1,nx,1,ny/), layout, Domain, name='2D domain', xhalo=nhalox, yhalo=nhaloy)
 120   call mpp_define_io_domain(Domain, io_layout)
 121   call mpp_get_compute_domain(Domain, isc, iec, jsc, jec)
 122   call mpp_get_data_domain(Domain, isd, ied, jsd, jed)
 123
 124   cdata = allocate_buffer(isc, iec, jsc, jec, nz, nw)
 125
 126   select case (test_case)
 127   case (test_normal)
 128     if (mpp_pe() .eq. mpp_root_pe()) print *, "Testing the normal send_data calls"
 129   case (test_halos)
 130     if (mpp_pe() .eq. mpp_root_pe()) print *, "Testing the send_data calls with halos"
 131     ddata = allocate_buffer(isd, ied, jsd, jed, nz, nw)
 132   case (test_openmp)
 133     if (mpp_pe() .eq. mpp_root_pe()) print *, "Testing the send_data calls with openmp blocks"
 134      call define_blocks ('testing_model', my_block, isc, iec, jsc, jec, kpts=0, &
 135                          nx_block=1, ny_block=4, message=message)
 136   end select
 137
 138   select case (mask_case)
 139   case (logical_mask)
 140     clmask = allocate_logical_mask(isc, iec, jsc, jec, nz, nw)
 141     if (mpp_pe() .eq. 0) clmask(isc, jsc, 1, :) = .False.
 142
 143     if (test_case .eq. test_halos) then
 144       dlmask = allocate_logical_mask(isd, ied, jsd, jed, nz, nw)
 145       if (mpp_pe() .eq. 0) dlmask(1+nhalox, 1+nhaloy, 1, :) = .False.
 146     endif
 147   case (real_mask)
 148     crmask = allocate_real_mask(isc, iec, jsc, jec, nz, nw)
 149     if (mpp_pe() .eq. 0) crmask(isc, jsc, 1, :) = 0_r8_kind
 150
 151     if (test_case .eq. test_halos) then
 152       drmask = allocate_real_mask(isd, ied, jsd, jed, nz, nw)
 153       if (mpp_pe() .eq. 0) drmask(1+nhalox, 1+nhaloy, 1, :) = 0_r8_kind
 154     endif
 155   end select
 156
 157
 158   !< Get the data domain indices (1 based)
 159   isd1 = isc-isd+1
 160   jsd1 = jsc-jsd+1
 161   ied1 = isd1 + iec-isc
 162   jed1 = jsd1 + jec-jsc
 163
 164   !< set up end time
 165   Time = set_date(2,1,1,0,0,0)
 166   Time_step = set_time (3600,0) !< 1 hour
 167   call diag_manager_set_time_end(set_date(2,nmonths+1,1,0,0,0))
 168
 169   !< Register the axis
 170   id_x  = diag_axis_init('x',  real((/ (i, i = 1,nx) /), kind=r8_kind),  'point_E', 'x', long_name='point_E', &
 171     Domain2=Domain)
 172   id_y  = diag_axis_init('y',  real((/ (i, i = 1,ny) /), kind=r8_kind),  'point_N', 'y', long_name='point_N', &
 173     Domain2=Domain)
 174   id_z  = diag_axis_init('z',  real((/ (i, i = 1,nz) /), kind=r8_kind),  'point_Z', 'z', long_name='point_Z')
 175   id_w  = diag_axis_init('w',  real((/ (i, i = 1,nw) /), kind=r8_kind),  'point_W', 'n', long_name='point_W')
 176
 177   missing_value = -666._r8_kind
 178   !< Register the fields
 179   id_var1 = register_diag_field  ('ocn_mod', 'var1', (/id_x/), Time, 'var1', &
 180     'mullions', missing_value = missing_value)
 181   id_var2 = register_diag_field  ('ocn_mod', 'var2', (/id_x, id_y/), Time, 'var2', &
 182     'mullions', missing_value = missing_value)
 183   id_var3 = register_diag_field  ('ocn_mod', 'var3', (/id_x, id_y, id_z/), Time, 'var3', &
 184     'mullions', missing_value = missing_value)
 185   id_var4 = register_diag_field  ('ocn_mod', 'var4', (/id_x, id_y, id_z, id_w/), Time, 'var4', &
 186     'mullions', missing_value = missing_value)
 187
 188   ! iterate through nmonths and each day, each hour
 189   do m = 1, nmonths + 1
 190     Time = set_date(2,m,1)
 191     ndays = days_in_month(Time)
 192     if (m .eq. nmonths + 1) then
 193       ! This it so  that is can run till (2 4 1 0 0 0)
 194       ndays = 1
 195       nhours = 0
 196     endif
 197     do d = 1, ndays
 198       do h = 0, nhours ! hours
 199         Time = set_date(2,m,d,hour=h)
 200
 201         call set_buffer(cdata, m, d, h)
 202
 203         select case(test_case)
 204         case (test_normal)
 205           select case (mask_case)
 206           case (no_mask)
 207             used = send_data(id_var1, cdata(:,1,1,1), Time)
 208             used = send_data(id_var2, cdata(:,:,1,1), Time)
 209             used = send_data(id_var3, cdata(:,:,:,1), Time)
 210             used = send_data(id_var4, cdata(:,:,:,:), Time)
 211           case (real_mask)
 212             used = send_data(id_var1, cdata(:,1,1,1), Time, rmask=crmask(:,1,1,1))
 213             used = send_data(id_var2, cdata(:,:,1,1), Time, rmask=crmask(:,:,1,1))
 214             used = send_data(id_var3, cdata(:,:,:,1), Time, rmask=crmask(:,:,:,1))
 215             used = send_data(id_var4, cdata(:,:,:,:), Time, rmask=crmask(:,:,:,:))
 216           case (logical_mask)
 217             used = send_data(id_var1, cdata(:,1,1,1), Time, mask=clmask(:,1,1,1))
 218             used = send_data(id_var2, cdata(:,:,1,1), Time, mask=clmask(:,:,1,1))
 219             used = send_data(id_var3, cdata(:,:,:,1), Time, mask=clmask(:,:,:,1))
 220             used = send_data(id_var4, cdata(:,:,:,:), Time, mask=clmask(:,:,:,:))
 221           end select
 222         case (test_halos)
 223           call set_buffer(ddata, m, d, h)
 224           select case (mask_case)
 225           case (no_mask)
 226             used = send_data(id_var1, ddata(:,1,1,1), Time)
 227             used = send_data(id_var2, ddata(:,:,1,1), Time, &
 228               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1)
 229             used = send_data(id_var3, ddata(:,:,:,1), Time, &
 230               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1)
 231             used = send_data(id_var4, ddata(:,:,:,:), Time, &
 232               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1)
 233           case (real_mask)
 234             used = send_data(id_var1, ddata(:,1,1,1), Time, &
 235               rmask=drmask(:,1,1,1))
 236             used = send_data(id_var2, ddata(:,:,1,1), Time, &
 237               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1, &
 238               rmask=drmask(:,:,1,1))
 239             used = send_data(id_var3, ddata(:,:,:,1), Time, &
 240               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1, &
 241               rmask=drmask(:,:,:,1))
 242             used = send_data(id_var4, ddata(:,:,:,:), Time, &
 243               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1, &
 244               rmask=drmask(:,:,:,:))
 245           case (logical_mask)
 246             used = send_data(id_var1, ddata(:,1,1,1), Time, &
 247               mask=dlmask(:,1,1,1))
 248             used = send_data(id_var2, ddata(:,:,1,1), Time, &
 249               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1, &
 250               mask=dlmask(:,:,1,1))
 251             used = send_data(id_var3, ddata(:,:,:,1), Time, &
 252               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1, &
 253               mask=dlmask(:,:,:,1))
 254             used = send_data(id_var4, ddata(:,:,:,:), Time, &
 255               is_in=isd1, ie_in=ied1, js_in=jsd1, je_in=jed1, &
 256               mask=dlmask(:,:,:,:))
 257           end select
 258         case (test_openmp)
 259           select case(mask_case)
 260           case (no_mask)
 261             used=send_data(id_var1, cdata(:, 1, 1, 1), time)
 262           case (logical_mask)
 263             used=send_data(id_var1, cdata(:, 1, 1, 1), time, &
 264                 mask=clmask(:, 1, 1, 1))
 265           case (real_mask)
 266             used=send_data(id_var1, cdata(:, 1, 1, 1), time, &
 267                 rmask=crmask(:, 1, 1, 1))
 268           end select
 269 !$OMP parallel do default(shared) private(iblock, isw, iew, jsw, jew, is1, ie1, js1, je1)
 270           do iblock=1, 4
 271             isw = my_block%ibs(iblock)
 272             jsw = my_block%jbs(iblock)
 273             iew = my_block%ibe(iblock)
 274             jew = my_block%jbe(iblock)
 275
 276           !--- indices for 1-based arrays ---
 277             is1 = isw-isc+1
 278             ie1 = iew-isc+1
 279             js1 = jsw-jsc+1
 280             je1 = jew-jsc+1
 281
 282             select case (mask_case)
 283             case (no_mask)
 284               used=send_data(id_var2, cdata(is1:ie1, js1:je1, 1, 1), time, is_in=is1, js_in=js1)
 285               used=send_data(id_var3, cdata(is1:ie1, js1:je1, :, 1), time, is_in=is1, js_in=js1)
 286               used=send_data(id_var4, cdata(is1:ie1, js1:je1, :, :), time, is_in=is1, js_in=js1)
 287             case (real_mask)
 288               used=send_data(id_var2, cdata(is1:ie1, js1:je1, 1, 1), time, is_in=is1, js_in=js1, &
 289                 rmask=crmask(is1:ie1, js1:je1, 1, 1))
 290               used=send_data(id_var3, cdata(is1:ie1, js1:je1, :, 1), time, is_in=is1, js_in=js1, &
 291                 rmask=crmask(is1:ie1, js1:je1, :, 1))
 292               used=send_data(id_var4, cdata(is1:ie1, js1:je1, :, :), time, is_in=is1, js_in=js1, &
 293                 rmask=crmask(is1:ie1, js1:je1, :, :))
 294             case (logical_mask)
 295               used=send_data(id_var2, cdata(is1:ie1, js1:je1, 1, 1), time, is_in=is1, js_in=js1, &
 296                 mask=clmask(is1:ie1, js1:je1, 1, 1))
 297               used=send_data(id_var3, cdata(is1:ie1, js1:je1, :, 1), time, is_in=is1, js_in=js1, &
 298                 mask=clmask(is1:ie1, js1:je1, :, 1))
 299               used=send_data(id_var4, cdata(is1:ie1, js1:je1, :, :), time, is_in=is1, js_in=js1, &
 300                 mask=clmask(is1:ie1, js1:je1, :, :))
 301             end select
 302           enddo
 303         end select
 304         call diag_send_complete(Time_step)
 305       enddo
 306     enddo
 307   enddo
 308
 309   call diag_manager_end(Time)
 310
 311   call fms_end
 312
 313   contains
 314
 315   !> @brief Allocate the logical mask based on the starting/ending indices
 316   !! @return logical mask initiliazed to .True.
 317   function allocate_logical_mask(is, ie, js, je, k, l) &
 318   result(buffer)
 319     integer, intent(in) :: is !< Starting x index
 320     integer, intent(in) :: ie !< Ending x index
 321     integer, intent(in) :: js !< Starting y index
 322     integer, intent(in) :: je !< Ending y index
 323     integer, intent(in) :: k  !< Number of points in the 4th dimension
 324     integer, intent(in) :: l  !< Number of points in the 5th dimension
 325
 326     logical, allocatable :: buffer(:,:,:,:)
 327
 328     allocate(buffer(is:ie, js:je, 1:k, 1:l))
 329     buffer = .True.
 330   end function allocate_logical_mask
 331
 332   !> @brief Allocate the real mask based on the starting/ending indices
 333   !! @returnreal mask initiliazed to 1_r8_kind
 334   function allocate_real_mask(is, ie, js, je, k, l) &
 335   result(buffer)
 336     integer, intent(in) :: is !< Starting x index
 337     integer, intent(in) :: ie !< Ending x index
 338     integer, intent(in) :: js !< Starting y index
 339     integer, intent(in) :: je !< Ending y index
 340     integer, intent(in) :: k  !< Number of points in the 4th dimension
 341     integer, intent(in) :: l  !< Number of points in the 5th dimension
 342     real(kind=r8_kind), allocatable :: buffer(:,:,:,:)
 343
 344     allocate(buffer(is:ie, js:je, 1:k, 1:l))
 345     buffer = 1.0_r8_kind
 346   end function allocate_real_mask
 347
 348
 349   !> @brief Set the buffer based on the time_index
 350   subroutine set_buffer(buffer, month, day, hour)
 351     real(kind=r8_kind), intent(inout) :: buffer(:,:,:,:) !< Output buffer
 352     integer,            intent(in)    :: month, day, hour !< Time index
 353
 354     buffer = hour ! month * 10000 + day * 100 + hour
 355
 356   end subroutine set_buffer
 357
 358 end program test_diag_diurnal