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