348 lines
12 KiB
Fortran
348 lines
12 KiB
Fortran
module module_satalite
|
|
use HDF5
|
|
use netcdf
|
|
implicit none
|
|
!include 'netcdf.inc'
|
|
|
|
public :: init_satalite,read_satalite
|
|
public :: allocate_satalite,deallocate_satalite
|
|
|
|
integer,public :: nscanl,nscanp,nchannel
|
|
integer,public,parameter :: allocated=1,deallocated=0,clear=3,cloudy=0
|
|
|
|
integer,public,allocatable :: landmask(:,:)
|
|
real(Kind=4),public,allocatable :: satzenith(:,:)
|
|
real(Kind=4),public,allocatable :: elevation(:,:) !unit: m
|
|
real(Kind=4),public,allocatable :: satazimuth(:,:)
|
|
real(Kind=4),public,allocatable :: lat(:,:)
|
|
real(kind=4),public,allocatable :: lon(:,:)
|
|
real(Kind=4),public,allocatable :: sunzenith(:,:)
|
|
real(Kind=4),public,allocatable :: sunazimuth(:,:)
|
|
real(Kind=4),public,allocatable :: obs(:,:,:)
|
|
real(Kind=8),public,allocatable :: time(:,:)
|
|
real(kind=4),public,allocatable :: cloudmask(:,:)
|
|
|
|
integer,private :: hdferror,alocstatus
|
|
integer(HID_T),private :: file_id
|
|
integer(HID_T),private :: dataset_id,dspace_id
|
|
integer(HSIZE_T),private :: maxdimsd(2),maxdimst(3)
|
|
integer(HSIZE_T),private,dimension(2) :: dims_geo
|
|
integer(HSIZE_T),private,dimension(3) :: dims_3d
|
|
integer,private :: ncid,ncerror,dimid,varid
|
|
contains
|
|
|
|
!======================================================================
|
|
subroutine init_satalite(satfilename)
|
|
implicit none
|
|
character*400 :: satfilename
|
|
|
|
call h5open_f(hdferror)
|
|
call h5fopen_f(trim(satfilename), H5F_ACC_RDONLY_F, file_id, hdferror)
|
|
print*,'satfile=',trim(satfilename)
|
|
|
|
call h5dopen_f(file_id, 'ES_RealLW',dataset_id, hdferror)
|
|
call h5dget_space_f(dataset_id, dspace_id, hdferror)
|
|
call h5sget_simple_extent_dims_f(dspace_id, dims_3d, maxdimsd,hdferror)
|
|
|
|
nscanp = dims_3d(1)
|
|
nscanl = dims_3d(2)
|
|
!nchannel = dims_3d(3) !modifed by zpx
|
|
nchannel=721
|
|
|
|
dims_geo(1) = nscanp
|
|
dims_geo(2) = nscanl
|
|
|
|
print*,'nscanp=',nscanp,'nscanl=',nscanl,'nchannel=',nchannel
|
|
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
call h5fclose_f(file_id,hdferror)
|
|
call h5close_f(hdferror)
|
|
|
|
return
|
|
end subroutine
|
|
|
|
!=======================================================================
|
|
subroutine allocate_satalite(alocstatus)
|
|
implicit none
|
|
integer(kind=4) :: alocstatus
|
|
|
|
allocate(satzenith(nscanp,nscanl))
|
|
allocate(satazimuth(nscanp, nscanl))
|
|
allocate(lat(nscanp, nscanl))
|
|
allocate(lon(nscanp, nscanl))
|
|
allocate(elevation(nscanp, nscanl))
|
|
allocate(sunzenith(nscanp, nscanl))
|
|
allocate(sunazimuth(nscanp, nscanl))
|
|
allocate(obs(nscanp, nscanl,nchannel))
|
|
allocate(time(nscanp, nscanl))
|
|
allocate(landmask(nscanp, nscanl))
|
|
allocate(cloudmask(nscanp, nscanl))
|
|
alocstatus=allocated
|
|
return
|
|
end subroutine allocate_satalite
|
|
|
|
!=====================================================================
|
|
subroutine deallocate_satalite(alocstatus)
|
|
implicit none
|
|
integer(kind=4) :: alocstatus
|
|
|
|
deallocate(satzenith)
|
|
deallocate(satazimuth)
|
|
deallocate(lat)
|
|
deallocate(lon)
|
|
deallocate(sunzenith)
|
|
deallocate(sunazimuth)
|
|
deallocate(elevation)
|
|
deallocate(obs)
|
|
deallocate(time)
|
|
deallocate(landmask)
|
|
deallocate(cloudmask)
|
|
alocstatus=deallocated
|
|
return
|
|
end subroutine deallocate_satalite
|
|
|
|
!========================================================================================================
|
|
subroutine read_satalite(lsmfilename,elefilename,geofilename,satfilename,obsfilename,clmfilename)
|
|
implicit none
|
|
|
|
character*400 :: lsmfilename,elefilename,satfilename,geofilename,obsfilename,clmfilename
|
|
character*20 :: times
|
|
integer :: idstart
|
|
integer :: second
|
|
integer(kind=8) :: timet
|
|
integer :: i,j,k,ncal
|
|
integer :: ny,nx,nxelev,nyelev
|
|
real(kind=4) :: scale,offset
|
|
real(kind=4) :: linetime
|
|
real(kind=4),allocatable :: landseamask(:,:)
|
|
real(kind=4),allocatable :: zsfc(:,:)
|
|
integer(kind=4),allocatable :: index_lone(:,:),index_late(:,:)
|
|
integer(kind=4),allocatable :: index_lonlsm(:,:),index_latlsm(:,:)
|
|
integer(kind=4),allocatable :: flag_universe(:,:)
|
|
real(kind=4),allocatable :: lon180(:,:)
|
|
integer(HSIZE_T),dimension(1) :: dims_1d
|
|
real(kind=4),allocatable :: waven(:)
|
|
real(kind=4),allocatable :: NOM(:,:,:)!,ES(:,:,:)
|
|
real(kind=4) :: Hanming_a=0.23
|
|
real(kind=4) :: w0,w1,w2,ES
|
|
integer(HSIZE_T) :: maxdims(1)
|
|
|
|
|
|
allocate(index_lone(nscanp,nscanl))
|
|
allocate(index_late(nscanp,nscanl))
|
|
allocate(index_lonlsm(nscanp,nscanl))
|
|
allocate(index_latlsm(nscanp,nscanl))
|
|
allocate(flag_universe(nscanp,nscanl))
|
|
allocate(lon180(nscanp,nscanl))
|
|
|
|
flag_universe = 0
|
|
|
|
!============================= geo ==================================
|
|
|
|
call h5open_f(hdferror)
|
|
call h5fopen_f(trim(satfilename), H5F_ACC_RDONLY_F, file_id, hdferror)
|
|
|
|
call h5dopen_f(file_id, 'IRLW_SatelliteZenith',dataset_id, hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_REAL,satzenith,dims_geo,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
|
|
call h5dopen_f(file_id,'IRLW_SatelliteAzimuth',dataset_id,hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_REAL,satazimuth,dims_geo,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
|
|
call h5dopen_f(file_id,'IRLW_SolarZenith',dataset_id,hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_REAL,sunzenith,dims_geo,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
|
|
call h5dopen_f(file_id,'IRLW_SolarAzimuth',dataset_id,hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_REAL,sunazimuth,dims_geo,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
|
|
call h5dopen_f(file_id,'TimeStamp',dataset_id,hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_DOUBLE,time,dims_geo,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
print*,'time=',time(1,1:5)
|
|
|
|
call h5dopen_f(file_id,'IRLW_Latitude',dataset_id,hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_REAL,lat,dims_geo,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
! print*,'lat=',lat(1:20,1)
|
|
call h5dopen_f(file_id,'IRLW_Longitude',dataset_id,hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_REAL,lon,dims_geo,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
lon180 = lon
|
|
where(lon<0 .and. lon>=-180)
|
|
lon = lon+360
|
|
endwhere
|
|
where(lon>360 .or. lon<-180)
|
|
flag_universe=1
|
|
endwhere
|
|
|
|
where(lat>90 .or. lat<-90)
|
|
flag_universe=1
|
|
endwhere
|
|
! print*,'lon=',lon(1,1:5)
|
|
!=========================== observation ========================================
|
|
|
|
allocate(NOM (nscanp,nscanl,nchannel))
|
|
! allocate(ES (nscanp,nscanl,nchannel))
|
|
allocate(waven (nchannel))
|
|
obs = -9999.9
|
|
dims_1d(1) = nchannel
|
|
|
|
call h5dopen_f(file_id,'IRLW_VaildWaveLength',dataset_id,hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_REAL,waven,dims_1d,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
|
|
call h5dopen_f(file_id,'ES_RealLW',dataset_id,hdferror)
|
|
call h5dread_f(dataset_id,H5T_NATIVE_REAL,NOM,dims_3d,hdferror)
|
|
call h5dclose_f(dataset_id,hdferror)
|
|
! print*,'NOM=',NOM(10:50,10:50,352)
|
|
! print*,'waven=',waven(352)!print*,'obs=',obs(10:50,10:50,352)
|
|
call h5fclose_f(file_id,hdferror)
|
|
call h5close_f(hdferror)
|
|
w0 = Hanming_a
|
|
w1 = 1.0-2.0*Hanming_a
|
|
w2 = Hanming_a
|
|
do i=1,nchannel
|
|
do j=1,nscanl
|
|
do k=1,nscanp
|
|
|
|
if (i==1)then
|
|
if (NOM(k,j,i)<0 .or. NOM(k,j,i)>=1000 .or. NOM(k,j,i+1)<0 .or. NOM(k,j,i+1)>=1000)then
|
|
obs(k,j,i) = -9999.9
|
|
else
|
|
ES = w1/(w1+w2)*NOM(k,j,i)+w2/(w1+w2)*NOM(k,j,i+1)
|
|
call radiance_to_BT(dble(ES),dble(waven(i)),obs(k,j,i))
|
|
endif
|
|
end if
|
|
|
|
if(i>1 .and. i<nchannel)then
|
|
if (NOM(k,j,i-1)<0 .or. NOM(k,j,i-1)>=1000 .or. NOM(k,j,i)<0 .or. NOM(k,j,i)>=1000 .or. NOM(k,j,i+1)<0 .or. NOM(k,j,i+1)>=1000)then
|
|
obs(k,j,i) = -9999.9
|
|
else
|
|
ES=w0*NOM(k,j,i-1)+w1*NOM(k,j,i)+w2*NOM(k,j,i+1)
|
|
call radiance_to_BT(dble(ES),dble(waven(i)),obs(k,j,i))
|
|
endif
|
|
endif
|
|
|
|
if(i == nchannel)then
|
|
if (NOM(k,j,i-1)<0 .or. NOM(k,j,i-1)>=1000 .or. NOM(k,j,i)<0 .or. NOM(k,j,i)>=1000)then
|
|
obs(k,j,i) = -9999.9
|
|
else
|
|
ES = w1/(w1+w2)*NOM(k,j,i)+w0/(w1+w0)*NOM(k,j,i-1)
|
|
call radiance_to_BT(dble(ES),dble(waven(i)),obs(k,j,i))
|
|
endif
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
print*,'obs=',obs(1,1,1:5)
|
|
|
|
! print*,'obs=',obs(10:50,10:50,352)
|
|
deallocate(NOM )
|
|
! deallocate(ES )
|
|
deallocate(waven )
|
|
!========================== landseamask =====================================
|
|
ncerror = nf90_open(trim(lsmfilename), nf90_nowrite, ncid)
|
|
ncerror = nf90_inq_dimid(ncid,'longitude',dimid)
|
|
ncerror = nf90_inquire_dimension(ncid,dimid,len = nx)
|
|
ncerror = nf90_inq_dimid(ncid,'latitude',dimid)
|
|
ncerror = nf90_inquire_dimension(ncid,dimid,len = ny)
|
|
|
|
allocate(landseamask(nx,ny))
|
|
ncerror = nf90_open(trim(lsmfilename), nf90_nowrite, ncid)
|
|
ncerror = nf90_inq_varid(ncid, 'lsm', varid)
|
|
ncerror = nf90_get_var(ncid,varid, landseamask)
|
|
ncerror = nf90_get_att(ncid,varid,'scale_factor',scale)
|
|
ncerror = nf90_get_att(ncid,varid,'add_offset',offset)
|
|
where(landseamask/=-32767)
|
|
landseamask = landseamask*scale+offset
|
|
endwhere
|
|
where(landseamask==-32767)
|
|
landseamask = -9999.9
|
|
endwhere
|
|
where(landseamask==0)
|
|
landseamask = 3
|
|
endwhere
|
|
where(landseamask>0 .and. landseamask<1)
|
|
landseamask = 2
|
|
endwhere
|
|
ncerror=nf90_close(ncid)
|
|
index_lonlsm = anint(lon/0.125+1)
|
|
index_latlsm = anint((90-lat)/0.125+1)
|
|
|
|
do i=1,nscanl
|
|
do j=1,nscanp
|
|
if (flag_universe(j,i)==0)then
|
|
landmask(j,i) = landseamask(index_lonlsm(j,i),index_latlsm(j,i))
|
|
else
|
|
landmask(j,i) = 4
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
deallocate(landseamask)
|
|
!======================== elevation =========================================
|
|
!elevation data: resolution=0.05,lat0=-89.975,lon0=-179.975)
|
|
ncerror = nf90_open(trim(elefilename), nf90_nowrite, ncid)
|
|
ncerror = nf90_inq_dimid(ncid,'lon',dimid)
|
|
ncerror = nf90_inquire_dimension(ncid,dimid,len = nxelev)
|
|
ncerror = nf90_inq_dimid(ncid,'lat',dimid)
|
|
ncerror = nf90_inquire_dimension(ncid,dimid,len = nyelev)
|
|
|
|
allocate(zsfc(nxelev,nyelev))
|
|
ncerror = nf90_open(trim(elefilename), nf90_nowrite, ncid)
|
|
ncerror = nf90_inq_varid(ncid, 'Band1', varid)
|
|
ncerror = nf90_get_var(ncid, varid, zsfc)
|
|
ncerror = nf90_close(ncid)
|
|
index_lone = anint((lon180+179.975)/0.05+1)
|
|
index_late = anint((lat+89.975)/0.05+1)
|
|
do i=1,nscanl
|
|
do j=1,nscanp
|
|
if (flag_universe(j,i)==0)then
|
|
elevation(j,i) = zsfc(index_lone(j,i),index_late(j,i))
|
|
else
|
|
elevation(j,i) = -9999.9
|
|
endif
|
|
enddo
|
|
enddo
|
|
deallocate(zsfc)
|
|
!============================= cloud mask ================================
|
|
if(trim(adjustl(clmfilename))/='None')then
|
|
ncerror = nf90_open(trim(clmfilename), nf90_nowrite, ncid)
|
|
ncerror = nf90_inq_varid(ncid, 'CLM', varid)
|
|
ncerror = nf90_get_var(ncid,varid, cloudmask)
|
|
ncerror = nf90_close(ncid)
|
|
else
|
|
cloudmask = cloudy
|
|
endif
|
|
print*,'-----------------After read satelite!--------------------------'
|
|
!========================== deallocate =====================================
|
|
|
|
deallocate(index_lone)
|
|
deallocate(index_late)
|
|
deallocate(index_lonlsm)
|
|
deallocate(index_latlsm)
|
|
deallocate(flag_universe)
|
|
deallocate(lon180)
|
|
! stop
|
|
return
|
|
end subroutine read_satalite
|
|
|
|
subroutine radiance_to_BT(radiance,waven,BT)
|
|
|
|
implicit none
|
|
real(kind=8) :: c_lights = 2.997924 !e8
|
|
real(kind=8) :: h_plank = 6.6262 !e-34
|
|
real(kind=8) :: k_bolzm = 1.3805 !e-23
|
|
real(kind=8) :: waven,radiance
|
|
real(kind=8) :: wavenumber
|
|
real(kind=4) :: BT
|
|
wavenumber = waven*100.0
|
|
radiance = radiance*wavenumber*wavenumber*1e-5
|
|
BT = real(h_plank*c_lights*wavenumber*1e-3/(k_bolzm*log(2*h_plank*c_lights**2*1e-3*(wavenumber*1e-3)**5/radiance+1)))
|
|
return
|
|
end subroutine radiance_to_BT
|
|
|
|
end module module_satalite |