22 subroutine get_amip_sst_ (Time, Interp, sst, err_msg, lon_model, lat_model)
23 type (time_type),
intent(in) :: Time
24 type (amip_interp_type),
target,
intent(inout) :: Interp
25 real(FMS_AMIP_INTERP_KIND_),
intent(out) :: sst(:,:)
26 character(len=*),
optional,
intent(out) :: err_msg
28 real(FMS_AMIP_INTERP_KIND_),
dimension(mobs,nobs) :: sice
29 real(FMS_AMIP_INTERP_KIND_),
allocatable,
save :: temp1(:,:), temp2(:,:)
31 integer :: year1, year2, month1, month2
32 real(FMS_AMIP_INTERP_KIND_) :: fmonth
33 type (date_type) :: Date1, Date2, Udate1, Udate2
35 type(time_type) :: Amip_Time
36 integer :: tod(3),dum(3)
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
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
53 logical,
parameter :: DEBUG = .false.
56 integer,
parameter :: lkind = fms_amip_interp_kind_
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
66 if ( use_ncep_sst .and. forecast_mode ) no_anom_sst = .false.
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))
76 if ( .not.use_daily )
then
79 if ( .not.
allocated(temp1) )
allocate (temp1(mobs,nobs))
80 if ( .not.
allocated(temp2) )
allocate (temp2(mobs,nobs))
83 call zonal_sst_ (amip_time, sice, temp1)
84 call horiz_interp (interp%Hintrp, temp1, sst)
91 call time_interp (amip_time, fmonth, year1, year2, month1, month2)
93 if (interp%use_climo)
then
96 if (interp%use_annual)
then
102 date1 = date_type( year1, month1, 0 )
103 date2 = date_type( year2, month2, 0 )
109 if (date1 /= interp%Date1)
then
111 if (date1 == interp%Date2)
then
112 interp%Date1 = interp%Date2
113 interp%DATA1_ = interp%DATA2_
114 temp1(:,:) = temp2(:,:)
116 call read_record_ (
'sst', date1, udate1, temp1)
117 call horiz_interp ( interp%Hintrp, temp1, interp%DATA1_)
118 call clip_data_ (
'sst', interp%DATA1_)
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_)
135 sst = interp%DATA1_ + fmonth * (interp%DATA2_ - interp%DATA1_)
145 call get_date(amip_time,jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5),jhctod(6))
147 write (*,200)
'JHC: use_daily = F, AMIP_Time: ',jhctod(1),jhctod(2),jhctod(3),jhctod(4),jhctod(5), &
149 write (*,300)
'JHC: use_daily = F, interped SST: ', sst(1,1),sst(5,5),sst(10,10)
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)
162 yr = jhctod(1); mo = jhctod(2); dy = jhctod(3)
164 write (yyyy,
'(i4)') jhctod(1)
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'
170 mobs_sst = 1440; nobs_sst = 720
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" )
176 the_file_exists = fms2_io_file_exists(ncfilename)
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)
182 if (
mpp_pe() == mpp_root_pe())
call mpp_error (
'amip_interp_mod', &
183 'Reading NetCDF formatted daily SST from: '//trim(ncfilename), note)
185 if(.not. open_file(fileobj, trim(ncfilename),
'read')) &
186 call error_mesg (
'get_amip_sst',
'Error in opening file '//trim(ncfilename), fatal)
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)
196 print *,
'JHC: nrecords = ', nrecords
197 print *,
'JHC: TIME = ', timeval
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)
208 if ( yr == ryr(k) .and. mo == rmo(k) .and. dy == rdy(k) ) ierr = 0
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
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)
227 if ( .not.
allocated(tempamip) ) &
228 &
allocate (tempamip(mobs_sst,nobs_sst))
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
238 print*,
'JHC: TFREEZE = ', real(tfreeze, fms_amip_interp_kind_)
241 print*, lbound(tempamip)
242 print*, ubound(tempamip)
243 write(*,300)
'JHC: tempamip : ', tempamip(100,100), tempamip(200,200), tempamip(300,300)
247 call horiz_interp ( interp%Hintrp2, tempamip_, sst )
248 call clip_data_ (
'sst', sst)
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
259 200
format(a35, 6(i5,1x))
260 300
format(a35, 3(f7.3,2x))
268 if ( do_sst_pert )
then
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
277 print*,
'mobs = ', mobs
278 print*,
'nobs = ', nobs
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)
294 end subroutine get_amip_sst_
297 subroutine get_amip_ice_ (Time, Interp, ice, err_msg)
298 type (time_type),
intent(in) :: Time
299 type (amip_interp_type),
target,
intent(inout) :: Interp
300 real(FMS_AMIP_INTERP_KIND_),
intent(out) :: ice(:,:)
301 character(len=*),
optional,
intent(out) :: err_msg
303 real(FMS_AMIP_INTERP_KIND_),
dimension(mobs,nobs) :: sice, temp
305 integer :: year1, year2, month1, month2
306 real(FMS_AMIP_INTERP_KIND_) :: fmonth
307 type (date_type) :: Date1, Date2, Udate1, Udate2
309 type(time_type) :: Amip_Time
310 integer :: tod(3),dum(3)
311 integer,
parameter :: lkind = fms_amip_interp_kind_
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
322 if (any(amip_date>0))
then
324 call get_date(time,dum(1),dum(2),dum(3),tod(1),tod(2),tod(3))
326 amip_time = set_date(amip_date(1),amip_date(2),amip_date(3),tod(1),tod(2),tod(3))
336 call zonal_sst_ (amip_time, sice, temp)
337 call horiz_interp ( interp%Hintrp, sice, ice )
345 call time_interp (amip_time, fmonth, year1, year2, month1, month2)
348 if (interp%use_climo)
then
351 if (interp%use_annual)
then
357 date1 = date_type( year1, month1, 0 )
358 date2 = date_type( year2, month2, 0 )
363 if (date1 /= interp%Date1)
then
365 if (date1 == interp%Date2)
then
366 interp%Date1 = interp%Date2
367 interp%DATA1_ = interp%DATA2_
372 call read_record_ (
'ice', date1, udate1, sice)
374 call horiz_interp ( interp%Hintrp, sice, interp%DATA1_)
375 call clip_data_ (
'ice', interp%DATA1_)
382 if (date2 /= interp%Date2)
then
385 call read_record_ (
'ice', date2, udate2, sice)
387 call horiz_interp ( interp%Hintrp, sice, interp%DATA2_)
388 call clip_data_ (
'ice', interp%DATA2_)
397 ice = interp%DATA1_ + fmonth * (interp%DATA2_ - interp%DATA1_)
400 end subroutine get_amip_ice_
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
410 type (amip_interp_type) :: Interp
412 if(.not.module_is_initialized)
call amip_interp_init
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
419 if ( date_out_of_range ==
'fail' .and. interp%use_climo ) &
420 call error_mesg (
'amip_interp_new_1d',
'use_climo mismatch', fatal)
422 if ( date_out_of_range ==
'fail' .and. interp%use_annual ) &
423 call error_mesg (
'amip_interp_new_1d',
'use_annual(climo) mismatch', fatal)
425 interp%Date1 = date_type( -99, -99, -99 )
426 interp%Date2 = date_type( -99, -99, -99 )
431 call horiz_interp_new ( interp%Hintrp, lon_bnd_, lat_bnd_, &
432 lon, lat, interp_method= interp_method )
434 allocate(interp%DATA1_ (
size(lon(:))-1,
size(lat(:))-1))
435 allocate(interp%DATA2_ (
size(lon(:))-1,
size(lat(:))-1))
437 interp%I_am_initialized = .true.
438 end function amip_interp_new_1d_
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
448 type (amip_interp_type) :: Interp
450 if(.not.module_is_initialized)
call amip_interp_init
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
457 if ( date_out_of_range ==
'fail' .and. interp%use_climo ) &
458 call error_mesg (
'amip_interp_new_2d',
'use_climo mismatch', fatal)
460 if ( date_out_of_range ==
'fail' .and. interp%use_annual ) &
461 call error_mesg (
'amip_interp_new_2d',
'use_annual(climo) mismatch', fatal)
463 interp%Date1 = date_type( -99, -99, -99 )
464 interp%Date2 = date_type( -99, -99, -99 )
469 call horiz_interp_new ( interp%Hintrp, lon_bnd_, lat_bnd_, &
470 lon, lat, interp_method = interp_method)
472 allocate(interp%DATA1_ (
size(lon,1),
size(lat,2)))
473 allocate(interp%DATA2_ (
size(lon,1),
size(lat,2)))
475 interp%I_am_initialized = .true.
476 end function amip_interp_new_2d_
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_
484 if(
allocated(lon_bnd))
deallocate(lon_bnd)
485 if(
allocated(lat_bnd))
deallocate(lat_bnd)
487 allocate(lon_bnd(mobs_sst+1))
488 allocate(lat_bnd(nobs_sst+1))
492 hpie = pi / 2._r8_kind
493 dlon = 4._r8_kind*hpie/real(mobs_sst, r8_kind)
498 lon_bnd(i) = wb + dlon * real(i-1, r8_kind)
500 lon_bnd(mobs_sst+1) = lon_bnd(1) + 4._r8_kind*hpie
502 dlat = 2._r8_kind*hpie/real(nobs_sst, r8_kind)
506 lat_bnd(nobs_sst+1) = hpie
508 lat_bnd(j) = sb + dlat * real(j-1, r8_kind)
510 end subroutine set_sst_grid_edges_daily_
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)
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
527 integer,
parameter :: lkind = fms_amip_interp_kind_
536 dx1 = 360._lkind/real(nx, fms_amip_interp_kind_)
537 dy1 = 180._lkind/real(ny, fms_amip_interp_kind_)
540 lon1(i) = 0.5_lkind*dx1 + real(i-1, fms_amip_interp_kind_)*dx1
543 lat1(j) = -90._lkind + 0.5_lkind*dy1 + real(j-1, fms_amip_interp_kind_)*dy1
546 dx2 = 360._lkind/real(n1, fms_amip_interp_kind_)
547 dy2 = 180._lkind/real(n2, fms_amip_interp_kind_)
550 lon2(i) = 0.5_lkind*dx2 + real(i-1, fms_amip_interp_kind_)*dx2
553 lat2(j) = -90._lkind + 0.5_lkind*dy2 + real(j-1, fms_amip_interp_kind_)*dy2
560 if ( yc<lat1(1) )
then
563 elseif ( yc>lat1(ny) )
then
568 if ( yc>=lat1(j0) .and. yc<=lat1(j0+1) )
then
571 b1 = (yc-lat1(jc)) / dy1
581 if ( xc>lon1(nx) )
then
583 a1 = (xc-lon1(nx)) / dx1
584 elseif ( xc<lon1(1) )
then
586 a1 = (xc+360._lkind-lon1(nx)) / dx1
589 if ( xc>=lon1(i0) .and. xc<=lon1(i0+1) )
then
592 a1 = (xc-lon1(i1)) / dx1
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')
605 c1 = (1._lkind-a1) * (1._lkind-b1)
606 c2 = a1 * (1._lkind-b1)
608 c4 = (1._lkind-a1) * b1
611 dat2(i,j) = c1*dat1(i1,jc) + c2*dat1(i2,jc) + c3*dat1(i2,jc+1) + c4*dat1(i1,jc+1)
616 end subroutine a2a_bilinear_
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)
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_
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'
649 if (verbose > 2 .and.
mpp_pe() == 0) &
650 print *,
'looking for date = ', date
653 if (
mpp_pe() == mpp_root_pe())
call mpp_error (
'amip_interp_mod', &
654 'Reading NetCDF formatted input data file: '//trim(ncfilename), note)
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)
666 yr = ryr(k); mo = rmo(k)
667 adate = date_type( yr, mo, 0)
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
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)
682 if (yr == 0 .or. mo == 0)
then
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)
694 if ( interp_oi_sst )
then
695 call fms2_io_read_data(fileobj, ncfieldname, tmp_dat, unlim_dim_level=k)
697 if ( mobs/=360 .or. nobs/=180 )
then
698 call a2a_bilinear_ (360, 180, tmp_dat, mobs, nobs, dat)
700 dat(:,:) = tmp_dat(:,:)
703 call fms2_io_read_data(fileobj, ncfieldname, dat, unlim_dim_level=k)
708 idat = nint(dat*100._lkind, i2_kind)
712 if (
type(1:3) ==
'ice') then
714 if (lowercase(trim(data_set)) /=
'amip2' .and. lowercase(trim(data_set)) /=
'hurrell')
then
715 where ( idat <= ice_crit )
723 else if (
type(1:3) ==
'sst') then
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_)
731 end subroutine read_record_
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_
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)
743 end subroutine clip_data_
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
750 integer,
parameter :: lkind = fms_amip_interp_kind_
760 tpi = 2.0_lkind*real(pi, fms_amip_interp_kind_)
762 fdate = fraction_of_year(time)
764 eps = sin( tpi*(fdate-real(tlag, fms_amip_interp_kind_)) ) * real(tann, fms_amip_interp_kind_)
768 ph = 0.5_lkind * real(lat_bnd(j)+lat_bnd(j+1), fms_amip_interp_kind_)
772 ts = real(teq, fms_amip_interp_kind_) - real(tdif, fms_amip_interp_kind_)*sph2 - eps*sph
778 where ( sst < real(tice_crit_k, fms_amip_interp_kind_) )
780 sst = real(tice_crit_k, fms_amip_interp_kind_)
784 end subroutine zonal_sst_
integer function mpp_pe()
Returns processor ID.