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