subroutine rdgrbsst(file_sst,mlat_sst,mlon_sst,&
sst_an,rlats_sst,rlons_sst,nlat_sst,nlon_sst)
!$$$ subprogram documentation block
! . . . .
! subprogram: rdgrbsst read SST analysis
! prgmmr: xu li org: np23 date: 2003-04-04
!
! abstract: read SST analysis (GRIB format) and save it as expanded and transposed array
!
! Subroutine rdgrbsst 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 sst grid
! with mean of adjacent row. This treatment is
! consistent with rdgesfc.f90 and rdgesig.f90
!
! input argument list:
! file_sst - file name of GRIB SST file
! mlat_sst,mlon_sst
!
! argument list defined by this reading:
! sst_an - SST field (for 0.5 x 0.5 resolution)
! sst_an(1,1) is at 90.00 deg. S, 0.25 deg. W (-90.00,-0.25)
! sst_an(1,2) is at 90.00 deg. S, 0.25 deg. E (-90.00,+0.25)
! sst_an(2,1) is at 89.75 deg. S, 0.25 deg. W (-89.75,-0.25)
! sst_an(2,2) is at 89.75 deg. S, 0.25 deg. E (-89.75,+0.25)
! sst_an(361,721) is at 89.75 deg. N, 0.25 deg. W (+89.75,359.75)
! sst_an(361,722) is at 89.75 deg. N, 0.25 deg. E (+89.75,360.25)
! sst_an(362,721) is at 90.00 deg. N, 0.25 deg. W (+90.00,359.75)
! sst_an(362,722) is at 90.00 deg. N, 0.25 deg. E (+90.00,360.25)
! Note: (1) The data is stored from north to south originally in GRIB format,
! but is stored from south to north with this reading routine
! (2) Two poles are added and their values are extrapolated in meridional direction
! (3) Two ends are added and their values (periodic) in zonal direction
! (4) The output (sst_an) dimension is changed as opposite order, e.g. (for 0.5 x 0.5)
! (722,362) ==> (362, 722)
! nlat_sst - latitudinal dimension of SST
! nlon_sst - longitudinal dimension of SST
! xsst0 - latitude of origin
! ysst0 - longitude of origin
! dres - lat/lon increment
!
! call subs: getgbh, getgb
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: h1000,half,deg2rad,zero
implicit none
! Declare passed variables and arrays
character(6) ,intent(in ) :: file_sst
integer(i_kind) ,intent(in ) :: mlat_sst,mlon_sst
integer(i_kind) ,intent( out) :: nlat_sst,nlon_sst
real(r_kind),dimension(mlat_sst) ,intent( out) :: rlats_sst
real(r_kind),dimension(mlon_sst) ,intent( out) :: rlons_sst
real(r_kind),dimension(mlat_sst,mlon_sst),intent( out) :: sst_an
! Declare local parameters
integer(i_kind),parameter:: lu_sst = 21 ! FORTRAN unit number of GRIB SST file
! Declare local variables and arrays
logical(1), allocatable, dimension(:) :: lb
integer(i_kind) iret,ni,nj
integer(i_kind) iyrstmp,imstmp,idstmp,mscan,kb1
integer(i_kind) jincdir,i,iincdir,kb2,kb3,kf,kg,k,j,jf
integer(i_kind),dimension(22):: jgds,kgds
integer(i_kind),dimension(25):: jpds,kpds
real(r_kind) xsst0,ysst0,dres,sums,sumn
real(r_kind), allocatable, dimension(:) :: f
real(r_kind), allocatable,dimension(:,:) :: sst ! SST analysis (RTG or others), for read
!******************************************************************************************
!
! Open SST analysis file (GRIB)
call baopenr(lu_sst,trim(file_sst),iret)
if (iret /= 0 ) then
write(6,*)'RDGRBSST: ***ERROR*** opening SST file'
call stop2(59)
endif
! Define SST variables for read
j=-1
jpds=-1
jgds=-1
jpds(5)=11 ! SST variable
jpds(6)=1 ! surface
call getgbh(lu_sst,0,j,jpds,jgds,kg,kf,k,kpds,kgds,iret)
nlat_sst = kgds(3) ! number points on longitude circle (360)
nlon_sst = kgds(2) ! number points on latitude circle (720)
! Allocate arrays
allocate(lb(nlat_sst*nlon_sst))
allocate(f(nlat_sst*nlon_sst))
jf=nlat_sst*nlon_sst
! Read in the analysis
call getgb(lu_sst,0,jf,j,jpds,jgds,kf,k,kpds,kgds,lb,f,iret)
if (iret /= 0) then
write(6,*)'RDGRBSST: ***ERROR*** reading sst analysis data record'
deallocate(lb,f)
call stop2(59)
endif
! Allocate SST analysis and the arrays for its dimensions
nlat_sst = nlat_sst + 2 ! Add two poles (90S & 90N)
nlon_sst = nlon_sst + 2 ! Add two buffer ends in zonal direction (-0.25 & 0.25)
if ( (nlat_sst>mlat_sst) .or. (nlon_sst>mlon_sst) ) then
write(6,*)'RDGRBSST: inconsistent dimensions. mlat_sst,mlon_sst=',&
mlat_sst,mlon_sst,' -versus- nlat_sst,nlon_sst=',nlat_sst,nlon_sst
deallocate(lb,f)
call stop2(60)
endif
allocate(sst(nlon_sst,nlat_sst))
xsst0 = -real(kgds(4))/h1000 ! latitude of origin
ysst0 = real(kgds(5))/h1000 ! longitude of origin
dres = real(kgds(9))/h1000 ! increment of latitude/longitude (the same here)
! Get lat_sst & lon_sst
do i = 2, nlat_sst - 1
rlats_sst(i) = (xsst0 + float(i-2)*dres)*deg2rad
enddo
rlats_sst(1) = -90.0_r_kind*deg2rad
rlats_sst(nlat_sst) = 90.0_r_kind*deg2rad
do j = 2, nlon_sst - 1
rlons_sst(j) = (ysst0 + float(j-2)*dres)*deg2rad
enddo
rlons_sst(1) = -half*dres*deg2rad ! 1
rlons_sst(nlon_sst) = (360._r_kind+half*dres)*deg2rad ! 722
! 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
sst(i+1,j+1)=f(k)
end do
! Set sst value at 90 S (1) & 90 N (362 = nlon_sst) to be mean of
! adjacent row
i = 0
sums=zero
sumn=zero
do j = 2, nlon_sst - 1
i = i + 1
sums = sums + sst(j,2)
sumn = sumn + sst(j,nlat_sst-1)
end do
sums = sums / float(i)
sumn = sumn / float(i)
do j = 2,nlon_sst-1
sst(j,1) = sums
sst(j,nlat_sst) = sumn
enddo
! Get sst value for added two ends in zonal direction
do i = 1, nlat_sst
sst(1,i) = sst(nlon_sst-1,i)
sst(nlon_sst,i) = sst(2,i)
enddo
! Transpose sst to sst_an (output)
do j = 1, nlon_sst
do i = 1, nlat_sst
sst_an(i,j) = sst(j,i)
enddo
enddo
iyrstmp=kpds(8) + (kpds(21)-1)*100
imstmp=kpds(9)
idstmp=kpds(10)
deallocate(lb,f)
return
end subroutine rdgrbsst