subroutine get_seaice(xlats_ij,xlons_ij,ny,nx,sice)
!
! Abstract: get sea ice analysis at atime and target resolution (nx,ny)
!
!
! input & output
!
real, dimension(nx*ny), intent(in) :: xlats_ij ! latitudes of target grids (nx*ny)
real, dimension(nx*ny), intent(in) :: xlons_ij ! latitudes of target grids (nx*ny)
real, dimension(nx,ny), intent(out) :: sice ! sea ice atime (nx,ny)
integer, intent(in) :: nx,ny
! local declare
integer, dimension(nx*ny) :: mask_ij ! mask at target grids (0 = water, 1 = land, 2 = sea ice) (nx*ny)
real, dimension(nx*ny) :: sice_ij ! sea ice analysis at target grids (nx*ny)
real, allocatable, dimension(:,:) :: sice0 ! sea ice analysis
integer, allocatable, dimension(:,:) :: smask ! surface mask of sea ice
real, allocatable, dimension(:) :: sxlats ! latitudes of sea ice analysis
real, allocatable, dimension(:) :: sxlons ! longitudes of sea ice analysis
integer :: nxs,nys,mon1,mon2,sfcflag,i,j
character (len=6), parameter :: fin_iceanl='iceanl' ! sea ice analysis file name
!
! get the dimensions of the sea ice analysis & allocate the related arrays
!
call get_dim_grb(fin_iceanl,nys,nxs)
allocate( sice0(nxs,nys),smask(nxs,nys),sxlats(nys),sxlons(nxs) )
!
! read sea ice analysis
!
call read_seaice_grb(fin_iceanl,nys,nxs,sxlats,sxlons,sice0)
!
! get sice (nx by ny lat/lon)
!
if ( nx == nxs .and. ny == nys ) then
sice(:,:) = sice0(:,:)
write(*,'(a,4I8,2F9.3)') 'same dimensions,nx,ny,nxs,nys : ',nx,ny,nxs,nys,minval(sice),maxval(sice)
else
sfcflag=0
smask=0
mask_ij=0
write(*,'(a,4I8)') 'different dimensions,nx,ny,nxs,nys : ',nx,ny,nxs,nys
call lalo_to_tile(sice0, smask, sxlats, sxlons, nys, nxs, &
sice_ij, mask_ij, xlats_ij,xlons_ij,ny, nx, &
sfcflag,0.0,0,0.0)
write(*,'(a,2F9.3)') 'done with lalo_to_tile for sice',minval(sice_ij),maxval(sice_ij)
sice(:,:) = reshape (sice_ij, (/nx,ny/) )
endif
end subroutine get_seaice
subroutine get_dim_grb(file_ice,mlat,mlon)
! . . . .
! abstract: get_grb_dim : get dimensions
! prgmmr: xu li org: np23 date: 2019-03-13
!
!
! program history log:
!
! input argument list:
! file - file name of GRIB SST file
! output
! mlat,mlon
! call subs: getgbh
!
!
!$$$
implicit none
! Declare passed variables and arrays
character (len=*), intent(in) :: file_ice
integer, intent(out) :: mlat,mlon
! Declare local parameters
integer,parameter:: lu = 22 ! FORTRAN unit number of GRIB SST file
integer :: iret,i,j,kf,kg,k
integer, dimension(22):: jgds,kgds
integer, dimension(25):: jpds,kpds
!************+******************************************************************************
!
! Open analysis file (GRIB)
write(*,*) 'get_dim_grb, file_ice : ',file_ice
call baopenr(lu,trim(file_ice),iret)
if (iret /= 0 ) then
write(6,*)'get_dim_grb: ***ERROR*** opening sea ice file',file_ice
stop
endif
! Define variables for read
j=-1
jpds=-1
jgds=-1
jpds(5)=91 ! sea ice variable
jpds(6)=102 ! surface
call getgbh(lu,0,j,jpds,jgds,kg,kf,k,kpds,kgds,iret)
mlat = kgds(3) ! number points on longitude circle (360)
mlon = kgds(2) ! number points on latitude circle (720)
write(*,*) 'sea ice nys, nxs : ',mlat, mlon
call baclose(lu,iret)
if (iret /= 0 ) then
write(6,*)'get_dim_grb: ***ERROR*** close sea ice file'
stop
endif
end subroutine get_dim_grb
subroutine read_seaice_grb(file_sice,mlat,mlon,xlats,xlons,sice)
!$$$ subprogram documentation block
! . . . .
! subprogram: rdgrbsice read grib1 SST analysis
! prgmmr: xu li org: np23 date: 2003-04-04
!
! abstract: read sea ice analysis (GRIB format) and save it as expanded and transposed array
!
! Subroutine rdgrbsice must be compiled with the NCEP W3 library
! and the BACIO library.
!
!
! program history log:
! 2003-04-04 xu li, bert katz
! 2005-04-18 treadon - fill southern and northern rows of sice grid
! with mean of adjacent row. This treatment is
! consistent with rdgesfc.f90 and rdgesig.f90
!
! input argument list:
! file_sice - file name of GRIB sea ice file
! mlat - dimension in y-direction
! mlon - dimension in x-direction
!
! call subs: getgbh, getgb
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
implicit none
! Declare passed variables and arrays
character(6), intent(in) :: file_sice
integer, intent(in) :: mlat,mlon
real, dimension(mlat), intent(out) :: xlats
real, dimension(mlon), intent(out) :: xlons
real, dimension(mlon,mlat), intent(out) :: sice
! Declare local parameters
integer,parameter:: lu_sice = 22 ! FORTRAN unit number of GRIB SST file
real, parameter :: deg2rad = 3.141593/180.0
! Declare local variables and arrays
logical(1), allocatable, dimension(:) :: lb
integer :: nlat_sice,nlon_sice,ilon,jlat
integer :: iret,ni,nj
integer :: mscan,kb1
integer :: jincdir,i,iincdir,kb2,kb3,kf,kg,k,j,jf
integer, dimension(22):: jgds,kgds
integer, dimension(25):: jpds,kpds
real :: x0,y0,dres
real, allocatable, dimension(:) :: f
!************+******************************************************************************
!
! Open SST analysis file (GRIB)
call baopenr(lu_sice,trim(file_sice),iret)
if (iret /= 0 ) then
write(6,*)'RDGRBICE: ***ERROR*** opening Sea Ice file'
stop
endif
! Define SST variables for read
j=-1
jpds=-1
jgds=-1
jpds(5)=91 ! Sea Ice fraction (1=ice; 0=no ice)
jpds(6)=102 ! surface
call getgbh(lu_sice,0,j,jpds,jgds,kg,kf,k,kpds,kgds,iret)
nlat_sice = kgds(3) ! number points on longitude circle (360)
nlon_sice = kgds(2) ! number points on latitude circle (720)
if ( mlat /= nlat_sice .or. mlon /= nlon_sice ) then
write(*,*) 'STOP inconsistent dimensions : ',nlat_sice, nlon_sice,mlat,mlon
stop
endif
!
! get xlats and xlons
!
dres = 180.0/real(mlat)
y0 = 0.5*dres-90.0
x0 = 0.5*dres
! Get lat_sst & lon_sst
do jlat = 1, mlat
xlats(jlat) = y0 + real(jlat-1)*dres
enddo
do ilon = 1, mlon
xlons(ilon) = (x0 + real(ilon-1)*dres)
enddo
! Allocate arrays
allocate(lb(nlat_sice*nlon_sice))
allocate(f(nlat_sice*nlon_sice))
jf=nlat_sice*nlon_sice
! Read in the analysis
call getgb(lu_sice,0,jf,j,jpds,jgds,kf,k,kpds,kgds,lb,f,iret)
if (iret /= 0) then
write(6,*)'RDGRBSICE: ***ERROR*** reading sice analysis data record'
deallocate(lb,f)
stop
endif
! Load dimensions and grid specs. Check for unusual values
ni=kgds(2) ! 720 for 0.5 x 0.5
nj=kgds(3) ! 360 for 0.5 x 0.5 resolution
mscan=kgds(11)
kb1=ibits(mscan,7,1) ! i scan direction
kb2=ibits(mscan,6,1) ! j scan direction
kb3=ibits(mscan,5,1) ! (i,j) or (j,i)
! Get i and j scanning directions from kb1 and kb2.
! 0 yields +1, 1 yields -1. +1 is west to east, -1 is east to west.
iincdir = 1-kb1*2
! 0 yields -1, 1 yields +1. +1 is south to north, -1 is north to south.
jincdir = kb2*2 - 1
do k=1,kf
! kb3 from scan mode indicates if i points are consecutive
! or if j points are consecutive
if(kb3==0)then ! (i,j)
i=(ni+1)*kb1+(mod(k-1,ni)+1)*iincdir
j=(nj+1)*(1-kb2)+jincdir*((k-1)/ni+1)
else ! (j,i)
j=(nj+1)*(1-kb2)+(mod(k-1,nj)+1)*jincdir
i=(ni+1)*kb1+iincdir*((k-1)/nj+1)
endif
sice(i,j)=f(k)
end do
deallocate(lb,f)
end subroutine read_seaice_grb