FMS  2025.04
Flexible Modeling System
amip_interp.inc
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 
19 ! modified by JHC
20 !> Retrieve sea surface temperature data and interpolated grid
21 subroutine get_amip_sst_ (Time, Interp, sst, err_msg, lon_model, lat_model)
22  type (time_type), intent(in) :: Time !< Time to interpolate
23  type (amip_interp_type), target, intent(inout) :: Interp !< Holds data for interpolation
24  real(FMS_AMIP_INTERP_KIND_), intent(out) :: sst(:,:) !< Sea surface temperature data
25  character(len=*), optional, intent(out) :: err_msg !< Holds error message string if present
26 
27  real(FMS_AMIP_INTERP_KIND_), dimension(mobs,nobs) :: sice
28  real(FMS_AMIP_INTERP_KIND_), allocatable, save :: temp1(:,:), temp2(:,:)
29 
30  integer :: year1, year2, month1, month2
31  real(FMS_AMIP_INTERP_KIND_) :: fmonth
32  type (date_type) :: Date1, Date2, Udate1, Udate2
33 
34  type(time_type) :: Amip_Time
35  integer :: tod(3),dum(3)
36 
37 ! add by JHC
38  real(FMS_AMIP_INTERP_KIND_), intent(in), dimension(:,:), optional :: lon_model, lat_model
39  real(FMS_AMIP_INTERP_KIND_) :: pert
40  integer :: i, j, mobs_sst, nobs_sst
41  integer :: jhctod(6)
42  type (time_type) :: Udate
43  character(len=4) :: yyyy
44  integer :: nrecords, ierr, k, yr, mo, dy
45  integer, dimension(:), allocatable :: ryr, rmo, rdy
46  character(len=30) :: time_unit
47  real(FMS_AMIP_INTERP_KIND_), dimension(:), allocatable :: timeval
48  character(len=FMS_FILE_LEN) :: ncfilename
49  type(FmsNetcdfFile_t) :: fileobj
50  logical :: the_file_exists
51 ! end add by JHC
52  logical, parameter :: DEBUG = .false. !> switch for debugging output
53  !> These are fms_io specific
54  integer :: iunit
55  integer, parameter :: lkind = fms_amip_interp_kind_
56 
57  if(present(err_msg)) err_msg = ''
58  if(.not.interp%I_am_initialized) then
59  if(fms_error_handler('get_amip_sst','The amip_interp_type variable is not initialized',err_msg)) return
60  endif
61 
62 !-----------------------------------------------------------------------
63 !----- compute zonally symetric sst ---------------
64 
65  if ( use_ncep_sst .and. forecast_mode ) no_anom_sst = .false.
66 
67  if (all(amip_date>0)) then
68  call get_date(time,dum(1),dum(2),dum(3),tod(1),tod(2),tod(3))
69  amip_time = set_date(amip_date(1),amip_date(2),amip_date(3),tod(1),tod(2),tod(3))
70  else
71  amip_time = time
72  endif
73 
74 ! add by JHC
75 if ( .not.use_daily ) then
76 ! end add by JHC
77 
78  if ( .not. allocated(temp1) ) allocate (temp1(mobs,nobs))
79  if ( .not. allocated(temp2) ) allocate (temp2(mobs,nobs))
80 
81  if (use_zonal) then
82  call zonal_sst_ (amip_time, sice, temp1)
83  call horiz_interp (interp%Hintrp, temp1, sst)
84  else
85 
86 !-----------------------------------------------------------------------
87 !---------- get new observed sea surface temperature -------------------
88 
89 ! ---- time interpolation for months -----
90  call time_interp (amip_time, fmonth, year1, year2, month1, month2)
91 ! ---- force climatology ----
92  if (interp%use_climo) then
93  year1=0; year2=0
94  endif
95  if (interp%use_annual) then
96  year1=0; year2=0
97  month1=0; month2=0
98  endif
99 ! ---------------------------
100 
101  date1 = date_type( year1, month1, 0 )
102  date2 = date_type( year2, month2, 0 )
103 
104 ! -- open/rewind file --
105  iunit = -1
106 !-----------------------------------------------------------------------
107 
108  if (date1 /= interp%Date1) then
109 ! ---- use Date2 for Date1 ----
110  if (date1 == interp%Date2) then
111  interp%Date1 = interp%Date2
112  interp%DATA1_ = interp%DATA2_
113  temp1(:,:) = temp2(:,:) ! SJL BUG fix: June 24, 2011
114  else
115  call read_record_ ('sst', date1, udate1, temp1)
116  call horiz_interp ( interp%Hintrp, temp1, interp%DATA1_)
117  call clip_data_ ('sst', interp%DATA1_)
118  interp%Date1 = date1
119  endif
120  endif
121 
122 !-----------------------------------------------------------------------
123 
124  if (date2 /= interp%Date2) then
125  call read_record_ ('sst', date2, udate2, temp2)
126  call horiz_interp ( interp%Hintrp, temp2, interp%DATA2_)
127  call clip_data_ ('sst', interp%DATA2_)
128  interp%Date2 = date2
129  endif
130 
131 !-----------------------------------------------------------------------
132 !---------- time interpolation (between months) of sst's ---------------
133 !-----------------------------------------------------------------------
134  sst = interp%DATA1_ + fmonth * (interp%DATA2_ - interp%DATA1_)
135 
136 !-------------------------------------------------------------------------------
137 ! SJL mods for NWP and TCSF ---
138 ! Nudging runs: (Note: NCEP SST updated only every 6-hr)
139 ! Compute SST anomaly from global SST datasets for subsequent forecast runs
140 !-------------------------------------------------------------------------------
141 
142 !! DEBUG CODE
143  if (debug) then
144  call get_date(amip_time,jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6))
145  if (mpp_pe() == 0) then
146  write (*,200) 'JHC: use_daily = F, AMIP_Time: ',jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5), &
147  & jhctod(6)
148  write (*,300) 'JHC: use_daily = F, interped SST: ', sst(1,1),sst(5,5),sst(10,10)
149  endif
150  endif
151 
152 
153  endif
154 
155 ! add by JHC
156 else
157  call get_date(amip_time,jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6))
158  if (mpp_pe() == mpp_root_pe()) write(*,200) 'amip_interp_mod: use_daily = T, Amip_Time = ',jhctod(1), &
159  & jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6)
160 
161  yr = jhctod(1); mo = jhctod(2); dy = jhctod(3)
162 
163  write (yyyy,'(i4)') jhctod(1)
164 
165  file_name_sst = 'INPUT/' // 'sst.day.mean.'//yyyy//'.v2.nc'
166  ncfilename = trim(file_name_sst)
167  time_unit = 'days since 1978-01-01 00:00:00'
168 
169  mobs_sst = 1440; nobs_sst = 720
170 
171  call set_sst_grid_edges_daily_ (mobs_sst, nobs_sst)
172  call horiz_interp_new ( interp%Hintrp2, lon_bnd_, lat_bnd_, &
173  lon_model, lat_model, interp_method="bilinear" )
174 
175  the_file_exists = fms2_io_file_exists(ncfilename)
176 
177  if ( (.NOT. the_file_exists) ) then
178  call mpp_error ('amip_interp_mod', &
179  'cannot find daily SST input data file: '//trim(ncfilename), note)
180  else
181  if (mpp_pe() == mpp_root_pe()) call mpp_error ('amip_interp_mod', &
182  'Reading NetCDF formatted daily SST from: '//trim(ncfilename), note)
183 
184  if(.not. open_file(fileobj, trim(ncfilename), 'read')) &
185  call error_mesg ('get_amip_sst', 'Error in opening file '//trim(ncfilename), fatal)
186 
187  call get_dimension_size(fileobj, 'TIME', nrecords)
188  if (nrecords < 1) call mpp_error('amip_interp_mod', &
189  'Invalid number of SST records in daily SST data file: '//trim(ncfilename), fatal)
190  allocate(timeval(nrecords), ryr(nrecords), rmo(nrecords), rdy(nrecords))
191  call fms2_io_read_data(fileobj, 'TIME', timeval)
192 !!! DEBUG CODE
193  if(debug) then
194  if (mpp_pe() == 0) then
195  print *, 'JHC: nrecords = ', nrecords
196  print *, 'JHC: TIME = ', timeval
197  endif
198  endif
199 
200  ierr = 1
201  do k = 1, nrecords
202 
203  udate = get_cal_time(timeval(k), time_unit, 'julian')
204  call get_date(udate,jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6))
205  ryr(k) = jhctod(1); rmo(k) = jhctod(2); rdy(k) = jhctod(3)
206 
207  if ( yr == ryr(k) .and. mo == rmo(k) .and. dy == rdy(k) ) ierr = 0
208  if (ierr==0) exit
209 
210  enddo
211 
212  if(debug) then
213  if (mpp_pe() == 0) then
214  print *, 'JHC: k =', k
215  print *, 'JHC: ryr(k) rmo(k) rdy(k)',ryr(k), rmo(k), rdy(k)
216  print *, 'JHC: yr mo dy ',yr, mo, dy
217  endif
218  endif
219 
220  if (ierr .ne. 0) call mpp_error('amip_interp_mod', &
221  'Model time is out of range not in SST data: '//trim(ncfilename), fatal)
222  endif ! if(file_exist(ncfilename))
223 
224 
225  !---- read NETCDF data ----
226  if ( .not. allocated(tempamip) ) &
227  & allocate (tempamip(mobs_sst,nobs_sst))
228 
229  if (the_file_exists) then
230  call fms2_io_read_data(fileobj, 'SST', tempamip, unlim_dim_level=k)
231  call close_file(fileobj)
232  tempamip = tempamip + tfreeze
233 
234 !!! DEBUG CODE
235  if(debug) then
236  if (mpp_pe() == 0) then
237  print*, 'JHC: TFREEZE = ', real(tfreeze, fms_amip_interp_kind_)
238  print*, lbound(sst)
239  print*, ubound(sst)
240  print*, lbound(tempamip)
241  print*, ubound(tempamip)
242  write(*,300) 'JHC: tempamip : ', tempamip(100,100), tempamip(200,200), tempamip(300,300)
243  endif
244  endif
245 
246  call horiz_interp ( interp%Hintrp2, tempamip_, sst )
247  call clip_data_ ('sst', sst)
248 
249  endif
250 
251  if(debug) then
252  if (mpp_pe() == 400) then
253  write(*,300)'JHC: use_daily = T, daily SST: ', sst(1,1),sst(5,5),sst(10,10)
254  print *,'JHC: use_daily = T, daily SST: ', sst
255  endif
256  endif
257 
258 200 format(a35, 6(i5,1x))
259 300 format(a35, 3(f7.3,2x))
260 
261 endif
262 ! end add by JHC
263 
264 ! add by JHC: add on non-zero sea surface temperature perturbation (namelist option)
265 ! This perturbation may be useful in accessing model sensitivities
266 
267  if ( do_sst_pert ) then
268 
269  if ( trim(sst_pert_type) == 'fixed' ) then
270  sst = sst + real(sst_pert, fms_amip_interp_kind_)
271  else if ( trim(sst_pert_type) == 'random' ) then
272  call random_seed()
273 
274  if(debug) then
275  if (mpp_pe() == 0) then
276  print*, 'mobs = ', mobs
277  print*, 'nobs = ', nobs
278  print*, lbound(sst)
279  print*, ubound(sst)
280  endif
281  endif
282 
283  do i = 1, size(sst,1)
284  do j = 1, size(sst,2)
285  call random_number(pert)
286  sst(i,j) = sst(i,j) + real(sst_pert, fms_amip_interp_kind_)*((pert-0.5_lkind)*2)
287  end do
288  end do
289  endif
290 
291  endif
292 ! end add by JHC
293  end subroutine get_amip_sst_
294 
295 !> AMIP interpolation for ice
296 subroutine get_amip_ice_ (Time, Interp, ice, err_msg)
297  type (time_type), intent(in) :: Time !< Time to interpolate
298  type (amip_interp_type), target, intent(inout) :: Interp !< Holds data for interpolation
299  real(FMS_AMIP_INTERP_KIND_), intent(out) :: ice(:,:) !< ice data
300  character(len=*), optional, intent(out) :: err_msg !< Holds error message string if present
301 
302  real(FMS_AMIP_INTERP_KIND_), dimension(mobs,nobs) :: sice, temp
303 
304  integer :: year1, year2, month1, month2
305  real(FMS_AMIP_INTERP_KIND_) :: fmonth
306  type (date_type) :: Date1, Date2, Udate1, Udate2
307 
308  type(time_type) :: Amip_Time
309  integer :: tod(3),dum(3)
310  integer, parameter :: lkind = fms_amip_interp_kind_
311 
312  if(present(err_msg)) err_msg = ''
313  if(.not.interp%I_am_initialized) then
314  if(fms_error_handler('get_amip_ice','The amip_interp_type variable is not initialized',err_msg)) return
315  endif
316 
317 !-----------------------------------------------------------------------
318 !----- compute zonally symetric sst ---------------
319 
320 
321  if (any(amip_date>0)) then
322 
323  call get_date(time,dum(1),dum(2),dum(3),tod(1),tod(2),tod(3))
324 
325  amip_time = set_date(amip_date(1),amip_date(2),amip_date(3),tod(1),tod(2),tod(3))
326 
327  else
328 
329  amip_time = time
330 
331  endif
332 
333 
334 if (use_zonal) then
335  call zonal_sst_ (amip_time, sice, temp)
336  call horiz_interp ( interp%Hintrp, sice, ice )
337 else
338 
339 !-----------------------------------------------------------------------
340 !---------- get new observed sea surface temperature -------------------
341 
342 ! ---- time interpolation for months -----
343 
344  call time_interp (amip_time, fmonth, year1, year2, month1, month2)
345 
346 ! ---- force climatology ----
347  if (interp%use_climo) then
348  year1=0; year2=0
349  endif
350  if (interp%use_annual) then
351  year1=0; year2=0
352  month1=0; month2=0
353  endif
354 ! ---------------------------
355 
356  date1 = date_type( year1, month1, 0 )
357  date2 = date_type( year2, month2, 0 )
358 
359  iunit = -1
360 !-----------------------------------------------------------------------
361 
362  if (date1 /= interp%Date1) then
363 ! ---- use Date2 for Date1 ----
364  if (date1 == interp%Date2) then
365  interp%Date1 = interp%Date2
366  interp%DATA1_ = interp%DATA2_
367  else
368 !-- SJL -------------------------------------------------------------
369 ! Can NOT use ncep_sst to determine sea_ice For seasonal forecast
370 ! Use climo sea ice for seasonal runs
371  call read_record_ ('ice', date1, udate1, sice)
372 !--------------------------------------------------------------------
373  call horiz_interp ( interp%Hintrp, sice, interp%DATA1_)
374  call clip_data_ ('ice', interp%DATA1_)
375  interp%Date1 = date1
376  endif
377  endif
378 
379 !-----------------------------------------------------------------------
380 
381  if (date2 /= interp%Date2) then
382 
383 !-- SJL -------------------------------------------------------------
384  call read_record_ ('ice', date2, udate2, sice)
385 !--------------------------------------------------------------------
386  call horiz_interp ( interp%Hintrp, sice, interp%DATA2_)
387  call clip_data_ ('ice', interp%DATA2_)
388  interp%Date2 = date2
389 
390  endif
391 
392 !-----------------------------------------------------------------------
393 !---------- time interpolation (between months) ------------------------
394 !-----------------------------------------------------------------------
395 
396  ice = interp%DATA1_ + fmonth * (interp%DATA2_ - interp%DATA1_)
397 
398 endif
399  end subroutine get_amip_ice_
400 
401  !> @return A newly created @ref amip_interp_type
402  function amip_interp_new_1d_ ( lon , lat , mask , use_climo, use_annual, &
403  interp_method ) result (Interp)
404  real(FMS_AMIP_INTERP_KIND_), intent(in), dimension(:) :: lon, lat
405  logical, intent(in), dimension(:,:) :: mask
406  character(len=*), intent(in), optional :: interp_method
407  logical, intent(in), optional :: use_climo, use_annual
408 
409  type (amip_interp_type) :: Interp
410 
411  if(.not.module_is_initialized) call amip_interp_init
412 
413  interp%use_climo = .false.
414  if (present(use_climo)) interp%use_climo = use_climo
415  interp%use_annual = .false.
416  if (present(use_annual)) interp%use_annual = use_annual
417 
418  if ( date_out_of_range == 'fail' .and. interp%use_climo ) &
419  call error_mesg ('amip_interp_new_1d', 'use_climo mismatch', fatal)
420 
421  if ( date_out_of_range == 'fail' .and. interp%use_annual ) &
422  call error_mesg ('amip_interp_new_1d', 'use_annual(climo) mismatch', fatal)
423 
424  interp%Date1 = date_type( -99, -99, -99 )
425  interp%Date2 = date_type( -99, -99, -99 )
426 
427 !-----------------------------------------------------------------------
428 ! ---- initialization of horizontal interpolation ----
429 
430  call horiz_interp_new ( interp%Hintrp, lon_bnd_, lat_bnd_, &
431  lon, lat, interp_method= interp_method )
432 
433  allocate(interp%DATA1_ (size(lon(:))-1,size(lat(:))-1))
434  allocate(interp%DATA2_ (size(lon(:))-1,size(lat(:))-1))
435 
436  interp%I_am_initialized = .true.
437  end function amip_interp_new_1d_
438 
439  !> @return A newly created @ref amip_interp_type
440  function amip_interp_new_2d_ ( lon , lat , mask , use_climo, use_annual, &
441  interp_method ) result (Interp)
442  real(FMS_AMIP_INTERP_KIND_), intent(in), dimension(:,:) :: lon, lat
443  logical, intent(in), dimension(:,:) :: mask
444  character(len=*), intent(in), optional :: interp_method
445  logical, intent(in), optional :: use_climo, use_annual
446 
447  type (amip_interp_type) :: Interp
448 
449  if(.not.module_is_initialized) call amip_interp_init
450 
451  interp%use_climo = .false.
452  if (present(use_climo)) interp%use_climo = use_climo
453  interp%use_annual = .false.
454  if (present(use_annual)) interp%use_annual = use_annual
455 
456  if ( date_out_of_range == 'fail' .and. interp%use_climo ) &
457  call error_mesg ('amip_interp_new_2d', 'use_climo mismatch', fatal)
458 
459  if ( date_out_of_range == 'fail' .and. interp%use_annual ) &
460  call error_mesg ('amip_interp_new_2d', 'use_annual(climo) mismatch', fatal)
461 
462  interp%Date1 = date_type( -99, -99, -99 )
463  interp%Date2 = date_type( -99, -99, -99 )
464 
465 !-----------------------------------------------------------------------
466 ! ---- initialization of horizontal interpolation ----
467 
468  call horiz_interp_new ( interp%Hintrp, lon_bnd_, lat_bnd_, &
469  lon, lat, interp_method = interp_method)
470 
471  allocate(interp%DATA1_ (size(lon,1),size(lat,2)))
472  allocate(interp%DATA2_ (size(lon,1),size(lat,2)))
473 
474  interp%I_am_initialized = .true.
475  end function amip_interp_new_2d_
476 
477 ! add by JHC
478  subroutine set_sst_grid_edges_daily_ (mobs_sst, nobs_sst)
479  integer :: i, j, mobs_sst, nobs_sst
480  real(FMS_AMIP_INTERP_KIND_) :: hpie, dlon, dlat, wb, sb
481  integer, parameter :: lkind = fms_amip_interp_kind_
482 
483  if(allocated(lon_bnd)) deallocate(lon_bnd)
484  if(allocated(lat_bnd)) deallocate(lat_bnd)
485 
486  allocate(lon_bnd(mobs_sst+1))
487  allocate(lat_bnd(nobs_sst+1))
488 
489 ! ---- compute grid edges (do only once) -----
490 
491  hpie = pi / 2._r8_kind
492  dlon = 4._r8_kind*hpie/real(mobs_sst, r8_kind)
493  wb = 0.0_r8_kind
494 
495  lon_bnd(1) = wb
496  do i = 2, mobs_sst+1
497  lon_bnd(i) = wb + dlon * real(i-1, r8_kind)
498  enddo
499  lon_bnd(mobs_sst+1) = lon_bnd(1) + 4._r8_kind*hpie
500 
501  dlat = 2._r8_kind*hpie/real(nobs_sst, r8_kind)
502  sb = -hpie
503 
504  lat_bnd(1) = sb
505  lat_bnd(nobs_sst+1) = hpie
506  do j = 2, nobs_sst
507  lat_bnd(j) = sb + dlat * real(j-1, r8_kind)
508  enddo
509  end subroutine set_sst_grid_edges_daily_
510 ! end add by JHC
511 
512  subroutine a2a_bilinear_ (nx, ny, dat1, n1, n2, dat2)
513  integer, intent(in) :: nx, ny
514  integer, intent(in) :: n1, n2
515  real(FMS_AMIP_INTERP_KIND_), intent(in) :: dat1(nx,ny)
516  real(FMS_AMIP_INTERP_KIND_), intent(out) :: dat2(n1,n2) !> output interpolated data
517 
518 ! local:
519  real(FMS_AMIP_INTERP_KIND_) :: lon1(nx), lat1(ny)
520  real(FMS_AMIP_INTERP_KIND_) :: lon2(n1), lat2(n2)
521  real(FMS_AMIP_INTERP_KIND_) :: dx1, dy1, dx2, dy2
522  real(FMS_AMIP_INTERP_KIND_) :: xc, yc
523  real(FMS_AMIP_INTERP_KIND_) :: a1, b1, c1, c2, c3, c4
524  integer :: i1, i2, jc, i0, j0, it, jt
525  integer :: i, j
526  integer, parameter :: lkind = fms_amip_interp_kind_
527 
528 
529 !-----------------------------------------------------------
530 ! * Interpolate from "FMS" 1x1 SST data grid to a finer grid
531 ! lon: 0.5, 1.5, ..., 359.5
532 ! lat: -89.5, -88.5, ... , 88.5, 89.5
533 !-----------------------------------------------------------
534 
535  dx1 = 360._lkind/real(nx, fms_amip_interp_kind_) !> INput Grid
536  dy1 = 180._lkind/real(ny, fms_amip_interp_kind_) !> INput Grid
537 
538  do i=1,nx
539  lon1(i) = 0.5_lkind*dx1 + real(i-1, fms_amip_interp_kind_)*dx1
540  enddo
541  do j=1,ny
542  lat1(j) = -90._lkind + 0.5_lkind*dy1 + real(j-1, fms_amip_interp_kind_)*dy1
543  enddo
544 
545  dx2 = 360._lkind/real(n1, fms_amip_interp_kind_) !> OutPut Grid:
546  dy2 = 180._lkind/real(n2, fms_amip_interp_kind_) !> OutPut Grid:
547 
548  do i=1,n1
549  lon2(i) = 0.5_lkind*dx2 + real(i-1, fms_amip_interp_kind_)*dx2
550  enddo
551  do j=1,n2
552  lat2(j) = -90._lkind + 0.5_lkind*dy2 + real(j-1, fms_amip_interp_kind_)*dy2
553  enddo
554 
555  jt = 1
556  do 5000 j=1,n2
557 
558  yc = lat2(j)
559  if ( yc<lat1(1) ) then
560  jc = 1
561  b1 = 0._lkind
562  elseif ( yc>lat1(ny) ) then
563  jc = ny-1
564  b1 = 1._lkind
565  else
566  do j0=jt,ny-1
567  if ( yc>=lat1(j0) .and. yc<=lat1(j0+1) ) then
568  jc = j0
569  jt = j0
570  b1 = (yc-lat1(jc)) / dy1
571  go to 222
572  endif
573  enddo
574  endif
575 222 continue
576 
577  it = 1
578  do i=1,n1
579  xc = lon2(i)
580  if ( xc>lon1(nx) ) then
581  i1 = nx; i2 = 1
582  a1 = (xc-lon1(nx)) / dx1
583  elseif ( xc<lon1(1) ) then
584  i1 = nx; i2 = 1
585  a1 = (xc+360._lkind-lon1(nx)) / dx1
586  else
587  do i0=it,nx-1
588  if ( xc>=lon1(i0) .and. xc<=lon1(i0+1) ) then
589  i1 = i0; i2 = i0+1
590  it = i0
591  a1 = (xc-lon1(i1)) / dx1
592  go to 111
593  endif
594  enddo
595  endif
596 111 continue
597 
598 ! Debug code:
599  if ( a1<-0.001_lkind .or. a1>1.001_lkind .or. b1<-0.001_lkind .or. b1>1.001_lkind ) then
600  write(*,*) i,j,a1, b1
601  call mpp_error(fatal,'a2a bilinear interpolation')
602  endif
603 
604  c1 = (1._lkind-a1) * (1._lkind-b1)
605  c2 = a1 * (1._lkind-b1)
606  c3 = a1 * b1
607  c4 = (1._lkind-a1) * b1
608 
609 ! Bilinear interpolation:
610  dat2(i,j) = c1*dat1(i1,jc) + c2*dat1(i2,jc) + c3*dat1(i2,jc+1) + c4*dat1(i1,jc+1)
611 
612  enddo !i-loop
613 
614 5000 continue ! j-loop
615  end subroutine a2a_bilinear_
616 
617  subroutine read_record_ (type, Date, Adate, dat)
618  character(len=*), intent(in) :: type
619  type (date_type), intent(in) :: Date
620  type (date_type), intent(inout) :: Adate
621  real(FMS_AMIP_INTERP_KIND_), intent(out) :: dat(mobs,nobs)
622  real(FMS_AMIP_INTERP_KIND_) :: tmp_dat(360,180)
623 
624  integer(I2_KIND) :: idat(mobs,nobs)
625  integer :: nrecords, yr, mo, dy, ierr, k
626  integer, dimension(:), allocatable :: ryr, rmo, rdy
627  character(len=FMS_FILE_LEN) :: ncfilename
628  character(len=NF90_MAX_NAME) :: ncfieldname
629  type(FmsNetcdfFile_t), pointer :: fileobj
630  integer, parameter :: lkind = fms_amip_interp_kind_
631 
632  !---- set file and field name for NETCDF data sets ----
633 
634  ncfieldname = 'sst'
635  if(type(1:3) == 'sst') then
636  ncfilename = trim(file_name_sst)
637  fileobj => fileobj_sst
638  else if(type(1:3) == 'ice') then
639  ncfilename = trim(file_name_ice)
640  fileobj => fileobj_ice
641  if (lowercase(trim(data_set)) == 'amip2' .or. &
642  lowercase(trim(data_set)) == 'hurrell' .or. &
643  lowercase(trim(data_set)) == 'daily') ncfieldname = 'ice' ! modified by JHC
644  endif
645 
646  dy = 0 ! only processing monthly data
647 
648  if (verbose > 2 .and. mpp_pe() == 0) &
649  print *, 'looking for date = ', date
650 
651  ! This code can handle amip1, reynolds, or reyoi type SST data files in netCDF format
652  if (mpp_pe() == mpp_root_pe()) call mpp_error ('amip_interp_mod', &
653  'Reading NetCDF formatted input data file: '//trim(ncfilename), note)
654 
655  call fms2_io_read_data (fileobj, 'nrecords', nrecords)
656  if (nrecords < 1) call mpp_error('amip_interp_mod', &
657  'Invalid number of SST records in SST datafile: '//trim(ncfilename), fatal)
658  allocate(ryr(nrecords), rmo(nrecords), rdy(nrecords))
659  call fms2_io_read_data(fileobj, 'yr', ryr)
660  call fms2_io_read_data(fileobj, 'mo', rmo)
661  call fms2_io_read_data(fileobj, 'dy', rdy)
662 
663  ierr = 1
664  do k = 1, nrecords
665  yr = ryr(k); mo = rmo(k)
666  adate = date_type( yr, mo, 0)
667  curr_date = adate
668  if (verbose > 2 .and. mpp_pe() == 0) &
669  print *, '....... checking ', adate
670  if (date == adate) ierr = 0
671  if (yr == 0 .and. mo == date%month) ierr = 0
672  if (ierr == 0) exit
673  enddo
674  if (ierr .ne. 0) call mpp_error('amip_interp_mod', &
675  'Model time is out of range not in SST data: '//trim(ncfilename), fatal)
676  deallocate(ryr, rmo, rdy)
677  !PRINT *, 'New SST data: ', k, yr, mo, dy, Date%year, Date%month, Date%day, ryr(1), rmo(1)
678 
679  !---- check if climatological data should be used ----
680 
681  if (yr == 0 .or. mo == 0) then
682  ierr = 0
683  if (date_out_of_range == 'fail' ) ierr = 1
684  if (date_out_of_range == 'initclimo' .and. &
685  date > date_end ) ierr = 1
686  if (ierr /= 0) call error_mesg &
687  ('read_record in amip_interp_mod', &
688  'climo data read when NO climo data requested', fatal)
689  endif
690 
691  !---- read NETCDF data ----
692 
693  if ( interp_oi_sst ) then
694  call fms2_io_read_data(fileobj, ncfieldname, tmp_dat, unlim_dim_level=k)
695 ! interpolate tmp_dat(360, 180) ---> dat(mobs,nobs) (to enable SST anom computation)
696  if ( mobs/=360 .or. nobs/=180 ) then
697  call a2a_bilinear_ (360, 180, tmp_dat, mobs, nobs, dat)
698  else
699  dat(:,:) = tmp_dat(:,:)
700  endif
701  else
702  call fms2_io_read_data(fileobj, ncfieldname, dat, unlim_dim_level=k)
703  endif
704  !TODO This assumes that the data is "packed" (has the scale_factor and add_offset attributes)
705  ! in fms2_io_read_data the data is unpacked (data_in_file*scale_factor + add_offset)
706  ! the line below "packs" the data again. This is needed for reproducibility
707  idat = nint(dat*100._lkind, i2_kind)
708 
709  !---- unpacking of data ----
710 
711  if (type(1:3) == 'ice') then
712  !---- create fractional [0,1] ice mask
713  if (lowercase(trim(data_set)) /= 'amip2' .and. lowercase(trim(data_set)) /= 'hurrell') then
714  where ( idat <= ice_crit )
715  dat = 1._lkind
716  elsewhere
717  dat = 0._lkind
718  endwhere
719  else
720  dat = dat*0.01_lkind
721  endif
722  else if (type(1:3) == 'sst') then
723  !---- unpack sst ----
724  if (lowercase(trim(data_set)) /= 'amip2' .and. lowercase(trim(data_set)) /= 'hurrell') then
725  dat = real(idat, fms_amip_interp_kind_)*0.01_lkind + real(tfreeze, fms_amip_interp_kind_)
726  endif
727  endif
728 
729  return
730  end subroutine read_record_
731 
732  subroutine clip_data_ (type, dat)
733  character(len=*), intent(in) :: type
734  real(FMS_AMIP_INTERP_KIND_), intent(inout) :: dat(:,:)
735  integer, parameter :: lkind = fms_amip_interp_kind_
736 
737  if (type(1:3) == 'ice') then
738  dat = min(max(dat,0.0_lkind), 1.0_lkind)
739  else if (type(1:3) == 'sst') then
740  dat = max(real(tice_crit_k, fms_amip_interp_kind_),dat)
741  endif
742  end subroutine clip_data_
743 
744 subroutine zonal_sst_ (Time, ice, sst)
745  type (time_type), intent(in) :: Time
746  real(FMS_AMIP_INTERP_KIND_), intent(out) :: ice(mobs,nobs), sst(mobs,nobs)
747  real(FMS_AMIP_INTERP_KIND_) :: tpi, fdate, eps, ph, sph, sph2, ts
748  integer :: j
749  integer, parameter :: lkind = fms_amip_interp_kind_
750 
751 ! namelist needed
752 !
753 ! teq = sst at equator
754 ! tdif = equator to pole sst difference
755 ! tann = amplitude of annual cycle
756 ! tlag = offset for time of year (for annual cycle)
757 !
758 
759  tpi = 2.0_lkind*real(pi, fms_amip_interp_kind_)
760 
761  fdate = fraction_of_year(time)
762 
763  eps = sin( tpi*(fdate-real(tlag, fms_amip_interp_kind_)) ) * real(tann, fms_amip_interp_kind_)
764 
765  do j = 1, nobs
766 
767  ph = 0.5_lkind * real(lat_bnd(j)+lat_bnd(j+1), fms_amip_interp_kind_)
768  sph = sin(ph)
769  sph2 = sph*sph
770 
771  ts = real(teq, fms_amip_interp_kind_) - real(tdif, fms_amip_interp_kind_)*sph2 - eps*sph
772 
773  sst(:,j) = ts
774 
775  enddo
776 
777  where ( sst < real(tice_crit_k, fms_amip_interp_kind_) )
778  ice = 1.0_lkind
779  sst = real(tice_crit_k, fms_amip_interp_kind_)
780  elsewhere
781  ice = 0.0_lkind
782  endwhere
783 end subroutine zonal_sst_
mpp_pe
integer function mpp_pe()
Returns processor ID.
Definition: mpp_util.inc:406