interpolator/interpolator.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 !> @defgroup interpolator_mod interpolator_mod
  20 !> @ingroup interpolator
  21 !> @brief A module to interpolate climatology data to model the grid.
  22 !> @author William Cooke <William.Cooke@noaa.gov>
  23
  24 module interpolator_mod
  25
  26 use mpp_mod,           only : mpp_error, &
  27                               FATAL,     &
  28                               mpp_pe,    &
  29                               mpp_init,  &
  30                               mpp_exit,  &
  31                               mpp_npes,  &
  32                               WARNING,   &
  33                               NOTE,      &
  34                               input_nml_file
  35 use mpp_domains_mod,   only : mpp_domains_init,      &
  36                               mpp_update_domains,    &
  37                               mpp_define_domains,    &
  38                               mpp_global_field,      &
  39                               domain2d,              &
  40                               mpp_define_layout,     &
  41                               mpp_get_compute_domain
  42 use diag_manager_mod,  only : diag_manager_init, get_base_time, &
  43                               register_diag_field, send_data, &
  44                               diag_axis_init
  45 use fms_mod,           only : lowercase, write_version_number, &
  46                               fms_init, &
  47                               mpp_root_pe, stdlog, &
  48                               check_nml_error
  49 use fms2_io_mod,       only : FmsNetcdfFile_t, fms2_io_file_exist => file_exists, dimension_exists, &
  50                               open_file, fms2_io_read_data=>read_data,    &
  51                               variable_exists, get_variable_num_dimensions, &
  52                               get_num_variables, get_dimension_size,   &
  53                               get_variable_units, get_variable_names,  &
  54                               get_time_calendar, close_file,           &
  55                               get_variable_dimension_names, get_variable_sense
  56
  57 use horiz_interp_mod,  only : horiz_interp_type, &
  58                               horiz_interp_new,  &
  59                               horiz_interp_init, &
  60                               assignment(=), &
  61                               horiz_interp,      &
  62                               horiz_interp_del
  63 use time_manager_mod,  only : time_type,   &
  64                               set_time,    &
  65                               set_date,    &
  66                               time_type_to_real, &
  67                               days_in_year, &
  68                               get_calendar_type, &
  69                               leap_year, &
  70                               JULIAN, NOLEAP, &
  71                               get_date, &
  72                               get_date_julian, set_date_no_leap, &
  73                               set_date_julian, get_date_no_leap, &
  74                               print_date, &
  75                               operator(+), &
  76                               operator(-), &
  77                               operator(*), &
  78                               operator(>), &
  79                               operator(<), &
  80                               assignment(=), &
  81                               decrement_time
  82 use time_interp_mod,   only : time_interp, YEAR
  83 use constants_mod,     only : grav, PI, SECONDS_PER_DAY
  84 use platform_mod,      only : r4_kind, r8_kind, r16_kind, FMS_PATH_LEN, FMS_FILE_LEN
  85
  86 !--------------------------------------------------------------------
  87
  88 implicit none
  89 private
  90
  91 !---------------------------------------------------------------------
  92 !-------  interfaces --------
  93
  94 public interpolator_init, &
  95        interpolator,      &
  96        interpolate_type_eq, &
  97        obtain_interpolator_time_slices, &
  98        unset_interpolator_time_flag, &
  99        interpolator_end,  &
 100        init_clim_diag,    &
 101        query_interpolator,&
 102        read_data
 103
 104 !> Interpolates a field to a model grid
 105 !!
 106 !> Example usage:
 107 !! ~~~~~~~~~~{.f90}
 108 !! call interpolator (sulfate, model_time, p_half, model_data, name, is, js, clim_units)
 109 !! call interpolator (o3, model_time, p_half, model_data, "ozone", is, js)
 110 !! ~~~~~~~~~~
 111 !!
 112 !! The first option is used to generate sulfate models.
 113 !!
 114 !! The sulfate data is set by
 115 !! ~~~~~~~~~~{.f90}
 116 !! type(interpolate_type), intent(inout) :: sulfate
 117 !! ~~~~~~~~~~
 118 !! The name of the model is set by
 119 !! ~~~~~~~~~~{.f90}
 120 !! character(len=*), intent(in) :: name
 121 !! ~~~~~~~~~~
 122 !! The units used in this model are outputted to
 123 !! ~~~~~~~~~~{.f90}
 124 !! character(len=*), intent(out), optional :: clim_units
 125 !! ~~~~~~~~~~
 126 !!
 127 !! The second option is generate ozone models.
 128 !!
 129 !! The ozone data is set by
 130 !! ~~~~~~~~~~{.f90}
 131 !! type(interpolate_type), intent(inout) :: o3
 132 !! ~~~~~~~~~~
 133 !!
 134 !! Both of these options use the following variables in the model.
 135 !!
 136 !! The time used in the model is set by
 137 !!
 138 !! ~~~~~~~~~~{.f90}
 139 !! type(time_type), intent(in) :: model_time
 140 !! ~~~~~~~~~~
 141 !! The model pressure field is set by
 142 !! ~~~~~~~~~~{.f90}
 143 !! real, intent(in), dimension(:,:,:) :: p_half
 144 !! ~~~~~~~~~~
 145 !!
 146 !! @param [inout] <clim_type> The interpolate type previously defined by a call to interpolator_init
 147 !! @param [in] <field_name> The name of a field that you wish to interpolate
 148 !! @param [in] <Time> The model time that you wish to interpolate to
 149 !! @param [in] <phalf> The half level model pressure field
 150 !! @param [in] <is> Index for the physics window
 151 !! @param [in] <js> Index for the physics window
 152 !! @param [out] <interp_data> The model fields with the interpolated climatology data
 153 !! @param [out] <clim_units> The units of field_name
 154 !> @ingroup interpolator_mod
 155 interface interpolator
 156    module procedure interpolator_4D_r4, interpolator_4D_r8
 157    module procedure interpolator_3D_r4, interpolator_3D_r8
 158    module procedure interpolator_2D_r4, interpolator_2D_r8
 159    module procedure interpolator_4D_no_time_axis_r4, interpolator_4D_no_time_axis_r8
 160    module procedure interpolator_3D_no_time_axis_r4, interpolator_3D_no_time_axis_r8
 161    module procedure interpolator_2D_no_time_axis_r4, interpolator_2D_no_time_axis_r8
 162 end interface interpolator
 163
 164 !> Private assignment override interface for interpolate type
 165 !> @ingroup interpolator_mod
 166 interface assignment(=)
 167    module procedure interpolate_type_eq
 168 end interface
 169
 170 interface interpolator_init
 171    module procedure interpolator_init_r4
 172    module procedure interpolator_init_r8
 173 end interface interpolator_init
 174
 175 interface fms2io_interpolator_init
 176    module procedure fms2io_interpolator_init_r4
 177    module procedure fms2io_interpolator_init_r8
 178 end interface fms2io_interpolator_init
 179
 180 interface get_axis_latlon_data
 181    module procedure get_axis_latlon_data_r4
 182    module procedure get_axis_latlon_data_r8
 183 end interface get_axis_latlon_data
 184
 185 interface get_axis_level_data
 186    module procedure get_axis_level_data_r4
 187    module procedure get_axis_level_data_r8
 188 end interface get_axis_level_data
 189
 190 interface cell_center2
 191    module procedure cell_center2_r4
 192    module procedure cell_center2_r8
 193 end interface cell_center2
 194
 195 interface cart_to_latlon
 196    module procedure cart_to_latlon_r4
 197    module procedure cart_to_latlon_r8
 198 end interface cart_to_latlon
 199
 200 interface latlon2xyz
 201    module procedure latlon2xyz_r4
 202    module procedure latlon2xyz_r8
 203 end interface latlon2xyz
 204
 205 interface diag_read_data
 206   module procedure diag_read_data_r4
 207   module procedure diag_read_data_r8
 208 end interface diag_read_data
 209
 210 interface read_data
 211    module procedure read_data_r4
 212    module procedure read_data_r8
 213 end interface read_data
 214
 215 interface read_data_no_time_axis
 216    module procedure read_data_no_time_axis_r4
 217    module procedure read_data_no_time_axis_r8
 218 end interface read_data_no_time_axis
 219
 220 interface interp_linear
 221    module procedure interp_linear_r4
 222    module procedure interp_linear_r8
 223 end interface interp_linear
 224
 225 !> Private interface for weighted scalar interpolation
 226 !!
 227 !> Example usage:
 228 !! ~~~~~~~~~~{.f90}
 229 !! call interp_weighted_scalar (pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:,:),interp_data(ilon,j,:,:))
 230 !! call interp_weighted_scalar (pclim, phalf(ilon,j,:),hinterp_data(ilon,j,:),interp_data(ilon,j,:))
 231 !! ~~~~~~~~~~
 232 !!
 233 !! @param [in] <grdin> Input grid
 234 !! @param [in] <grdout> Output grid
 235 !! @param [in] <datin> Input data
 236 !! @param [out] <datout> Output data
 237 !> @ingroup interpolator_mod
 238 interface interp_weighted_scalar
 239    module procedure interp_weighted_scalar_1D_r4, interp_weighted_scalar_1D_r8
 240    module procedure interp_weighted_scalar_2D_r4, interp_weighted_scalar_2d_r8
 241 end interface interp_weighted_scalar
 242
 243 !---------------------------------------------------------------------
 244 !----------- version number for this module --------------------------
 245
 246 ! Include variable "version" to be written to log file.
 247 #include<file_version.h>
 248
 249 !> Redundant climatology data between fields
 250 !> @ingroup interpolate_type
 251
 252 type, private :: interpolate_r4_type
 253 logical :: is_allocated = .false.
 254 real(r4_kind), allocatable :: lat(:)               !< No description
 255 real(r4_kind), allocatable :: lon(:)               !< No description
 256 real(r4_kind), allocatable :: latb(:)              !< No description
 257 real(r4_kind), allocatable :: lonb(:)              !< No description
 258 real(r4_kind), allocatable :: levs(:)              !< No description
 259 real(r4_kind), allocatable :: halflevs(:)          !< No description
 260 real(r4_kind), allocatable :: data5d(:,:,:,:,:)  !< (nlatmod,nlonmod,nlevclim,size(time_init,2),nfields)
 261 real(r4_kind), allocatable :: pmon_pyear(:,:,:,:)           !< No description
 262 real(r4_kind), allocatable :: pmon_nyear(:,:,:,:)           !< No description
 263 real(r4_kind), allocatable :: nmon_nyear(:,:,:,:)           !< No description
 264 real(r4_kind), allocatable :: nmon_pyear(:,:,:,:)           !< No description
 265 real(r4_kind) :: tweight      !< No description
 266 real(r4_kind) :: tweight1     !< The time weight between the climatology years
 267 real(r4_kind) :: tweight2     !< No description
 268 real(r4_kind) :: tweight3     !< The time weight between the month
 269 end type interpolate_r4_type
 270
 271 type, private :: interpolate_r8_type
 272 logical :: is_allocated = .false.
 273 real(r8_kind), allocatable :: lat(:)               !< No description
 274 real(r8_kind), allocatable :: lon(:)               !< No description
 275 real(r8_kind), allocatable :: latb(:)              !< No description
 276 real(r8_kind), allocatable :: lonb(:)              !< No description
 277 real(r8_kind), allocatable :: levs(:)              !< No description
 278 real(r8_kind), allocatable :: halflevs(:)          !< No description
 279 real(r8_kind), allocatable :: data5d(:,:,:,:,:)  !< (nlatmod,nlonmod,nlevclim,size(time_init,2),nfields)
 280 real(r8_kind), allocatable :: pmon_pyear(:,:,:,:)           !< No description
 281 real(r8_kind), allocatable :: pmon_nyear(:,:,:,:)           !< No description
 282 real(r8_kind), allocatable :: nmon_nyear(:,:,:,:)           !< No description
 283 real(r8_kind), allocatable :: nmon_pyear(:,:,:,:)           !< No description
 284 real(r8_kind) :: tweight          !< No description
 285 real(r8_kind) :: tweight1     !< The time weight between the climatology years
 286 real(r8_kind) :: tweight2     !< No description
 287 real(r8_kind) :: tweight3     !< The time weight between the month
 288 end type interpolate_r8_type
 289
 290 type, public  :: interpolate_type
 291 private
 292 !Redundant data between fields
 293 !All climatology data
 294 type(interpolate_r4_type) :: r4_type
 295 type(interpolate_r8_type) :: r8_type
 296 type(horiz_interp_type)  :: interph                         !< No description
 297 type(time_type), allocatable :: time_slice(:) !< An array of the times within the climatology.
 298 type(FmsNetcdfFile_t)    :: fileobj       ! object that stores opened file information
 299 character(len=FMS_PATH_LEN) :: file_name     !< Climatology filename
 300 integer                  :: TIME_FLAG     !< Linear or seaonal interpolation?
 301 integer                  :: level_type    !< Pressure or Sigma level
 302 integer                  :: is,ie,js,je       !< No description
 303 integer                  :: vertical_indices !< direction of vertical
 304                                               !! data axis
 305 logical                  :: climatological_year !< Is data for year = 0000?
 306
 307 !Field specific data  for nfields
 308 character(len=64), allocatable :: field_name(:)    !< name of this field
 309 logical,           allocatable :: has_level(:)     !< indicate if the variable has level dimension
 310 integer,           allocatable :: time_init(:,:)   !< second index is the number of time_slices being
 311                                                        !! kept. 2 or ntime.
 312 integer,           allocatable :: mr(:)            !< Flag for conversion of climatology to mixing ratio.
 313 integer,           allocatable :: out_of_bounds(:) !< Flag for when surface pressure is out of bounds.
 314 !++lwh
 315 integer,           allocatable :: vert_interp(:)   !< Flag for type of vertical interpolation.
 316 !--lwh
 317 !integer                    :: indexm, indexp, climatology
 318 integer,dimension(:),  allocatable :: indexm        !< No description
 319 integer,dimension(:),  allocatable :: indexp        !< No description
 320 integer,dimension(:),  allocatable :: climatology   !< No description
 321
 322 type(time_type), allocatable :: clim_times(:,:) !< No description
 323 logical :: separate_time_vary_calc              !< No description
 324 integer :: itaum     !< No description
 325 integer :: itaup     !< No description
 326
 327 end type interpolate_type
 328
 329 !> @addtogroup interpolator_mod
 330 !> @{
 331
 332 logical            :: module_is_initialized = .false.
 333 logical            :: clim_diag_initialized = .false.
 334
 335 integer :: ndim      !< No description
 336 integer :: nvar      !< No description
 337 integer :: ntime     !< No description
 338 integer :: nlat      !< No description
 339 integer :: nlatb     !< No description
 340 integer :: nlon      !< No description
 341 integer :: nlonb     !< No description
 342 integer :: nlev      !< No description
 343 integer :: nlevh     !< No description
 344 integer ::          len, ntime_in, num_fields               !< No description
 345
 346 ! pletzer real, allocatable :: time_in(:)
 347 ! sjs real, allocatable :: climdata(:,:,:), climdata2(:,:,:)
 348
 349 character(len=64) :: units                              !< No description
 350 integer           :: sense                                        !< No description
 351
 352 integer, parameter :: max_diag_fields = 30                    !< No description
 353
 354 ! flags to indicate direction of vertical axis in  data file
 355 integer, parameter :: INCREASING_DOWNWARD = 1, INCREASING_UPWARD = -1          !< Flags to indicate direction
 356                                                                                !! of vertical axis in  data file
 357 !++lwh
 358 ! Flags to indicate whether the time interpolation should be linear or some other scheme for seasonal data.
 359 ! NOTIME indicates that data file has no time axis.
 360 integer, parameter :: LINEAR = 1, SEASONAL = 2, BILINEAR = 3, NOTIME = 4     !< Flags to indicate whether the time
 361                                                                              !! interpolation should be linear or some
 362                                                                              !! other scheme for seasonal data.
 363                                                                              !! NOTIME indicates
 364                                                                              !! that data file has no time axis.
 365
 366 ! Flags to indicate where climatology pressure levels are pressure or sigma levels
 367 integer, parameter :: PRESSURE = 1, SIGMA = 2          !< Flags to indicate where climatology pressure
 368                                                        !! levels are pressure or sigma levels
 369
 370 ! Flags to indicate whether the climatology units are mixing ratio (kg/kg) or column integral (kg/m2).
 371 ! Vertical interpolation scheme requires mixing ratio at this time.
 372 integer, parameter :: NO_CONV = 1, KG_M2 = 2          !< Flags to indicate whether the climatology units
 373                                                       !! are mixing ratio (kg/kg) or column integral (kg/m2).
 374                                                       !! Vertical interpolation scheme requires mixing ratio at
 375                                                       !! this time.
 376
 377 ! Flags to indicate what to do when the model surface pressure exceeds the  climatology surface pressure level.
 378 integer, parameter, public :: CONSTANT = 1, ZERO = 2  !< Flags to indicate what to do when the model surface
 379                                                       !! pressure exceeds the climatology surface pressure level.
 380
 381 ! Flags to indicate the type of vertical interpolation
 382 integer, parameter, public :: INTERP_WEIGHTED_P = 10, INTERP_LINEAR_P = 20, INTERP_LOG_P = 30 !< Flags to indicate
 383                                                                              !! the type of vertical interpolation
 384 !--lwh
 385
 386 real(r8_kind), parameter :: TPI = (2.0_r8_kind*PI) ! 4.*acos(0.)
 387 real(r8_kind), parameter :: DTR = TPI/360._r8_kind
 388
 389
 390
 391 integer :: num_clim_diag = 0                                            !< No description
 392 character(len=64) :: climo_diag_name(max_diag_fields)                   !< No description
 393 integer :: climo_diag_id(max_diag_fields), hinterp_id(max_diag_fields)  !< No description
 394 real(r8_kind) ::  missing_value = -1.e10_r8_kind                        !< No description
 395 ! sjs integer :: itaum, itaup
 396
 397 #ifdef ENABLE_QUAD_PRECISION
 398 ! Higher precision (kind=16) for grid geometrical factors:
 399  integer, parameter:: f_p = r16_kind    !< Higher precision (kind=16) for grid geometrical factors
 400 #else
 401 ! 64-bit precision (kind=8)
 402  integer, parameter:: f_p = r8_kind     !< 64-bit precision (kind=8)
 403 #endif
 404
 405 logical :: read_all_on_init = .false.          !< No description
 406 integer :: verbose = 0                         !< No description
 407 logical :: conservative_interp = .true.        !< No description
 408 logical :: retain_cm3_bug = .false.            !< No description
 409 logical :: use_mpp_io = .false. !< Set to true to use mpp_io, otherwise fms2io is used
 410
 411 namelist /interpolator_nml/    &
 412                              read_all_on_init, verbose, conservative_interp, retain_cm3_bug, use_mpp_io
 413
 414 contains
 415
 416 !> @brief Assignment overload routine for interpolate_type, to be used
 417 !! through the assignment interface
 418 subroutine interpolate_type_eq (Out, In)
 419
 420 type(interpolate_type), intent(in) :: In
 421 type(interpolate_type), intent(inout) :: Out
 422
 423      if(In%r4_type%is_allocated) then
 424         if (allocated(In%r4_type%lat))      Out%r4_type%lat      =  In%r4_type%lat
 425         if (allocated(In%r4_type%lon))      Out%r4_type%lon      =  In%r4_type%lon
 426         if (allocated(In%r4_type%latb))     Out%r4_type%latb     =  In%r4_type%latb
 427         if (allocated(In%r4_type%lonb))     Out%r4_type%lonb     =  In%r4_type%lonb
 428         if (allocated(In%r4_type%levs))     Out%r4_type%levs     =  In%r4_type%levs
 429         if (allocated(In%r4_type%halflevs)) Out%r4_type%halflevs =  In%r4_type%halflevs
 430      else if(In%r8_type%is_allocated) then
 431         if (allocated(In%r8_type%lat))      Out%r8_type%lat      =  In%r8_type%lat
 432         if (allocated(In%r8_type%lon))      Out%r8_type%lon      =  In%r8_type%lon
 433         if (allocated(In%r8_type%latb))     Out%r8_type%latb     =  In%r8_type%latb
 434         if (allocated(In%r8_type%lonb))     Out%r8_type%lonb     =  In%r8_type%lonb
 435         if (allocated(In%r8_type%levs))     Out%r8_type%levs     =  In%r8_type%levs
 436         if (allocated(In%r8_type%halflevs)) Out%r8_type%halflevs =  In%r8_type%halflevs
 437      end if
 438
 439      Out%interph = In%interph
 440      if (allocated(In%time_slice)) Out%time_slice =  In%time_slice
 441      Out%file_name  = In%file_name
 442      Out%time_flag  = In%time_flag
 443      Out%level_type = In%level_type
 444      Out%is = In%is
 445      Out%ie = In%ie
 446      Out%js = In%js
 447      Out%je = In%je
 448      Out%vertical_indices = In%vertical_indices
 449      Out%climatological_year = In%climatological_year
 450      if (allocated(In%has_level    )) Out%has_level     =  In%has_level
 451      if (allocated(In%field_name   )) Out%field_name    =  In%field_name
 452      if (allocated(In%time_init    )) Out%time_init     =  In%time_init
 453      if (allocated(In%mr           )) Out%mr            =  In%mr
 454      if (allocated(In%out_of_bounds)) Out%out_of_bounds =  In%out_of_bounds
 455      if (allocated(In%vert_interp  )) Out%vert_interp   =  In%vert_interp
 456      if(In%r4_type%is_allocated) then
 457         if (allocated(In%r4_type%data5d       )) Out%r4_type%data5d        =  In%r4_type%data5d
 458         if (allocated(In%r4_type%pmon_pyear   )) Out%r4_type%pmon_pyear    =  In%r4_type%pmon_pyear
 459         if (allocated(In%r4_type%pmon_nyear   )) Out%r4_type%pmon_nyear    =  In%r4_type%pmon_nyear
 460         if (allocated(In%r4_type%nmon_nyear   )) Out%r4_type%nmon_nyear    =  In%r4_type%nmon_nyear
 461         if (allocated(In%r4_type%nmon_pyear   )) Out%r4_type%nmon_pyear    =  In%r4_type%nmon_pyear
 462      else if(In%r8_type%is_allocated) then
 463         if (allocated(In%r8_type%data5d       )) Out%r8_type%data5d        =  In%r8_type%data5d
 464         if (allocated(In%r8_type%pmon_pyear   )) Out%r8_type%pmon_pyear    =  In%r8_type%pmon_pyear
 465         if (allocated(In%r8_type%pmon_nyear   )) Out%r8_type%pmon_nyear    =  In%r8_type%pmon_nyear
 466         if (allocated(In%r8_type%nmon_nyear   )) Out%r8_type%nmon_nyear    =  In%r8_type%nmon_nyear
 467         if (allocated(In%r8_type%nmon_pyear   )) Out%r8_type%nmon_pyear    =  In%r8_type%nmon_pyear
 468      end if
 469      if (allocated(In%indexm       )) Out%indexm        =  In%indexm
 470      if (allocated(In%indexp       )) Out%indexp        =  In%indexp
 471      if (allocated(In%climatology  )) Out%climatology   =  In%climatology
 472      if (allocated(In%clim_times   )) Out%clim_times    =  In%clim_times
 473      Out%separate_time_vary_calc = In%separate_time_vary_calc
 474      if(In%r4_type%is_allocated) then
 475         Out%r4_type%tweight  = In%r4_type%tweight
 476         Out%r4_type%tweight1 = In%r4_type%tweight1
 477         Out%r4_type%tweight2 = In%r4_type%tweight2
 478         Out%r4_type%tweight3 = In%r4_type%tweight3
 479      else if(In%r8_type%is_allocated) then
 480         Out%r8_type%tweight  = In%r8_type%tweight
 481         Out%r8_type%tweight1 = In%r8_type%tweight1
 482         Out%r8_type%tweight2 = In%r8_type%tweight2
 483         Out%r8_type%tweight3 = In%r8_type%tweight3
 484      end if
 485      Out%itaum = In%itaum
 486      Out%itaup = In%itaup
 487
 488      Out%r4_type%is_allocated = In%r4_type%is_allocated
 489      Out%r8_type%is_allocated = Out%r8_type%is_allocated
 490
 491 end subroutine interpolate_type_eq
 492
 493
 494
 495
 496 !#######################################################################
 497 !
 498 !---------------------------------------------------------------------
 499 !> @brief check_climo_units checks the units that the climatology
 500 !!        data is using. This is needed to allow for conversion of
 501 !!        datasets to mixing ratios which is what the vertical
 502 !!        interpolation scheme requires. The default is to assume no
 503 !!        conversion is needed. If the units are those of a column
 504 !!        burden (kg/m2) then conversion to mixing ratio is flagged.
 505 !!
 506 !! @param [in] <units> The units which you will be checking
 507 function check_climo_units(units)
 508 ! Function to check the units that the climatology data is using.
 509 ! This is needed to allow for conversion of datasets to mixing ratios which is what the
 510 ! vertical interpolation scheme requires
 511 ! The default is to assume no conversion is needed.
 512 ! If the units are those of a column burden (kg/m2) then conversion to mixing ratio is flagged.
 513 !
 514 character(len=*), intent(in) :: units
 515
 516 integer :: check_climo_units
 517
 518 check_climo_units = NO_CONV
 519 select case(chomp(units))
 520   case('kg/m2', 'kg/m^2', 'kg/m**2', 'kg m^-2', 'kg m**-2')
 521      check_climo_units = KG_M2
 522   case('molecules/cm2/s', 'molecule/cm2/s', 'molec/cm2/s')
 523      check_climo_units = KG_M2
 524   case('kg/m2/s')
 525      check_climo_units = KG_M2
 526 end select
 527
 528 end function check_climo_units
 529 !
 530 !#######################################################################
 531 !
 532 !---------------------------------------------------------------------
 533 !> @brief init_clim_diag is a routine to register diagnostic fields
 534 !!        for the climatology file. This routine calculates the domain
 535 !!        decompostion of the climatology fields for later export
 536 !!        through send_data. The ids created here are for column
 537 !!        burdens that will diagnose the vertical interpolation
 538 !!        routine.
 539 !!
 540 !! @param [inout] <clim_type> The interpolate type containing the
 541 !!                      names of the fields in the climatology file
 542 !! @param [in] <mod_axes> The axes of the model
 543 !! @param [in] <init_time> The model initialization time
 544 !!
 545 !! @throw FATAL, "init_clim_diag : You must call interpolator_init before calling init_clim_diag"
 546 !! @throw FATAL, "init_clim_diag : Trying to set up too many diagnostic fields for the climatology data"
 547 subroutine init_clim_diag(clim_type, mod_axes, init_time)
 548 !
 549 ! Routine to register diagnostic fields for the climatology file.
 550 ! This routine calculates the domain decompostion of the climatology fields
 551 ! for later export through send_data.
 552 ! The ids created here are for column burdens that will diagnose the vertical interpolation routine.
 553 ! climo_diag_id : 'module_name = climo' is intended for use with the model vertical resolution.
 554 ! hinterp_id    : 'module_name = 'hinterp' is intended for use with the climatology vertical resolution.
 555
 556 ! INTENT INOUT :
 557 !    clim_type : The interpolate type containing the names of the fields in the climatology file.
 558 !
 559 ! INTENT IN    :
 560 !   mod_axes   : The axes of the model.
 561 !   init_time  : The model initialization time.
 562 !
 563 type(interpolate_type), intent(inout)  :: clim_type
 564 integer               , intent(in)     :: mod_axes(:)
 565 type(time_type)       , intent(in)     :: init_time
 566
 567 integer :: axes(2),nxd,nyd,ndivs,i
 568 type(domain2d) :: domain
 569 integer :: domain_layout(2), iscomp, iecomp,jscomp,jecomp
 570
 571 if(clim_type%r4_type%is_allocated) call init_clim_diag_r4(clim_type, mod_axes, init_time)
 572 if(clim_type%r8_type%is_allocated) call init_clim_diag_r8(clim_type, mod_axes, init_time)
 573
 574 end subroutine init_clim_diag
 575
 576
 577
 578 !
 579 !---------------------------------------------------------------------
 580 !> @brief obtain_interpolator_time_slices makes sure that the
 581 !!        appropriate time slices are available for interpolation on
 582 !!        this time step.
 583 !!
 584 !! @param [inout] <clim_type> The interpolate type previously defined
 585 !!                      by a call to interpolator_init
 586 !! @param [in] <Time> The model time that you wish to interpolate to
 587 !!
 588 !! @throw FATAL "interpolator_timeslice 1:  file="
 589 !! @throw FATAL "interpolator_timeslice 2:  file="
 590 !! @throw FATAL "interpolator_timeslice 3:  file="
 591 !! @throw FATAL "interpolator_timeslice 4:  file="
 592 !! @throw FATAL "interpolator_timeslice 5:  file="
 593 !! @throw FATAL "interpolator_timeslice : No data from the previous
 594 !!                    climatology time but we have the next time. How did
 595 !!                    this happen?"
 596 subroutine obtain_interpolator_time_slices (clim_type, Time)
 597
 598 !  Makes sure that appropriate time slices are available for interpolation
 599 !  on this time step
 600 !
 601 ! INTENT INOUT
 602 !   clim_type   : The interpolate type previously defined by a call to interpolator_init
 603 !
 604 ! INTENT IN
 605 !   Time        : The model time that you wish to interpolate to.
 606
 607 type(interpolate_type), intent(inout)  :: clim_type
 608 type(time_type)       , intent(in)  :: Time
 609
 610 integer :: taum, taup
 611 integer :: modyear, modmonth, modday, modhour, modminute, modsecond
 612 integer :: climyear, climmonth, climday, climhour, climminute, climsecond
 613 integer :: year1, month1, day, hour, minute, second
 614 integer :: climatology, m
 615 type(time_type) :: t_prev, t_next
 616 type(time_type), dimension(2) :: month
 617 integer :: indexm, indexp, yearm, yearp
 618 integer :: i, n
 619 character(len=256) :: err_msg
 620
 621 if(clim_type%r4_type%is_allocated) call obtain_interpolator_time_slices_r4(clim_type, Time)
 622 if(clim_type%r8_type%is_allocated) call obtain_interpolator_time_slices_r8(clim_type, Time)
 623
 624
 625 end subroutine obtain_interpolator_time_slices
 626
 627
 628 !#####################################################################
 629 !
 630 !---------------------------------------------------------------------
 631 !> @brief unset_interpolator_time_flag sets a flag in clim_type to
 632 !!        false.
 633 !!
 634 !! @param [inout] <clim_type> The interpolate type containing the names of the fields in the climatology file
 635 subroutine unset_interpolator_time_flag (clim_type)
 636
 637 type(interpolate_type), intent(inout) :: clim_type
 638
 639
 640       clim_type%separate_time_vary_calc = .false.
 641
 642
 643 end subroutine unset_interpolator_time_flag
 644
 645
 646 !#####################################################################
 647 !
 648 !---------------------------------------------------------------------
 649 !> @brief interpolator_end receives interpolate data as input
 650 !!        and deallocates its memory.
 651 !!
 652 !! @param [inout] <clim_type> The interpolate type whose components will be deallocated
 653 subroutine interpolator_end(clim_type)
 654 ! Subroutine to deallocate the interpolate type clim_type.
 655 !
 656 ! INTENT INOUT
 657 !  clim_type : allocate type whose components will be deallocated.
 658 !
 659 type(interpolate_type), intent(inout) :: clim_type
 660 integer :: logunit
 661
 662 logunit=stdlog()
 663 if ( mpp_pe() == mpp_root_pe() ) then
 664    write (logunit,'(/,(a))') 'Exiting interpolator, have a nice day ...'
 665 end if
 666
 667 if(clim_type%r4_type%is_allocated) then
 668    if (allocated (clim_type%r4_type%lat     )) deallocate(clim_type%r4_type%lat)
 669    if (allocated (clim_type%r4_type%lon     )) deallocate(clim_type%r4_type%lon)
 670    if (allocated (clim_type%r4_type%latb    )) deallocate(clim_type%r4_type%latb)
 671    if (allocated (clim_type%r4_type%lonb    )) deallocate(clim_type%r4_type%lonb)
 672    if (allocated (clim_type%r4_type%levs    )) deallocate(clim_type%r4_type%levs)
 673    if (allocated (clim_type%r4_type%halflevs)) deallocate(clim_type%r4_type%halflevs)
 674    if (allocated (clim_type%r4_type%data5d  )) deallocate(clim_type%r4_type%data5d)
 675 else if(clim_type%r8_type%is_allocated) then
 676    if (allocated (clim_type%r8_type%lat     )) deallocate(clim_type%r8_type%lat)
 677    if (allocated (clim_type%r8_type%lon     )) deallocate(clim_type%r8_type%lon)
 678    if (allocated (clim_type%r8_type%latb    )) deallocate(clim_type%r8_type%latb)
 679    if (allocated (clim_type%r8_type%lonb    )) deallocate(clim_type%r8_type%lonb)
 680    if (allocated (clim_type%r8_type%levs    )) deallocate(clim_type%r8_type%levs)
 681    if (allocated (clim_type%r8_type%halflevs)) deallocate(clim_type%r8_type%halflevs)
 682    if (allocated (clim_type%r8_type%data5d))   deallocate(clim_type%r8_type%data5d)
 683 end if
 684
 685 if (allocated (clim_type%time_slice)) deallocate(clim_type%time_slice)
 686 if (allocated (clim_type%field_name)) deallocate(clim_type%field_name)
 687 if (allocated (clim_type%time_init )) deallocate(clim_type%time_init)
 688 if (allocated (clim_type%has_level))  deallocate(clim_type%has_level)
 689 if (allocated (clim_type%mr        )) deallocate(clim_type%mr)
 690 if (allocated (clim_type%out_of_bounds )) deallocate(clim_type%out_of_bounds)
 691 if (allocated (clim_type%vert_interp ))   deallocate(clim_type%vert_interp)
 692 if (allocated(clim_type%indexm)) deallocate(clim_type%indexm)
 693 if (allocated(clim_type%indexp)) deallocate(clim_type%indexp)
 694 if (allocated(clim_type%clim_times))  deallocate(clim_type%clim_times)
 695 if (allocated(clim_type%climatology)) deallocate(clim_type%climatology)
 696
 697 call horiz_interp_del(clim_type%interph)
 698
 699 if(clim_type%r4_type%is_allocated) then
 700    if (allocated(clim_type%r4_type%pmon_pyear)) then
 701       deallocate(clim_type%r4_type%pmon_pyear)
 702       deallocate(clim_type%r4_type%pmon_nyear)
 703       deallocate(clim_type%r4_type%nmon_nyear)
 704       deallocate(clim_type%r4_type%nmon_pyear)
 705    end if
 706 else if(clim_type%r8_type%is_allocated) then
 707    if (allocated(clim_type%r8_type%pmon_pyear)) then
 708       deallocate(clim_type%r8_type%pmon_pyear)
 709       deallocate(clim_type%r8_type%pmon_nyear)
 710       deallocate(clim_type%r8_type%nmon_nyear)
 711       deallocate(clim_type%r8_type%nmon_pyear)
 712    end if
 713 endif
 714
 715 clim_type%r4_type%is_allocated=.false.
 716 clim_type%r8_type%is_allocated=.false.
 717
 718 !! RSH mod
 719 if( .not.(clim_type%TIME_FLAG .eq. LINEAR  .and. read_all_on_init) &
 720    .and. (clim_type%TIME_FLAG.ne.NOTIME) ) then
 721 !     read_all_on_init)) .or. clim_type%TIME_FLAG .eq. BILINEAR  ) then
 722      call close_file(clim_type%fileobj)
 723 endif
 724
 725
 726 module_is_initialized = .false.
 727
 728 end subroutine interpolator_end
 729 !
 730 !#######################################################################
 731 !
 732 !++lwh
 733 !
 734 !---------------------------------------------------------------------
 735 !> @brief query_interpolator receives an interpolate type as input
 736 !!        and returns the number of fields and field names.
 737 !!
 738 !! @param [in] <clim_type> The interpolate type which contains the data
 739 !! @param [out] <nfields> OPTIONAL: No description
 740 !! @param [out] <field_names> OPTIONAL: No description
 741 subroutine query_interpolator( clim_type, nfields, field_names )
 742 !
 743 ! Query an interpolate_type variable to find the number of fields and field names.
 744 !
 745 type(interpolate_type), intent(in)                    :: clim_type
 746 integer, intent(out), optional                        :: nfields
 747 character(len=*), dimension(:), intent(out), optional :: field_names
 748
 749 if( present( nfields ) )     nfields     = SIZE( clim_type%field_name(:) )
 750 if( present( field_names ) ) field_names = clim_type%field_name
 751
 752 end subroutine query_interpolator
 753 !--lwh
 754 !
 755 !#######################################################################
 756 !
 757 !---------------------------------------------------------------------
 758 !> @brief chomp receives a string from NetCDF files and removes
 759 !!        CHAR(0) from the end of this string.
 760 !!
 761 !! @param [in] <string> The string from the NetCDF file
 762 function chomp(string)
 763 !
 764 ! A function to remove CHAR(0) from the end of strings read from NetCDF files.
 765 !
 766 character(len=*), intent(in) :: string
 767 character(len=64) :: chomp
 768
 769 integer :: len
 770
 771 len = len_trim(string)
 772 if (string(len:len) == CHAR(0)) len = len -1
 773
 774 chomp = string(:len)
 775
 776 end function chomp
 777 !
 778 !########################################################################
 779
 780 #include "interpolator_r4.fh"
 781 #include "interpolator_r8.fh"
 782
 783 end module interpolator_mod
 784
 785 !> @}
 786 ! close documentation grouping
 787 !
 788 !#######################################################################