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