amip_interp/amip_interp.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 !> @defgroup amip_interp_mod amip_interp_mod
  21 !> @ingroup amip_interp
  22 !> @brief Provides observed sea surface temperature and ice mask data sets that have been
  23 !! interpolated onto your model's grid.
  24 !!
  25 !> @author Bruce Wyman
  26 !!
  27 !> When using these routines three possible data sets are available:
  28 !!
  29 !! 1. AMIP @link http://www-pcmdi.llnl.gov/amip @endlink from Jan 1979 to Jan 1989 (2 deg x 2 deg)
  30 !! 2. Reynolds OI @link amip_interp.rey_oi.txt @endlink from Nov 1981 to Jan 1999 (1 deg x 1 deg)
  31 !! 3. Reynolds EOF @link ftp://podaac.jpl.nasa.gov/pub/sea_surface_temperature/reynolds/rsst/doc/rsst.html
  32 !! @endlink from Jan 1950 to Dec 1998 (2 deg x 2 deg)
  33 !!
  34 !! All original data are observed monthly means. This module
  35 !! interpolates linearly in time between pairs of monthly means.
  36 !! Horizontal interpolation is done using the horiz_interp module.
  37 !!
  38 !! When a requested date falls outside the range of dates available
  39 !! a namelist option allows for use of the climatological monthly
  40 !! mean values which are computed from all of the data in a particular
  41 !! data set. \n
  42 !! \n AMIP 1:\n
  43 !!   from Jan 1979 to Jan 1989 (2 deg x 2 deg).\n\n
  44 !! Reynolds OI:\n
  45 !!   from Nov 1981 to Jan 1999 (1 deg x 1 deg)\n
  46 !!   The analysis uses in situ and satellite SST's plus
  47 !!   SST's simulated by sea-ice cover.\n\n
  48 !! Reynold's EOF:\n
  49 !!   from Jan 1950 to Dec 1998 (2 deg x 2 deg)\n
  50 !!   NCEP Reynolds Historical Reconstructed Sea Surface Temperature
  51 !!   The analysis uses both in-situ SSTs and satellite derived SSTs
  52 !!   from the NOAA Advanced Very High Resolution Radiometer.
  53 !!   In-situ data is used from 1950 to 1981, while both AVHRR derived
  54 !!   satellite SSTs and in-situ data are used from 1981 to the
  55 !!   end of 1998.
  56 !!
  57 !> @note The data set used by this module have been reformatted as 32-bit IEEE.
  58 !!   The data values are packed into 16-bit integers.
  59 !!
  60 !!   The data sets are read from the following files:
  61 !!
  62 !!         amip1           INPUT/amip1_sst.data
  63 !!         reynolds_io     INPUT/reyoi_sst.data
  64 !!         reynolds_eof    INPUT/reynolds_sst.data
  65
  66
  67 !> @file
  68 !> File for amip_interp_mod
  69
  70 !> @addtogroup amip_interp_mod
  71 !> @{
  72 module amip_interp_mod
  73
  74 use  time_interp_mod, only: time_interp, fraction_of_year
  75
  76 use time_manager_mod, only: time_type, operator(+), operator(>), &
  77                              get_date, set_time, set_date
  78
  79 ! add by JHC
  80 use get_cal_time_mod, only: get_cal_time
  81
  82 ! end add by JHC
  83
  84 use  horiz_interp_mod, only: horiz_interp_init, horiz_interp,  &
  85                              horiz_interp_new, horiz_interp_del, &
  86                              horiz_interp_type, assignment(=)
  87
  88 use           fms_mod, only: error_mesg, write_version_number,  &
  89                              NOTE, WARNING, FATAL, stdlog, check_nml_error, &
  90                              mpp_pe, lowercase, mpp_root_pe,    &
  91                              NOTE, mpp_error, fms_error_handler
  92
  93 use     constants_mod, only: TFREEZE, pi
  94 use      platform_mod, only: R4_KIND, I2_KIND
  95 use mpp_mod,           only: input_nml_file
  96 use fms2_io_mod,       only: FmsNetcdfFile_t, fms2_io_file_exists=>file_exists, open_file, close_file, &
  97                              get_dimension_size, fms2_io_read_data=>read_data
  98
  99 implicit none
 100 private
 101
 102 !-----------------------------------------------------------------------
 103 !----------------- Public interfaces -----------------------------------
 104
 105 public amip_interp_init, get_amip_sst, get_amip_ice, amip_interp_new, &
 106        amip_interp_del, amip_interp_type, assignment(=)
 107
 108 !-----------------------------------------------------------------------
 109 !----------------- Public Data -----------------------------------
 110 integer :: i_sst = 1200
 111 integer :: j_sst = 600
 112 real, parameter:: big_number = 1.E30
 113 logical :: forecast_mode = .false.
 114 real, allocatable, dimension(:,:) ::  sst_ncep, sst_anom
 115
 116 public i_sst, j_sst, sst_ncep, sst_anom, forecast_mode, use_ncep_sst
 117
 118 !-----------------------------------------------------------------------
 119 !--------------------- private below here ------------------------------
 120
 121 !  ---- version number -----
 122
 123 ! Include variable "version" to be written to log file.
 124 #include<file_version.h>
 125
 126    real, allocatable:: temp1(:,:), temp2(:,:)
 127 ! add by JHC
 128    real, allocatable, dimension(:,:) :: tempamip
 129 ! end add by JHC
 130 !-----------------------------------------------------------------------
 131 !------ private defined data type --------
 132
 133 !> @}
 134
 135 !> @brief Private data type for representing a calendar date
 136 !> @ingroup amip_interp_mod
 137 type date_type
 138    sequence
 139    integer :: year, month, day
 140 end type
 141
 142 !> Assignment overload to allow native assignment between amip_interp_type variables.
 143 !> @ingroup amip_interp_mod
 144 interface assignment(=)
 145   module procedure  amip_interp_type_eq
 146 end interface
 147
 148 !> Private logical equality overload for amip_interp_type
 149 !> @ingroup amip_interp_mod
 150 interface operator (==)
 151    module procedure date_equals
 152 end interface
 153
 154 !> Private logical inequality overload for amip_interp_type
 155 !> @ingroup amip_interp_mod
 156 interface operator (/=)
 157    module procedure date_not_equals
 158 end interface
 159
 160 !> Private logical greater than overload for amip_interp_type
 161 !> @ingroup amip_interp_mod
 162 interface operator (>)
 163    module procedure date_gt
 164 end interface
 165
 166
 167 !> Initializes data needed for the horizontal
 168 !! interpolation between the sst data and model grid.
 169 !!
 170 !> The returned variable of type amip_interp_type is needed when
 171 !! calling get_amip_sst and get_amip_ice.
 172 !!
 173 !> @param lon
 174 !!     Longitude in radians of the model's grid box edges (1d lat/lon grid case)
 175 !!     or at grid box mid-point (2d case for arbitrary grids).
 176 !> @param lat
 177 !!     Latitude in radians of the model's grid box edges (1d lat/lon grid case)
 178 !!     or at grid box mid-point (2d case for arbitrary grids).
 179 !> @param mask
 180 !!     A mask for the model grid.
 181 !> @param use_climo
 182 !!     Flag the specifies that monthly mean climatological values will be used.
 183 !> @param use_annual
 184 !!     Flag the specifies that the annual mean climatological
 185 !!              will be used.  If both use_annual = use_climo = true,
 186 !!              then use_annual = true will be used.
 187 !> @param interp_method
 188 !!     specify the horiz_interp scheme. = "conservative" means conservative scheme,
 189 !!     = "bilinear" means  bilinear interpolation.
 190 !!
 191 !> @return interp, a defined data type variable needed when calling get_amip_sst and get_amip_ice.
 192 !!
 193 !! \n Example usage:
 194 !!
 195 !!                 Interp = amip_interp_new ( lon, lat, mask, use_climo, use_annual, interp_method )
 196 !!
 197 !! This function may be called to initialize multiple variables
 198 !! of type amip_interp_type.  However, there currently is no
 199 !! call to release the storage used by this variable.
 200 !!
 201 !! The size of input augment mask must be a function of the size
 202 !! of input augments lon and lat. The first and second dimensions
 203 !! of mask must equal (size(lon,1)-1, size(lat,2)-1).
 204 !!
 205 !> @throws "FATAL: the value of the namelist parameter DATA_SET being used is not allowed"
 206 !! Check the value of namelist variable DATA_SET.
 207 !!
 208 !> @throws "FATAL: requested input data set does not exist"
 209 !! The data set requested is valid but the data does not exist in
 210 !! the INPUT subdirectory. You may have requested amip2 data which
 211 !! has not been officially set up.
 212 !! See the section on DATA SETS to properly set the data up.
 213 !!
 214 !> @throws "FATAL: use_climo mismatch"
 215 !! The namelist variable date_out_of_range = 'fail' and the amip_interp_new
 216 !! argument use_climo = true.  This combination is not allowed.
 217 !!
 218 !> @throws "FATAL: use_annual(climo) mismatch"
 219 !! The namelist variable date_out_of_range = 'fail' and the amip_interp_new
 220 !! argument use_annual = true.  This combination is not allowed.
 221 !!
 222 !> @ingroup amip_interp_mod
 223 interface amip_interp_new
 224    module procedure amip_interp_new_1d
 225    module procedure amip_interp_new_2d
 226 end interface
 227
 228
 229 !----- public data type ------
 230
 231 !> @brief Contains information needed by the interpolation module (exchange_mod) and buffers data.
 232 !> @ingroup amip_interp_mod
 233 type amip_interp_type
 234    private
 235    type (horiz_interp_type) :: Hintrp, Hintrp2 ! add by JHC
 236    real, pointer            ::    data1(:,:) =>NULL(), &
 237                                   data2(:,:) =>NULL()
 238    type (date_type)         ::    Date1,       Date2
 239    logical                  :: use_climo, use_annual
 240    logical                  :: I_am_initialized=.false.
 241 end type
 242
 243 !> @addtogroup amip_interp_mod
 244 !> @{
 245 !-----------------------------------------------------------------------
 246 !  ---- resolution/grid variables ----
 247
 248    integer :: mobs, nobs
 249    real, allocatable :: lon_bnd(:), lat_bnd(:)
 250
 251 !  ---- global unit & date ----
 252
 253    integer, parameter :: maxc = 128
 254    integer :: unit
 255    character(len=maxc) :: file_name_sst, file_name_ice
 256    type(FmsNetcdfFile_t), target :: fileobj_sst, fileobj_ice
 257
 258    type (date_type) :: Curr_date = date_type( -99, -99, -99 )
 259    type (date_type) :: Date_end  = date_type( -99, -99, -99 )
 260
 261    real             :: tice_crit_k
 262    integer(I2_KIND) ::  ice_crit
 263
 264    logical :: module_is_initialized = .false.
 265
 266 !-----------------------------------------------------------------------
 267 !---- namelist ----
 268
 269  character(len=24) :: data_set = 'amip1' !< use 'amip1', 'amip2', 'reynolds_eof'
 270                                          !! 'reynolds_oi', 'hurrell', or 'daily',
 271                                          !! when "use_daily=.T."
 272                                          ! add by JHC
 273
 274  character(len=16) :: date_out_of_range = 'fail' !< use 'fail', 'initclimo', or 'climo'
 275
 276  real    :: tice_crit    = -1.80 !<  in degC or degK
 277  integer :: verbose      = 0     !<  0 <= verbose <= 3
 278
 279  logical :: use_zonal    = .false. !< parameters for prescribed zonal sst option
 280  real :: teq  = 305. !< parameters for prescribed zonal sst option
 281  real :: tdif = 50. !< parameters for prescribed zonal sst option
 282  real :: tann = 20. !< parameters for prescribed zonal sst option
 283  real :: tlag = 0.875 !< parameters for prescribed zonal sst option
 284
 285
 286  integer :: amip_date(3)=(/-1,-1,-1/) !< amip date for repeating single day (rsd) option
 287
 288  real :: sst_pert = 0. !< global temperature perturbation used for sensitivity experiments
 289
 290  character(len=6) :: sst_pert_type = 'fixed'  !< use 'random' or 'fixed'
 291  logical :: do_sst_pert = .false.
 292  logical :: use_daily = .false. !< if '.true.', give 'data_set = 'daily''
 293
 294  logical :: use_ncep_sst = .false. !< SJL: During nudging:   use_ncep_sst = .T.;  no_anom_sst = .T.
 295                                    !!      during forecast:  use_ncep_sst = .T.;  no_anom_sst = .F.
 296  logical ::  no_anom_sst = .true.  !< SJL: During nudging:   use_ncep_sst = .T.;  no_anom_sst = .T.
 297                                    !!      during forecast:  use_ncep_sst = .T.;  no_anom_sst = .F.
 298  logical :: use_ncep_ice = .false. !< For seasonal forecast: use_ncep_ice = .F.
 299  logical :: interp_oi_sst = .false. !< changed to false for regular runs
 300  logical :: use_mpp_io = .false. !< Set to .true. to use mpp_io, otherwise fms2io is used
 301
 302 !> @page amip_interp_nml @ref amip_interp_mod Namelist
 303 !!
 304 !> @var character(len=24) data_set
 305 !! Name/type of SST data that will be used.
 306 !!        Possible values (case-insensitive) are:
 307 !!                          1) amip1
 308 !!                          2) reynolds_eof
 309 !!                          3) reynolds_oi
 310 !!        See the @ref amip_interp_oi page for more information
 311 !! @var character(len=16) date_out_of_range
 312 !!     Controls the use of climatological monthly mean data when
 313 !!     the requested date falls outside the range of the data set.<BR/>
 314 !!     Possible values are:
 315 !!     <PRE>
 316 !!   fail      - program will fail if requested date is prior
 317 !!               to or after the data set period.
 318 !!   initclimo - program uses climatological requested data is
 319 !!               prior to data set period and will fail if
 320 !!               requested date is after data set period.
 321 !!   climo     - program uses climatological data anytime.
 322 !!    </PRE>
 323 !! @var real tice_crit
 324 !!     Freezing point of sea water in degC or degK. Defaults to -1.80
 325 !! @var integer verbose
 326 !!     Controls printed output, 0 <= verbose <= 3, default=0
 327 !!     additional parameters for controlling zonal prescribed sst ----
 328 !!     these parameters only have an effect when use_zonal=.true. ----
 329 !! @var logical use_zonal
 330 !!     Flag to selected zonal sst or data set. Default=.false.
 331 !! @var real teq
 332 !!     sst at the equator. Default=305
 333 !! @var real tdif
 334 !!     Equator to pole sst difference. Default=50
 335 !! @var real tann
 336 !!     Amplitude of annual cycle. Default=20
 337 !! @var real tlag
 338 !!     Offset for time of year (for annual cycle). Default=0.875
 339 !! @var integer amip_date
 340 !!     Single calendar date in integer "(year,month,day)" format
 341 !!     that is used only if set with year>0, month>0, day>0.
 342 !!     If used, model calendar date is replaced by this date,
 343 !!     but model time of day is still used to determine ice/sst.
 344 !!     Used for repeating-single-day (rsd) experiments.
 345 !!     Default=/-1,-1,-1/
 346 !! @var real sst_pert
 347 !!     Temperature perturbation in degrees Kelvin added onto the SST.
 348 !!                The perturbation is globally-uniform (even near sea-ice).
 349 !!                It is only used when abs(sst_pert) > 1.e-4.  SST perturbation runs
 350 !!                may be useful in accessing model sensitivities.
 351 !!     Default=0.
 352  namelist /amip_interp_nml/ use_ncep_sst, no_anom_sst, use_ncep_ice,  tice_crit, &
 353                             interp_oi_sst, data_set, date_out_of_range,          &
 354                             use_zonal, teq, tdif, tann, tlag, amip_date,         &
 355                             ! add by JHC
 356                             sst_pert, sst_pert_type, do_sst_pert,                &
 357                             use_daily,                                           &
 358                             ! end add by JHC
 359                             verbose, i_sst, j_sst, forecast_mode,                &
 360                             use_mpp_io
 361
 362 !-----------------------------------------------------------------------
 363
 364 contains
 365
 366 ! modified by JHC
 367 !> Retrieve sea surface temperature data and interpolated grid
 368 subroutine get_amip_sst (Time, Interp, sst, err_msg, lon_model, lat_model)
 369
 370    type (time_type),         intent(in)    :: Time !< Time to interpolate
 371    type (amip_interp_type),  intent(inout) :: Interp !< Holds data for interpolation
 372    real,                     intent(out)   ::  sst(:,:) !< Sea surface temperature data
 373    character(len=*), optional, intent(out) :: err_msg !< Holds error message string if present
 374
 375    real, dimension(mobs,nobs) :: sice
 376
 377     integer :: year1, year2, month1, month2
 378     real    :: fmonth
 379     type (date_type) :: Date1, Date2, Udate1, Udate2
 380
 381     type(time_type) :: Amip_Time
 382     integer :: tod(3),dum(3)
 383
 384 ! add by JHC
 385     real,    intent(in), dimension(:,:), optional :: lon_model, lat_model
 386     real :: pert
 387     integer :: i, j, mobs_sst, nobs_sst
 388     integer :: jhctod(6)
 389     type (time_type) :: Udate
 390     character(len=4) :: yyyy
 391     integer :: nrecords, ierr, k, yr, mo, dy
 392     integer, dimension(:), allocatable :: ryr, rmo, rdy
 393     character(len=30) :: time_unit
 394     real, dimension(:), allocatable :: timeval
 395     character(len=maxc) :: ncfilename
 396     type(FmsNetcdfFile_t) :: fileobj
 397     logical :: the_file_exists
 398 ! end add by JHC
 399     logical, parameter :: DEBUG = .false. !> switch for debugging output
 400     !> These are fms_io specific
 401     integer :: unit
 402
 403     if(present(err_msg)) err_msg = ''
 404     if(.not.Interp%I_am_initialized) then
 405       if(fms_error_handler('get_amip_sst','The amip_interp_type variable is not initialized',err_msg)) return
 406     endif
 407
 408 !-----------------------------------------------------------------------
 409 !----- compute zonally symetric sst ---------------
 410
 411     if ( use_ncep_sst .and. forecast_mode ) no_anom_sst = .false.
 412
 413     if (all(amip_date>0)) then
 414        call get_date(Time,dum(1),dum(2),dum(3),tod(1),tod(2),tod(3))
 415        Amip_Time = set_date(amip_date(1),amip_date(2),amip_date(3),tod(1),tod(2),tod(3))
 416     else
 417        Amip_Time = Time
 418     endif
 419
 420 ! add by JHC
 421 if ( .not.use_daily ) then
 422 ! end add by JHC
 423
 424    if ( .not. allocated(temp1) ) allocate (temp1(mobs,nobs))
 425    if ( .not. allocated(temp2) ) allocate (temp2(mobs,nobs))
 426
 427    if (use_zonal) then
 428       call zonal_sst (Amip_Time, sice, temp1)
 429       call horiz_interp ( Interp%Hintrp, temp1, sst )
 430    else
 431
 432 !-----------------------------------------------------------------------
 433 !---------- get new observed sea surface temperature -------------------
 434
 435 ! ---- time interpolation for months -----
 436      call time_interp (Amip_Time, fmonth, year1, year2, month1, month2)
 437 ! ---- force climatology ----
 438      if (Interp % use_climo) then
 439          year1=0; year2=0
 440      endif
 441      if (Interp % use_annual) then
 442           year1=0;  year2=0
 443          month1=0; month2=0
 444      endif
 445 ! ---------------------------
 446
 447      Date1 = date_type( year1, month1, 0 )
 448      Date2 = date_type( year2, month2, 0 )
 449
 450 !  -- open/rewind file --
 451      unit = -1
 452 !-----------------------------------------------------------------------
 453
 454
 455       if (Date1 /= Interp % Date1) then
 456 !       ---- use Date2 for Date1 ----
 457           if (Date1 == Interp % Date2) then
 458               Interp % Date1 = Interp % Date2
 459               Interp % data1 = Interp % data2
 460               temp1(:,:) = temp2(:,:)   ! SJL BUG fix: June 24, 2011
 461           else
 462               call read_record ('sst', Date1, Udate1, temp1)
 463               if ( use_ncep_sst .and. (.not. no_anom_sst) ) then
 464                    temp1(:,:) = temp1(:,:) + sst_anom(:,:)
 465               endif
 466               call horiz_interp ( Interp%Hintrp, temp1, Interp%data1 )
 467               call clip_data ('sst', Interp%data1)
 468              Interp % Date1 = Date1
 469           endif
 470       endif
 471
 472 !-----------------------------------------------------------------------
 473
 474       if (Date2 /= Interp % Date2) then
 475           call read_record ('sst', Date2, Udate2, temp2)
 476           if ( use_ncep_sst .and. (.not. no_anom_sst) ) then
 477                temp2(:,:) = temp2(:,:) + sst_anom(:,:)
 478           endif
 479           call horiz_interp ( Interp%Hintrp, temp2, Interp%data2 )
 480           call clip_data ('sst', Interp%data2)
 481           Interp % Date2 = Date2
 482       endif
 483
 484 !-----------------------------------------------------------------------
 485 !---------- time interpolation (between months) of sst's ---------------
 486 !-----------------------------------------------------------------------
 487     sst = Interp % data1 + fmonth * (Interp % data2 - Interp % data1)
 488
 489 !-------------------------------------------------------------------------------
 490 ! SJL mods for NWP and TCSF ---
 491 !      Nudging runs: (Note: NCEP SST updated only every 6-hr)
 492 !      Compute SST anomaly from global SST datasets for subsequent forecast runs
 493 !-------------------------------------------------------------------------------
 494     if ( use_ncep_sst .and. no_anom_sst ) then
 495          sst_anom(:,:) = sst_ncep(:,:) - (temp1(:,:) + fmonth*(temp2(:,:) - temp1(:,:)) )
 496          call horiz_interp ( Interp%Hintrp, sst_ncep, sst )
 497          call clip_data ('sst', sst)
 498     endif
 499
 500 !! DEBUG CODE
 501     if (DEBUG) then
 502           call get_date(Amip_Time,jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6))
 503           if (mpp_pe() == 0) then
 504              write (*,200) 'JHC: use_daily = F, AMIP_Time: ',jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5), &
 505                    & jhctod(6)
 506              write (*,300) 'JHC: use_daily = F, interped SST: ', sst(1,1),sst(5,5),sst(10,10)
 507           endif
 508     endif
 509
 510
 511   endif
 512
 513 ! add by JHC
 514 else
 515     call get_date(Amip_Time,jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6))
 516      if (mpp_pe() == mpp_root_pe()) write(*,200) 'amip_interp_mod: use_daily = T, Amip_Time = ',jhctod(1), &
 517         & jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6)
 518
 519     yr = jhctod(1); mo = jhctod(2); dy = jhctod(3)
 520
 521     write (yyyy,'(i4)') jhctod(1)
 522
 523     file_name_sst = 'INPUT/' // 'sst.day.mean.'//yyyy//'.v2.nc'
 524     ncfilename = trim(file_name_sst)
 525     time_unit = 'days since 1978-01-01 00:00:00'
 526
 527     mobs_sst = 1440;  nobs_sst = 720
 528
 529     call set_sst_grid_edges_daily(mobs_sst, nobs_sst)
 530     call horiz_interp_new ( Interp%Hintrp2, lon_bnd, lat_bnd, &
 531                              lon_model, lat_model, interp_method="bilinear" )
 532
 533     the_file_exists = fms2_io_file_exists(ncfilename)
 534
 535     if ( (.NOT. the_file_exists)  ) then
 536         call mpp_error ('amip_interp_mod', &
 537              'cannot find daily SST input data file: '//trim(ncfilename), NOTE)
 538     else
 539         if (mpp_pe() == mpp_root_pe()) call mpp_error ('amip_interp_mod', &
 540              'Reading NetCDF formatted daily SST from: '//trim(ncfilename), NOTE)
 541
 542             if(.not. open_file(fileobj, trim(ncfilename), 'read')) &
 543                 call error_mesg ('get_amip_sst', 'Error in opening file '//trim(ncfilename), FATAL)
 544
 545             call get_dimension_size(fileobj, 'TIME', nrecords)
 546             if (nrecords < 1) call mpp_error('amip_interp_mod', &
 547                            'Invalid number of SST records in daily SST data file: '//trim(ncfilename), FATAL)
 548             allocate(timeval(nrecords), ryr(nrecords), rmo(nrecords), rdy(nrecords))
 549             call fms2_io_read_data(fileobj, 'TIME', timeval)
 550 !!! DEBUG CODE
 551         if(DEBUG) then
 552           if (mpp_pe() == 0) then
 553              print *, 'JHC: nrecords = ', nrecords
 554              print *, 'JHC: TIME = ', timeval
 555           endif
 556         endif
 557
 558         ierr = 1
 559         do k = 1, nrecords
 560
 561           Udate = get_cal_time (timeval(k), time_unit, 'julian')
 562           call get_date(Udate,jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6))
 563           ryr(k) = jhctod(1); rmo(k) = jhctod(2); rdy(k) = jhctod(3)
 564
 565           if ( yr == ryr(k) .and. mo == rmo(k) .and. dy == rdy (k) ) ierr = 0
 566           if (ierr==0) exit
 567
 568         enddo
 569
 570         if(DEBUG) then
 571           if (mpp_pe() == 0) then
 572             print *, 'JHC: k =', k
 573             print *, 'JHC: ryr(k) rmo(k) rdy(k)',ryr(k), rmo(k), rdy(k)
 574             print *, 'JHC:  yr     mo     dy   ',yr, mo, dy
 575           endif
 576         endif
 577
 578         if (ierr .ne. 0) call mpp_error('amip_interp_mod', &
 579                          'Model time is out of range not in SST data: '//trim(ncfilename), FATAL)
 580     endif ! if(file_exist(ncfilename))
 581
 582
 583    !---- read NETCDF data ----
 584      if ( .not. allocated(tempamip) ) allocate (tempamip(mobs_sst,nobs_sst))
 585
 586      if (the_file_exists) then
 587           call fms2_io_read_data(fileobj, 'SST', tempamip, unlim_dim_level=k)
 588           call close_file(fileobj)
 589           tempamip = tempamip + TFREEZE
 590
 591 !!! DEBUG CODE
 592           if(DEBUG) then
 593             if (mpp_pe() == 0) then
 594               print*, 'JHC: TFREEZE = ', TFREEZE
 595               print*, lbound(sst)
 596               print*, ubound(sst)
 597               print*, lbound(tempamip)
 598               print*, ubound(tempamip)
 599               write(*,300) 'JHC: tempamip : ', tempamip(100,100), tempamip(200,200), tempamip(300,300)
 600             endif
 601           endif
 602
 603           call horiz_interp ( Interp%Hintrp2, tempamip, sst )
 604           call clip_data ('sst', sst)
 605
 606      endif
 607
 608     if(DEBUG) then
 609       if (mpp_pe() == 400) then
 610         write(*,300)'JHC: use_daily = T, daily SST: ', sst(1,1),sst(5,5),sst(10,10)
 611         print *,'JHC: use_daily = T, daily SST: ', sst
 612       endif
 613     endif
 614
 615 200 format(a35, 6(i5,1x))
 616 300 format(a35, 3(f7.3,2x))
 617
 618 endif
 619 ! end add by JHC
 620
 621 ! add by JHC: add on non-zero sea surface temperature perturbation (namelist option)
 622 !             This perturbation may be useful in accessing model sensitivities
 623
 624  if ( do_sst_pert ) then
 625
 626       if ( trim(sst_pert_type) == 'fixed' ) then
 627           sst = sst + sst_pert
 628       else if ( trim(sst_pert_type) == 'random' ) then
 629           call random_seed()
 630
 631        if(DEBUG) then
 632          if (mpp_pe() == 0) then
 633              print*, 'mobs = ', mobs
 634              print*, 'nobs = ', nobs
 635              print*, lbound(sst)
 636              print*, ubound(sst)
 637           endif
 638        endif
 639
 640           do i = 1, size(sst,1)
 641           do j = 1, size(sst,2)
 642              call random_number(pert)
 643              sst (i,j) = sst (i,j) + sst_pert*((pert-0.5)*2)
 644           end do
 645           end do
 646       endif
 647
 648   endif
 649 ! end add by JHC
 650
 651 !-----------------------------------------------------------------------
 652
 653  end subroutine get_amip_sst
 654
 655 !> AMIP interpolation for ice
 656 subroutine get_amip_ice (Time, Interp, ice, err_msg)
 657
 658    type (time_type),         intent(in)    :: Time !< Time to interpolate
 659    type (amip_interp_type),  intent(inout) :: Interp !< Holds data for interpolation
 660    real,                     intent(out)   :: ice(:,:) !< ice data
 661    character(len=*), optional, intent(out) :: err_msg !< Holds error message string if present
 662
 663     real, dimension(mobs,nobs) :: sice, temp
 664
 665     integer :: year1, year2, month1, month2
 666     real    :: fmonth
 667     type (date_type) :: Date1, Date2, Udate1, Udate2
 668
 669     type(time_type) :: Amip_Time
 670     integer :: tod(3),dum(3)
 671
 672     if(present(err_msg)) err_msg = ''
 673     if(.not.Interp%I_am_initialized) then
 674       if(fms_error_handler('get_amip_ice','The amip_interp_type variable is not initialized',err_msg)) return
 675     endif
 676
 677 !-----------------------------------------------------------------------
 678 !----- compute zonally symetric sst ---------------
 679
 680
 681     if (any(amip_date>0)) then
 682
 683        call get_date(Time,dum(1),dum(2),dum(3),tod(1),tod(2),tod(3))
 684
 685        Amip_Time = set_date(amip_date(1),amip_date(2),amip_date(3),tod(1),tod(2),tod(3))
 686
 687     else
 688
 689        Amip_Time = Time
 690
 691     endif
 692
 693
 694 if (use_zonal) then
 695    call zonal_sst (Amip_Time, sice, temp)
 696    call horiz_interp ( Interp%Hintrp, sice, ice )
 697 else
 698
 699 !-----------------------------------------------------------------------
 700 !---------- get new observed sea surface temperature -------------------
 701
 702 ! ---- time interpolation for months -----
 703
 704    call time_interp (Amip_Time, fmonth, year1, year2, month1, month2)
 705
 706 ! ---- force climatology ----
 707    if (Interp % use_climo) then
 708        year1=0; year2=0
 709    endif
 710    if (Interp % use_annual) then
 711         year1=0;  year2=0
 712        month1=0; month2=0
 713    endif
 714 ! ---------------------------
 715
 716    Date1 = date_type( year1, month1, 0 )
 717    Date2 = date_type( year2, month2, 0 )
 718
 719    unit = -1
 720 !-----------------------------------------------------------------------
 721
 722     if (Date1 /= Interp % Date1) then
 723 !       ---- use Date2 for Date1 ----
 724         if (Date1 == Interp % Date2) then
 725             Interp % Date1 = Interp % Date2
 726             Interp % data1 = Interp % data2
 727         else
 728 !-- SJL -------------------------------------------------------------
 729 ! Can NOT use ncep_sst to determine sea_ice For seasonal forecast
 730 ! Use climo sea ice for seasonal runs
 731             if ( use_ncep_sst .and. use_ncep_ice ) then
 732                where ( sst_ncep <= (TFREEZE+tice_crit) )
 733                    sice = 1.
 734                elsewhere
 735                    sice = 0.
 736                endwhere
 737             else
 738                call read_record ('ice', Date1, Udate1, sice)
 739             endif
 740 !--------------------------------------------------------------------
 741             call horiz_interp ( Interp%Hintrp, sice, Interp%data1 )
 742             call clip_data ('ice', Interp%data1)
 743             Interp % Date1 = Date1
 744         endif
 745     endif
 746
 747 !-----------------------------------------------------------------------
 748
 749     if (Date2 /= Interp % Date2) then
 750
 751 !-- SJL -------------------------------------------------------------
 752             if ( use_ncep_sst .and. use_ncep_ice ) then
 753                where ( sst_ncep <= (TFREEZE+tice_crit) )
 754                    sice = 1.
 755                elsewhere
 756                    sice = 0.
 757                endwhere
 758             else
 759                call read_record ('ice', Date2, Udate2, sice)
 760             endif
 761 !--------------------------------------------------------------------
 762         call horiz_interp ( Interp%Hintrp, sice, Interp%data2 )
 763         call clip_data ('ice', Interp%data2)
 764         Interp % Date2 = Date2
 765
 766     endif
 767
 768 !-----------------------------------------------------------------------
 769 !---------- time interpolation (between months) ------------------------
 770 !-----------------------------------------------------------------------
 771
 772    ice = Interp % data1 + fmonth * (Interp % data2 - Interp % data1)
 773
 774 endif
 775
 776 !-----------------------------------------------------------------------
 777
 778  end subroutine get_amip_ice
 779
 780 !#######################################################################
 781
 782  !> @return A newly created @ref amip_interp_type
 783  function amip_interp_new_1d ( lon , lat , mask , use_climo, use_annual, &
 784                                 interp_method ) result (Interp)
 785
 786  real,    intent(in), dimension(:)   :: lon, lat
 787  logical, intent(in), dimension(:,:) :: mask
 788  character(len=*), intent(in), optional       :: interp_method
 789  logical, intent(in), optional       :: use_climo, use_annual
 790
 791    type (amip_interp_type) :: Interp
 792
 793    if(.not.module_is_initialized) call amip_interp_init
 794
 795    Interp % use_climo  = .false.
 796    if (present(use_climo)) Interp % use_climo  = use_climo
 797    Interp % use_annual = .false.
 798    if (present(use_annual)) Interp % use_annual  = use_annual
 799
 800    if ( date_out_of_range == 'fail' .and. Interp%use_climo ) &
 801       call error_mesg ('amip_interp_new_1d', 'use_climo mismatch', FATAL)
 802
 803    if ( date_out_of_range == 'fail' .and. Interp%use_annual ) &
 804       call error_mesg ('amip_interp_new_1d', 'use_annual(climo) mismatch', FATAL)
 805
 806    Interp % Date1 = date_type( -99, -99, -99 )
 807    Interp % Date2 = date_type( -99, -99, -99 )
 808
 809 !-----------------------------------------------------------------------
 810 !   ---- initialization of horizontal interpolation ----
 811
 812     call horiz_interp_new ( Interp%Hintrp, lon_bnd, lat_bnd, &
 813                              lon, lat, interp_method= interp_method )
 814
 815     allocate ( Interp % data1 (size(lon(:))-1,size(lat(:))-1), &
 816                Interp % data2 (size(lon(:))-1,size(lat(:))-1)  )
 817
 818     Interp%I_am_initialized = .true.
 819
 820    end function amip_interp_new_1d
 821
 822  !> @return A newly created @ref amip_interp_type
 823  function amip_interp_new_2d ( lon , lat , mask , use_climo, use_annual, &
 824                                 interp_method ) result (Interp)
 825
 826  real,    intent(in), dimension(:,:)   :: lon, lat
 827  logical, intent(in), dimension(:,:) :: mask
 828  character(len=*), intent(in), optional :: interp_method
 829  logical, intent(in), optional       :: use_climo, use_annual
 830
 831    type (amip_interp_type) :: Interp
 832
 833    if(.not.module_is_initialized) call amip_interp_init
 834
 835    Interp % use_climo  = .false.
 836    if (present(use_climo)) Interp % use_climo  = use_climo
 837    Interp % use_annual = .false.
 838    if (present(use_annual)) Interp % use_annual  = use_annual
 839
 840    if ( date_out_of_range == 'fail' .and. Interp%use_climo ) &
 841       call error_mesg ('amip_interp_new_2d', 'use_climo mismatch', FATAL)
 842
 843    if ( date_out_of_range == 'fail' .and. Interp%use_annual ) &
 844       call error_mesg ('amip_interp_new_2d', 'use_annual(climo) mismatch', FATAL)
 845
 846    Interp % Date1 = date_type( -99, -99, -99 )
 847    Interp % Date2 = date_type( -99, -99, -99 )
 848
 849 !-----------------------------------------------------------------------
 850 !   ---- initialization of horizontal interpolation ----
 851
 852    call horiz_interp_new ( Interp%Hintrp, lon_bnd, lat_bnd, &
 853                            lon, lat, interp_method = interp_method)
 854
 855    allocate ( Interp % data1 (size(lon,1),size(lat,2)), &
 856               Interp % data2 (size(lon,1),size(lat,2)))
 857
 858    Interp%I_am_initialized = .true.
 859
 860    end function amip_interp_new_2d
 861
 862 !#######################################################################
 863
 864  !> initialize @ref amip_interp_mod for use
 865  subroutine amip_interp_init()
 866
 867    integer :: unit,io,ierr
 868
 869 !-----------------------------------------------------------------------
 870
 871     call horiz_interp_init
 872
 873 !   ---- read namelist ----
 874
 875     read (input_nml_file, amip_interp_nml, iostat=io)
 876     ierr = check_nml_error(io,'amip_interp_nml')
 877
 878 !  ----- write namelist/version info -----
 879     call write_version_number("AMIP_INTERP_MOD", version)
 880
 881     unit = stdlog ( )
 882     if (mpp_pe() == 0) then
 883         write (unit,nml=amip_interp_nml)
 884     endif
 885
 886     if (use_mpp_io) then
 887             !! USE_MPP_IO_WARNING
 888             call mpp_error ('amip_interp_mod', &
 889              'MPP_IO is no longer supported.  Please remove use_mpp_io from amip_interp_nml',&
 890               FATAL)
 891     endif
 892     if ( .not. use_ncep_sst ) interp_oi_sst = .false.
 893
 894 !   ---- freezing point of sea water in deg K ---
 895
 896     tice_crit_k = tice_crit
 897     if ( tice_crit_k < 200. ) tice_crit_k = tice_crit_k + TFREEZE
 898     ice_crit = nint((tice_crit_k-TFREEZE)*100., I2_KIND)
 899
 900 !   ---- set up file dependent variable ----
 901 !   ----   global file name   ----
 902 !   ----   grid box edges     ----
 903 !   ---- initialize zero size grid if not pe 0 ------
 904
 905     if (lowercase(trim(data_set)) == 'amip1') then
 906         file_name_sst = 'INPUT/' // 'amip1_sst.data'
 907         file_name_ice = 'INPUT/' // 'amip1_sst.data'
 908         mobs = 180;  nobs = 91
 909         call set_sst_grid_edges_amip1
 910         if (mpp_pe() == 0) &
 911         call error_mesg ('amip_interp_init', 'using AMIP 1 sst', NOTE)
 912         Date_end = date_type( 1989, 1, 0 )
 913     else if (lowercase(trim(data_set)) == 'amip2') then
 914         file_name_sst = 'INPUT/' // 'amip2_sst.data'
 915         file_name_ice = 'INPUT/' // 'amip2_ice.data'
 916         mobs = 360;  nobs = 180
 917         call set_sst_grid_edges_oi
 918 !       --- specfied min for amip2 ---
 919         tice_crit_k = 271.38
 920         if (mpp_pe() == 0) &
 921         call error_mesg ('amip_interp_init', 'using AMIP 2 sst', NOTE)
 922         Date_end = date_type( 1996, 3, 0 )
 923     else if (lowercase(trim(data_set)) == 'hurrell') then
 924         file_name_sst = 'INPUT/' // 'hurrell_sst.data'
 925         file_name_ice = 'INPUT/' // 'hurrell_ice.data'
 926         mobs = 360;  nobs = 180
 927         call set_sst_grid_edges_oi
 928 !       --- specfied min for hurrell ---
 929         tice_crit_k = 271.38
 930         if (mpp_pe() == 0) &
 931         call error_mesg ('amip_interp_init', 'using HURRELL sst', NOTE)
 932         Date_end = date_type( 2011, 8, 16 ) ! updated by JHC
 933 ! add by JHC
 934     else if (lowercase(trim(data_set)) == 'daily') then
 935         file_name_sst = 'INPUT/' // 'hurrell_sst.data'
 936         file_name_ice = 'INPUT/' // 'hurrell_ice.data'
 937         mobs = 360;  nobs = 180
 938         call set_sst_grid_edges_oi
 939         if (mpp_pe() == 0) &
 940         call error_mesg ('amip_interp_init', 'using AVHRR daily sst', NOTE)
 941         Date_end = date_type( 2011, 8, 16 )
 942 ! end add by JHC
 943     else if (lowercase(trim(data_set)) == 'reynolds_eof') then
 944         file_name_sst = 'INPUT/' // 'reynolds_sst.data'
 945         file_name_ice = 'INPUT/' // 'reynolds_sst.data'
 946         mobs = 180;  nobs = 90
 947         call set_sst_grid_edges_oi
 948         if (mpp_pe() == 0) &
 949         call error_mesg ('amip_interp_init',  &
 950              'using NCEP Reynolds Historical Reconstructed SST', NOTE)
 951         Date_end = date_type( 1998, 12, 0 )
 952     else if (lowercase(trim(data_set)) == 'reynolds_oi') then
 953         file_name_sst = 'INPUT/' // 'reyoi_sst.data'
 954         file_name_ice = 'INPUT/' // 'reyoi_sst.data'
 955 !--- Added by SJL ----------------------------------------------
 956         if ( use_ncep_sst ) then
 957              mobs = i_sst;  nobs = j_sst
 958             if (.not. allocated (sst_ncep)) then
 959                 allocate (sst_ncep(i_sst,j_sst))
 960                 sst_ncep(:,:) = big_number
 961             endif
 962             if (.not. allocated (sst_anom)) then
 963                 allocate (sst_anom(i_sst,j_sst))
 964                 sst_anom(:,:) = big_number
 965             endif
 966         else
 967              mobs = 360;    nobs = 180
 968         endif
 969 !--- Added by SJL ----------------------------------------------
 970         call set_sst_grid_edges_oi
 971         if (mpp_pe() == 0) &
 972         call error_mesg ('amip_interp_init', 'using Reynolds OI SST', &
 973                                                                 NOTE)
 974         Date_end = date_type( 1999, 1, 0 )
 975     else
 976         call error_mesg ('amip_interp_init', 'the value of the &
 977         &namelist parameter DATA_SET being used is not allowed', FATAL)
 978     endif
 979
 980     if (verbose > 1 .and. mpp_pe() == 0) &
 981               print *, 'ice_crit,tice_crit_k=',ice_crit,tice_crit_k
 982
 983 !  --- check existence of sst data file ??? ---
 984     file_name_sst = trim(file_name_sst)//'.nc'
 985     file_name_ice = trim(file_name_ice)//'.nc'
 986
 987     if (.not. fms2_io_file_exists(trim(file_name_sst)) ) then
 988         call error_mesg ('amip_interp_init', &
 989              'file '//trim(file_name_sst)//' does not exist', FATAL)
 990     endif
 991     if (.not. fms2_io_file_exists(trim(file_name_ice)) ) then
 992         call error_mesg ('amip_interp_init', &
 993              'file '//trim(file_name_ice)//' does not exist', FATAL)
 994     endif
 995
 996     if (.not. open_file(fileobj_sst, trim(file_name_sst), 'read')) &
 997         call error_mesg ('amip_interp_init', 'Error in opening file '//trim(file_name_sst), FATAL)
 998     if (.not. open_file(fileobj_ice, trim(file_name_ice), 'read')) &
 999         call error_mesg ('amip_interp_init', 'Error in opening file '//trim(file_name_ice), FATAL)
1000     module_is_initialized = .true.
1001
1002  end subroutine amip_interp_init
1003
1004 !#######################################################################
1005
1006 !> Frees data associated with a amip_interp_type variable. Should be used for any
1007 !! variables initialized via @ref amip_interp_new.
1008 !> @param[inout] Interp A defined data type variable initialized by amip_interp_new and used
1009 !! when calling get_amip_sst and get_amip_ice.
1010    subroutine amip_interp_del (Interp)
1011    type (amip_interp_type), intent(inout) :: Interp
1012      if(associated(Interp%data1)) deallocate(Interp%data1)
1013      if(associated(Interp%data2)) deallocate(Interp%data2)
1014      if(allocated(lon_bnd))       deallocate(lon_bnd)
1015      if(allocated(lat_bnd))       deallocate(lat_bnd)
1016      call horiz_interp_del ( Interp%Hintrp )
1017
1018      Interp%I_am_initialized = .false.
1019
1020    end subroutine amip_interp_del
1021
1022 !#######################################################################
1023
1024    subroutine set_sst_grid_edges_amip1
1025
1026    integer :: i, j
1027    real    :: hpie, dlon, dlat, wb, sb
1028
1029       allocate ( lon_bnd(mobs+1), lat_bnd(nobs+1) )
1030
1031 ! ---- compute grid edges (do only once) -----
1032
1033       hpie = 0.5*pi
1034
1035       dlon = 4.*hpie/float(mobs);  wb = -0.5*dlon
1036       do i = 1, mobs+1
1037           lon_bnd(i) = wb + dlon * float(i-1)
1038       enddo
1039           lon_bnd(mobs+1) = lon_bnd(1) + 4.*hpie
1040
1041       dlat = 2.*hpie/float(nobs-1);  sb = -hpie + 0.5*dlat
1042       lat_bnd(1) = -hpie;  lat_bnd(nobs+1) = hpie
1043       do j = 2, nobs
1044           lat_bnd(j) = sb + dlat * float(j-2)
1045       enddo
1046
1047    end subroutine set_sst_grid_edges_amip1
1048
1049 !#######################################################################
1050    subroutine set_sst_grid_edges_oi
1051
1052    integer :: i, j
1053    real    :: hpie, dlon, dlat, wb, sb
1054
1055 ! add by JHC
1056       if(allocated(lon_bnd))       deallocate(lon_bnd)
1057       if(allocated(lat_bnd))       deallocate(lat_bnd)
1058 ! end add by JHC
1059       allocate ( lon_bnd(mobs+1), lat_bnd(nobs+1) )
1060
1061 ! ---- compute grid edges (do only once) -----
1062
1063       hpie = 0.5*pi
1064
1065       dlon = 4.*hpie/float(mobs);  wb = 0.0
1066           lon_bnd(1) = wb
1067       do i = 2, mobs+1
1068           lon_bnd(i) = wb + dlon * float(i-1)
1069       enddo
1070           lon_bnd(mobs+1) = lon_bnd(1) + 4.*hpie
1071
1072       dlat = 2.*hpie/float(nobs);  sb = -hpie
1073       lat_bnd(1) = sb;  lat_bnd(nobs+1) = hpie
1074       do j = 2, nobs
1075           lat_bnd(j) = sb + dlat * float(j-1)
1076       enddo
1077
1078    end subroutine set_sst_grid_edges_oi
1079 !#######################################################################
1080 ! add by JHC
1081    subroutine set_sst_grid_edges_daily(mobs_sst, nobs_sst)
1082
1083    integer :: i, j, mobs_sst, nobs_sst
1084    real    :: hpie, dlon, dlat, wb, sb
1085
1086       if(allocated(lon_bnd))       deallocate(lon_bnd)
1087       if(allocated(lat_bnd))       deallocate(lat_bnd)
1088       allocate ( lon_bnd(mobs_sst+1), lat_bnd(nobs_sst+1) )
1089
1090 ! ---- compute grid edges (do only once) -----
1091
1092       hpie = 0.5*pi
1093
1094       dlon = 4.*hpie/float(mobs_sst);  wb = 0.0
1095           lon_bnd(1) = wb
1096       do i = 2, mobs_sst+1
1097           lon_bnd(i) = wb + dlon * float(i-1)
1098       enddo
1099           lon_bnd(mobs_sst+1) = lon_bnd(1) + 4.*hpie
1100
1101       dlat = 2.*hpie/float(nobs_sst);  sb = -hpie
1102       lat_bnd(1) = sb;  lat_bnd(nobs_sst+1) = hpie
1103       do j = 2, nobs_sst
1104           lat_bnd(j) = sb + dlat * float(j-1)
1105       enddo
1106
1107    end subroutine set_sst_grid_edges_daily
1108 ! end add by JHC
1109 !#######################################################################
1110
1111
1112    subroutine a2a_bilinear(nx, ny, dat1, n1, n2, dat2)
1113    integer, intent(in):: nx, ny
1114    integer, intent(in):: n1, n2
1115    real, intent(in) :: dat1(nx,ny)
1116    real, intent(out):: dat2(n1,n2)      !> output interpolated data
1117
1118 ! local:
1119   real:: lon1(nx), lat1(ny)
1120   real:: lon2(n1), lat2(n2)
1121   real:: dx1, dy1, dx2, dy2
1122   real:: xc, yc
1123   real:: a1, b1, c1, c2, c3, c4
1124   integer i1, i2, jc, i0, j0, it, jt
1125   integer i,j
1126
1127
1128 !-----------------------------------------------------------
1129 ! * Interpolate from "FMS" 1x1 SST data grid to a finer grid
1130 !                     lon: 0.5, 1.5, ..., 359.5
1131 !                     lat: -89.5, -88.5, ... , 88.5, 89.5
1132 !-----------------------------------------------------------
1133
1134   dx1 = 360./real(nx) !> INput Grid
1135   dy1 = 180./real(ny) !> INput Grid
1136
1137   do i=1,nx
1138      lon1(i) = 0.5*dx1 + real(i-1)*dx1
1139   enddo
1140   do j=1,ny
1141      lat1(j) = -90. + 0.5*dy1 + real(j-1)*dy1
1142   enddo
1143
1144   dx2 = 360./real(n1) !> OutPut Grid:
1145   dy2 = 180./real(n2) !> OutPut Grid:
1146
1147   do i=1,n1
1148      lon2(i) = 0.5*dx2 + real(i-1)*dx2
1149   enddo
1150   do j=1,n2
1151      lat2(j) = -90. + 0.5*dy2 + real(j-1)*dy2
1152   enddo
1153
1154   jt = 1
1155   do 5000 j=1,n2
1156
1157      yc = lat2(j)
1158      if ( yc<lat1(1) ) then
1159             jc = 1
1160             b1 = 0.
1161      elseif ( yc>lat1(ny) ) then
1162             jc = ny-1
1163             b1 = 1.
1164      else
1165           do j0=jt,ny-1
1166           if ( yc>=lat1(j0) .and. yc<=lat1(j0+1) ) then
1167                jc = j0
1168                jt = j0
1169                b1 = (yc-lat1(jc)) / dy1
1170                go to 222
1171           endif
1172           enddo
1173      endif
1174 222  continue
1175
1176      it = 1
1177      do i=1,n1
1178         xc = lon2(i)
1179        if ( xc>lon1(nx) ) then
1180             i1 = nx;     i2 = 1
1181             a1 = (xc-lon1(nx)) / dx1
1182        elseif ( xc<lon1(1) ) then
1183             i1 = nx;     i2 = 1
1184             a1 = (xc+360.-lon1(nx)) / dx1
1185        else
1186             do i0=it,nx-1
1187             if ( xc>=lon1(i0) .and. xc<=lon1(i0+1) ) then
1188                i1 = i0;  i2 = i0+1
1189                it = i0
1190                a1 = (xc-lon1(i1)) / dx1
1191                go to 111
1192             endif
1193             enddo
1194        endif
1195 111    continue
1196
1197 ! Debug code:
1198        if ( a1<-0.001 .or. a1>1.001 .or.  b1<-0.001 .or. b1>1.001 ) then
1199             write(*,*) i,j,a1, b1
1200             call mpp_error(FATAL,'a2a bilinear interpolation')
1201        endif
1202
1203        c1 = (1.-a1) * (1.-b1)
1204        c2 =     a1  * (1.-b1)
1205        c3 =     a1  *     b1
1206        c4 = (1.-a1) *     b1
1207
1208 ! Bilinear interpolation:
1209        dat2(i,j) = c1*dat1(i1,jc) + c2*dat1(i2,jc) + c3*dat1(i2,jc+1) + c4*dat1(i1,jc+1)
1210
1211      enddo   !i-loop
1212
1213 5000 continue   ! j-loop
1214
1215    end subroutine a2a_bilinear
1216
1217 !#######################################################################
1218
1219 !> @brief Returns the size (i.e., number of longitude and latitude
1220 !!         points) of the observed data grid.
1221 !! @throws FATAL have not called amip_interp_new
1222 !!     Must call amip_interp_new before get_sst_grid_size.
1223    subroutine get_sst_grid_size (nlon, nlat)
1224
1225    integer, intent(out) :: nlon !> The number of longitude points (first dimension) in the
1226                                 !! observed data grid.  For AMIP 1 nlon = 180, and the Reynolds nlon = 360.
1227    integer, intent(out) :: nlat !> The number of latitude points (second dimension) in the
1228                                 !! observed data grid.  For AMIP 1 nlon = 91, and the Reynolds nlon = 180.
1229
1230       if ( .not.module_is_initialized ) call amip_interp_init
1231
1232       nlon = mobs;  nlat = nobs
1233
1234    end subroutine get_sst_grid_size
1235
1236 !#######################################################################
1237
1238 !> @brief Returns the grid box boundaries of the observed data grid.
1239 !!
1240 !! @throws FATAL, have not called amip_interp_new
1241 !!     Must call amip_interp_new before get_sst_grid_boundary.
1242 !!
1243 !! @throws FATAL, invalid argument dimensions
1244 !!     The size of the output argument arrays do not agree with
1245 !!     the size of the observed data. See the documentation for
1246 !!     interfaces get_sst_grid_size and get_sst_grid_boundary.
1247    subroutine get_sst_grid_boundary (blon, blat, mask)
1248
1249    real,    intent(out) :: blon(:) !> The grid box edges (in radians) for longitude points of the
1250                                    !! observed data grid. The size of this argument must be nlon+1.
1251    real,    intent(out) :: blat(:) !> The grid box edges (in radians) for latitude points of the
1252                                    !! observed data grid. The size of this argument must be nlat+1.
1253    logical, intent(out) :: mask(:,:)
1254
1255       if ( .not.module_is_initialized ) call amip_interp_init
1256
1257 ! ---- check size of argument(s) ----
1258
1259       if (size(blon(:)) /= mobs+1 .or. size(blat(:)) /= nobs+1)   &
1260       call error_mesg ('get_sst_grid_boundary in amip_interp_mod',  &
1261                        'invalid argument dimensions', FATAL)
1262
1263 ! ---- return grid box edges -----
1264
1265       blon = lon_bnd
1266       blat = lat_bnd
1267
1268 ! ---- masking (data exists at all points) ----
1269
1270       mask = .true.
1271
1272
1273    end subroutine get_sst_grid_boundary
1274
1275 !#######################################################################
1276
1277    subroutine read_record (type, Date, Adate, dat)
1278
1279      character(len=*), intent(in)  :: type
1280      type (date_type), intent(in)  :: Date
1281      type (date_type), intent(inout) :: Adate
1282      real,             intent(out) :: dat(mobs,nobs)
1283      real :: tmp_dat(360,180)
1284
1285      integer(I2_KIND) :: idat(mobs,nobs)
1286      integer :: nrecords, yr, mo, dy, ierr, k
1287      integer, dimension(:), allocatable :: ryr, rmo, rdy
1288      character(len=maxc) :: ncfilename, ncfieldname
1289      type(FmsNetcdfFile_t), pointer :: fileobj
1290
1291     !---- set file and field name for NETCDF data sets ----
1292
1293         ncfieldname = 'sst'
1294      if(type(1:3) == 'sst') then
1295         ncfilename = trim(file_name_sst)
1296         fileobj => fileobj_sst
1297      else if(type(1:3) == 'ice') then
1298         ncfilename = trim(file_name_ice)
1299         fileobj => fileobj_ice
1300         if (lowercase(trim(data_set)) == 'amip2' .or. &
1301             lowercase(trim(data_set)) == 'hurrell' .or. &
1302             lowercase(trim(data_set)) == 'daily') ncfieldname = 'ice' ! modified by JHC
1303      endif
1304
1305      dy = 0 ! only processing monthly data
1306
1307      if (verbose > 2 .and. mpp_pe() == 0)  &
1308           print *, 'looking for date = ', Date
1309
1310      ! This code can handle amip1, reynolds, or reyoi type SST data files in netCDF format
1311      if (mpp_pe() == mpp_root_pe()) call mpp_error ('amip_interp_mod', &
1312           'Reading NetCDF formatted input data file: '//trim(ncfilename), NOTE)
1313
1314         call fms2_io_read_data (fileobj, 'nrecords', nrecords)
1315         if (nrecords < 1) call mpp_error('amip_interp_mod', &
1316                            'Invalid number of SST records in SST datafile: '//trim(ncfilename), FATAL)
1317         allocate(ryr(nrecords), rmo(nrecords), rdy(nrecords))
1318         call fms2_io_read_data(fileobj, 'yr', ryr)
1319         call fms2_io_read_data(fileobj, 'mo', rmo)
1320         call fms2_io_read_data(fileobj, 'dy', rdy)
1321
1322      ierr = 1
1323      do k = 1, nrecords
1324        yr = ryr(k);  mo = rmo(k)
1325        Adate = date_type( yr, mo, 0)
1326        Curr_date = Adate
1327        if (verbose > 2 .and. mpp_pe() == 0)  &
1328              print *, '....... checking   ', Adate
1329        if (Date == Adate) ierr = 0
1330        if (yr == 0 .and. mo == Date%month) ierr = 0
1331        if (ierr == 0) exit
1332      enddo
1333      if (ierr .ne. 0) call mpp_error('amip_interp_mod', &
1334                       'Model time is out of range not in SST data: '//trim(ncfilename), FATAL)
1335         deallocate(ryr, rmo, rdy)
1336        !PRINT *, 'New SST data: ', k, yr, mo, dy, Date%year, Date%month, Date%day, ryr(1), rmo(1)
1337
1338    !---- check if climatological data should be used ----
1339
1340      if (yr == 0 .or. mo == 0) then
1341         ierr = 0
1342         if (date_out_of_range == 'fail' )               ierr = 1
1343         if (date_out_of_range == 'initclimo' .and.  &
1344              Date > Date_end )   ierr = 1
1345         if (ierr /= 0) call error_mesg &
1346              ('read_record in amip_interp_mod', &
1347              'climo data read when NO climo data requested', FATAL)
1348      endif
1349
1350    !---- read NETCDF data ----
1351
1352      if ( interp_oi_sst ) then
1353           call fms2_io_read_data(fileobj, ncfieldname, tmp_dat, unlim_dim_level=k)
1354 !     interpolate tmp_dat(360, 180) ---> dat(mobs,nobs) (to enable SST anom computation)
1355           if ( mobs/=360 .or. nobs/=180 ) then
1356                call a2a_bilinear(360, 180, tmp_dat, mobs, nobs, dat)
1357           else
1358                dat(:,:) = tmp_dat(:,:)
1359           endif
1360      else
1361           call fms2_io_read_data(fileobj, ncfieldname, dat, unlim_dim_level=k)
1362      endif
1363      idat =  nint(dat, I2_KIND) ! reconstruct packed data for reproducibility
1364
1365    !---- unpacking of data ----
1366
1367      if (type(1:3) == 'ice') then
1368         !---- create fractional [0,1] ice mask
1369         if (lowercase(trim(data_set)) /= 'amip2' .and. lowercase(trim(data_set)) /= 'hurrell') then
1370                where ( idat <= ice_crit )
1371                    dat = 1.
1372                elsewhere
1373                    dat = 0.
1374                endwhere
1375         else
1376            dat = dat*0.01
1377         endif
1378      else if (type(1:3) == 'sst') then
1379         !---- unpack sst ----
1380         if (lowercase(trim(data_set)) /= 'amip2' .and. lowercase(trim(data_set)) /= 'hurrell') then
1381                dat = real(idat)*0.01 + TFREEZE
1382         endif
1383      endif
1384
1385      return
1386
1387    end subroutine read_record
1388
1389 !#######################################################################
1390
1391    subroutine clip_data (type, dat)
1392
1393    character(len=*), intent(in)    :: type
1394    real,             intent(inout) :: dat(:,:)
1395
1396    if (type(1:3) == 'ice') then
1397        dat = min(max(dat,0.0),1.0)
1398    else if (type(1:3) == 'sst') then
1399        dat = max(tice_crit_k,dat)
1400    endif
1401
1402    end subroutine clip_data
1403
1404 !#######################################################################
1405
1406 !> @return logical answer
1407 function date_equals (Left, Right) result (answer)
1408 type (date_type), intent(in) :: Left, Right
1409 logical :: answer
1410
1411    if (Left % year  == Right % year  .and.  &
1412        Left % month == Right % month .and.  &
1413        Left % day   == Right % day ) then
1414            answer = .true.
1415    else
1416            answer = .false.
1417    endif
1418
1419 end function date_equals
1420
1421 !#######################################################################
1422
1423 !> @return logical answer
1424 function date_not_equals (Left, Right) result (answer)
1425 type (date_type), intent(in) :: Left, Right
1426 logical :: answer
1427
1428    if (Left % year  == Right % year  .and.  &
1429        Left % month == Right % month .and.  &
1430        Left % day   == Right % day ) then
1431            answer = .false.
1432    else
1433            answer = .true.
1434    endif
1435
1436 end function date_not_equals
1437
1438 !#######################################################################
1439
1440 !> @return logical answer
1441 function date_gt (Left, Right) result (answer)
1442 type (date_type), intent(in) :: Left, Right
1443 logical :: answer
1444 integer :: i, dif(3)
1445
1446    dif(1) = Left%year  - Right%year
1447    dif(2) = Left%month - Right%month
1448    dif(3) = Left%day   - Right%day
1449    answer = .false.
1450    do i = 1, 3
1451      if (dif(i) == 0) cycle
1452      if (dif(i)  < 0) exit
1453      if (dif(i)  > 0) then
1454          answer = .true.
1455          exit
1456      endif
1457    enddo
1458
1459 end function date_gt
1460
1461 !#######################################################################
1462
1463 subroutine print_dates (Time, Date1, Udate1,  &
1464                               Date2, Udate2, fmonth)
1465
1466    type (time_type), intent(in) :: Time
1467    type (date_type), intent(in) :: Date1, Udate1, Date2, Udate2
1468    real,             intent(in) :: fmonth
1469
1470    integer :: year, month, day, hour, minute, second
1471
1472    call get_date (Time, year, month, day, hour, minute, second)
1473
1474    write (*,10) year,month,day, hour,minute,second
1475    write (*,20) fmonth
1476    write (*,30) Date1, Udate1
1477    write (*,40) Date2, Udate2
1478
1479 10 format (/,' date(y/m/d h:m:s) = ',i4,2('/',i2.2),1x,2(i2.2,':'),i2.2)
1480 20 format (' fmonth = ',f9.7)
1481 30 format (' date1(y/m/d) = ',i4,2('/',i2.2),6x, &
1482                     'used = ',i4,2('/',i2.2),6x  )
1483 40 format (' date2(y/m/d) = ',i4,2('/',i2.2),6x, &
1484                     'used = ',i4,2('/',i2.2),6x  )
1485
1486 end subroutine print_dates
1487
1488 !#######################################################################
1489
1490 subroutine zonal_sst (Time, ice, sst)
1491
1492    type (time_type), intent(in)  :: Time
1493    real,             intent(out) :: ice(mobs,nobs), sst(mobs,nobs)
1494
1495    real    :: tpi, fdate, eps, ph, sph, sph2, ts
1496    integer :: j
1497
1498 ! namelist needed
1499 !
1500 !  teq  = sst at equator
1501 !  tdif = equator to pole sst difference
1502 !  tann = amplitude of annual cycle
1503 !  tlag = offset for time of year (for annual cycle)
1504 !
1505
1506     tpi = 2.0*pi
1507
1508     fdate = fraction_of_year (Time)
1509
1510     eps = sin( tpi*(fdate-tlag) ) * tann
1511
1512     do j = 1, nobs
1513
1514         ph  = 0.5*(lat_bnd(j)+lat_bnd(j+1))
1515        sph  = sin(ph)
1516        sph2 = sph*sph
1517
1518        ts = teq - tdif*sph2 - eps*sph
1519
1520        sst(:,j) = ts
1521
1522     enddo
1523
1524     where ( sst < tice_crit_k )
1525        ice = 1.0
1526        sst = tice_crit_k
1527     elsewhere
1528        ice  = 0.0
1529     endwhere
1530
1531
1532 end subroutine zonal_sst
1533
1534 !#######################################################################
1535
1536 subroutine amip_interp_type_eq(amip_interp_out, amip_interp_in)
1537     type(amip_interp_type), intent(inout) :: amip_interp_out
1538     type(amip_interp_type), intent(in)    :: amip_interp_in
1539
1540     if(.not.amip_interp_in%I_am_initialized) then
1541       call mpp_error(FATAL,'amip_interp_type_eq: amip_interp_type variable on right hand side is unassigned')
1542     endif
1543
1544     amip_interp_out%Hintrp     =  amip_interp_in%Hintrp
1545     amip_interp_out%data1      => amip_interp_in%data1
1546     amip_interp_out%data2      => amip_interp_in%data2
1547     amip_interp_out%Date1      =  amip_interp_in%Date1
1548     amip_interp_out%Date2      =  amip_interp_in%Date2
1549     amip_interp_out%Date1      =  amip_interp_in%Date1
1550     amip_interp_out%Date2      =  amip_interp_in%Date2
1551     amip_interp_out%use_climo  =  amip_interp_in%use_climo
1552     amip_interp_out%use_annual =  amip_interp_in%use_annual
1553     amip_interp_out%I_am_initialized = .true.
1554
1555 end subroutine amip_interp_type_eq
1556
1557 !#######################################################################
1558
1559 end module amip_interp_mod
1560 !> @}
1561 ! <INFO>
1562
1563 !   <FUTURE>
1564 !     Add AMIP 2 data set.
1565 !
1566 !     Other data sets (or extend current data sets).
1567 !   </FUTURE>
1568
1569 ! </INFO>