346 lines
12 KiB
Fortran
346 lines
12 KiB
Fortran
program simulate
|
|
|
|
Use parkind1, Only : jpim, jprb, jplm
|
|
Use rttov_types, Only : &
|
|
& rttov_coefs, &
|
|
& rttov_options, &
|
|
& rttov_options_scatt, &
|
|
& rttov_scatt_coef, &
|
|
& rttov_chanprof
|
|
implicit none
|
|
!USE MPI
|
|
INCLUDE 'mpif.h'
|
|
|
|
#include "rttov_init_coefs.interface"
|
|
#include "rttov_read_coefs.interface"
|
|
#include "rttov_dealloc_coefs.interface"
|
|
#include "rttov_errorreport.interface"
|
|
#include "rttov_read_scattcoeffs.interface"
|
|
#include "rttov_dealloc_scattcoeffs.interface"
|
|
|
|
integer(kind=4) :: ilevel,iprofile,nprofile,nsurf,nvar
|
|
integer :: itime
|
|
real(kind=4),allocatable,dimension(:) :: pressure,T,q,cfrac,CO2,O3,CH4,N2O,CO,rawdata
|
|
real(kind=4),allocatable,dimension(:,:) :: cloud
|
|
real(kind=4) :: lat0,lon0,psfc
|
|
integer(kind=4) :: landseamask,qcmark
|
|
real(kind=4) :: seaice,clmnwp
|
|
integer :: i,j,k,jch
|
|
real(kind=4) :: t2m,q2m,u10m,v10m,ts,dem
|
|
character*400 :: bkgfile
|
|
Real(Kind=jprb) :: zenith
|
|
Real(Kind=jprb) :: azimuth
|
|
Real(Kind=jprb) :: sunzang
|
|
Real(Kind=jprb) :: sunazang
|
|
Real(Kind=jprb) :: lat
|
|
real(kind=jprb) :: lon
|
|
integer(Kind=jpim) :: year
|
|
integer(Kind=jpim) :: month
|
|
integer(Kind=jpim) :: day
|
|
Integer(Kind=jpim) :: nlevel
|
|
real :: nrecord
|
|
real(kind=4) :: v1,v2,resolution
|
|
Integer(Kind=jpim) :: nchannels,nchan,channelbegin,nchanprof,channelend,channelbin
|
|
integer :: x,y
|
|
integer :: platformid,satid,insid
|
|
integer :: wnposition
|
|
character(len=400) :: coefffile
|
|
Character(len=400) :: emisfilepath
|
|
Character(len=400) :: zeemancoef_file
|
|
Character(len=400) :: scattcoef_file
|
|
Character(len=400) :: temppath,outpath,outfile,latlonfile,outfilename,latlonfilename
|
|
Character (len=10) :: spectype
|
|
Character*10 :: file_id
|
|
Real(kind=4),allocatable,dimension(:) :: BT,obs
|
|
integer(kind=4),allocatable,dimension(:) :: quality
|
|
real(kind=4),allocatable,dimension(:) :: emissivity
|
|
integer :: narg=0,mpinum,mpiid,mpirow,ierror,s_row,e_row,stat
|
|
integer,allocatable,dimension(:) :: pronum
|
|
logical,dimension(2) :: lexist
|
|
character*300 :: configfile
|
|
character*10 :: str
|
|
real(kind=8) :: time
|
|
|
|
Type(rttov_options) :: opts ! Options structure
|
|
Type(rttov_coefs) :: coefs ! Coefficients structure
|
|
Type(rttov_chanprof), Allocatable :: chanprof(:)
|
|
Integer(Kind=jpim) :: errorstatus
|
|
Integer(Kind=jpim), Allocatable :: input_chan(:)
|
|
|
|
namelist /L1info/ nchannels,channelbegin,platformid,satid,insid,coefffile,scattcoef_file,emisfilepath,zeemancoef_file
|
|
namelist /nwpfileinfo/ spectype,bkgfile,nprofile,nlevel
|
|
namelist /output/ temppath,outpath,outfilename
|
|
|
|
narg=IARGC()
|
|
if (narg==0) then
|
|
print* , 'Must have configure file'
|
|
stop
|
|
end if
|
|
CALL get_command_argument(1, configfile)
|
|
open(1,File=configfile,Status='Old')
|
|
Read(1,NML=L1info)
|
|
Read(1,NML=nwpfileinfo)
|
|
Read(1,NML=output)
|
|
close(1)
|
|
!inquire(file=trim(outfile),EXIST=lexist)
|
|
|
|
|
|
call MPI_INIT(ierror)
|
|
call MPI_COMM_SIZE(MPI_COMM_WORLD,mpinum,ierror)
|
|
call MPI_COMM_RANK(MPI_COMM_WORLD,mpiid,ierror)
|
|
|
|
nvar = 20+3*nlevel+7*(nlevel-1) + nchannels + nchannels + 4
|
|
|
|
open(100,file=trim(bkgfile),access='DIRECT',recl=nvar,status='OLD')
|
|
|
|
allocate(pressure(nlevel))
|
|
allocate(CO2(nlevel))
|
|
allocate(N2O(nlevel))
|
|
allocate(CO(nlevel))
|
|
allocate(CH4(nlevel))
|
|
allocate(rawdata(nvar))
|
|
|
|
read(100,rec=1)rawdata
|
|
|
|
nrecord = int(rawdata(1))
|
|
year = int(rawdata(2))
|
|
month = int(rawdata(3))
|
|
day = int(rawdata(4))
|
|
CO2(:) = rawdata(4+1:4+nlevel)
|
|
CH4(:) = rawdata(4+nlevel+1:4+2*nlevel)
|
|
N2O(:) = rawdata(4+2*nlevel+1:4+3*nlevel)
|
|
CO(:) = rawdata(4+3*nlevel+1:4+4*nlevel)
|
|
pressure(:) = rawdata(4+4*nlevel+1:4+5*nlevel)
|
|
|
|
nprofile = nrecord
|
|
allocate(pronum(mpinum))
|
|
pronum(:)=nprofile/mpinum
|
|
do i=1,mpinum
|
|
if (mod(nprofile,mpinum) >0 .and. i<=mod(nprofile,mpinum)) then
|
|
pronum(i)=pronum(i)+1
|
|
end if
|
|
enddo
|
|
mpirow=sum(pronum)
|
|
s_row=sum(pronum(1:mpiid+1))-pronum(mpiid+1)+1
|
|
e_row=sum(pronum(1:mpiid+1))
|
|
if (e_row>nprofile) then
|
|
e_row=nprofile
|
|
end if
|
|
|
|
write(file_id,'(I5)')mpiid+1
|
|
outfile=trim(adjustl(temppath))//'/temp'//trim(adjustl(file_id))//'_'//trim(adjustl(outfilename))
|
|
print*,'===',trim(outfile)
|
|
|
|
open(mpiid+300,file=trim(outfile),access='DIRECT',recl=nchannels*2,status='REPLACE')
|
|
!print*,'s_row=',s_row,'e_row=',e_row,'mpiid=',mpiid
|
|
|
|
if(mpiid==0)then
|
|
print*,'the number of record is ',nrecord
|
|
print*,'pronum=',pronum
|
|
endif
|
|
|
|
!============================== read coef ==========================================
|
|
opts % interpolation % addinterp = .TRUE. ! Allow interpolation of input profile
|
|
opts % interpolation % interp_mode = 1 ! Set interpolation method
|
|
|
|
opts % rt_all % addrefrac = .TRUE. ! Include refraction in path calc
|
|
opts % rt_ir % addaerosl = .FALSE. ! Don't include aerosol effects
|
|
opts % rt_ir % addclouds = .TRUE. ! Include cloud effects
|
|
opts % rt_ir % ir_sea_emis_model = 2
|
|
opts % rt_ir % ir_scatt_model = 2 ! Scattering model for emission source term:
|
|
! 1 => DOM; 2 => Chou-scaling
|
|
opts % rt_ir % vis_scatt_model = 2 ! Scattering model for solar source term:
|
|
! 1 => DOM; 2 => single-scattering
|
|
opts % rt_ir % grid_box_avg_cloud = .TRUE.
|
|
|
|
opts % config % apply_reg_limits = .TRUE.
|
|
opts % config % do_checkinput = .TRUE.
|
|
opts % config % verbose = .TRUE. ! Enable printing of warnings
|
|
|
|
|
|
!======================mixed gas mixing ratio(ppmv)===================================
|
|
|
|
|
|
! channelbegin = 1_jpim
|
|
channelbin = nchannels
|
|
Allocate(chanprof(channelbin)) ! Note these array sizes nchan can vary per profile
|
|
Allocate(input_chan(channelbin))
|
|
|
|
DO jch = channelbegin, channelbegin+channelbin-1
|
|
chanprof(jch-channelbegin+1)%prof = 1_jpim
|
|
chanprof(jch-channelbegin+1)%chan = jch-channelbegin+1 !input_chan(jch)
|
|
input_chan(jch-channelbegin+1) = jch
|
|
End Do
|
|
call rttov_read_coefs (errorstatus, coefs, opts, &
|
|
channels = input_chan(1:channelbin), &
|
|
channels_rec = input_chan(1:channelbin),&
|
|
file_coef=trim(coefffile), &
|
|
file_sccld=trim(scattcoef_file))
|
|
|
|
! call rttov_init_coefs (errorstatus, opts, coefs)
|
|
! print*,'~~~~~simulate~~~~~~~~~~~~~',coefs%coef%fmv_chn,coefs%coef%fmv_ori_nchn,coefs%coef%ff_cwn
|
|
|
|
do iprofile=s_row,e_row !5307,5307 572,572!4300,4300, s_row,e_row
|
|
!print*,'iprofile=',iprofile
|
|
!do iprofile=500,505
|
|
!do iprofile=500,1005
|
|
!do iprofile=1000,1005
|
|
!do iprofile=1,10
|
|
allocate(T(nlevel))
|
|
allocate(q(nlevel))
|
|
allocate(cloud(6,nlevel-1))
|
|
allocate(cfrac(nlevel))
|
|
allocate(O3(nlevel))
|
|
allocate(BT(nchannels))
|
|
allocate(quality(nchannels))
|
|
allocate(obs(nchannels))
|
|
allocate(emissivity(nchannels))
|
|
read(100,rec=iprofile+1,IOSTAT=stat)(rawdata(k),k=1,nvar-4),time
|
|
!print*,'iprofile=',iprofile
|
|
x = int(rawdata(1 ))
|
|
y = int(rawdata(2 ))
|
|
lat = rawdata(3 )
|
|
lon = rawdata(4 )
|
|
azimuth = rawdata(5 )
|
|
zenith = rawdata(6 )
|
|
sunazang = rawdata(7 )
|
|
sunzang = rawdata(8 )
|
|
psfc = rawdata(9 )
|
|
q2m = rawdata(10)
|
|
t2m = rawdata(11)
|
|
u10m = rawdata(12)
|
|
v10m = rawdata(13)
|
|
dem = rawdata(14)
|
|
seaice = rawdata(15)
|
|
landseamask = int(rawdata(16))
|
|
ts = rawdata(17)
|
|
T(:) = rawdata(17+1:17+nlevel)
|
|
q(:) = rawdata(17+nlevel+1:17+2*nlevel)
|
|
O3(:) = rawdata(17+2*nlevel+1:17+3*nlevel)
|
|
qcmark = int(rawdata(17+3*nlevel+1))
|
|
do j=1,6
|
|
cloud(j,:) = rawdata(17+3*nlevel+2+(j-1)*(nlevel-1):17+3*nlevel+1+j*(nlevel-1))
|
|
enddo
|
|
cfrac(:) = rawdata(17+3*nlevel+6*(nlevel-1)+1+1:17+3*nlevel+7*(nlevel-1)+1)
|
|
obs(:) = rawdata(17+3*nlevel+7*(nlevel-1)+1+1:17+3*nlevel+7*(nlevel-1)+1+nchannels)
|
|
clmnwp = rawdata(17+3*nlevel+7*(nlevel-1)+1+nchannels+2)
|
|
emissivity(:) = rawdata(17+3*nlevel+7*(nlevel-1)+1+nchannels+2+1:17+3*nlevel+7*(nlevel-1)+1+nchannels+2+7)
|
|
|
|
! print*,emissivity(:)
|
|
! print*,cloud(1,:)
|
|
! print*,cloud(2,:)
|
|
! print*,cloud(3,:)
|
|
! print*,cloud(4,:)
|
|
! print*,cloud(5,:)
|
|
! print*,cloud(6,:)
|
|
! print*,'landseamask',landseamask
|
|
! print*,'x,y',x,y
|
|
! print*,'T',T
|
|
! print*,'q',q
|
|
! print*,'O3',O3
|
|
! print*,'psfc',psfc
|
|
! print*,'sunzang',sunzang
|
|
! print*,'sunazang',sunazang
|
|
! print*,'zenith',zenith
|
|
! print*,'azimuth',azimuth
|
|
! print*,'q2m',q2m
|
|
! print*,'dem',dem
|
|
! print*,cfrac
|
|
! ! print*,indexid
|
|
! print*,'lat=',lat,'lon=',lon
|
|
! print*,'clmnwp=',clmnwp
|
|
call RTTOV_fwd( nchannels , &
|
|
channelbegin , &
|
|
platformid , &
|
|
satid , &
|
|
insid , &
|
|
coefffile , &
|
|
scattcoef_file , &
|
|
emisfilepath , &
|
|
zeemancoef_file, &
|
|
spectype , &
|
|
nlevel , &
|
|
T , &
|
|
q , &
|
|
cloud , &
|
|
cfrac , &
|
|
CO2 , &
|
|
O3 , &
|
|
CH4 , &
|
|
N2O , &
|
|
CO , &
|
|
t2m , &
|
|
q2m , &
|
|
u10m , &
|
|
v10m , &
|
|
ts , &
|
|
psfc , &
|
|
dem , &
|
|
pressure , &
|
|
lat , &
|
|
lon , &
|
|
landseamask , &
|
|
year , &
|
|
month , &
|
|
zenith , &
|
|
azimuth , &
|
|
sunzang , &
|
|
sunazang , &
|
|
emissivity(:) , &
|
|
coefs , &
|
|
quality(:) , &
|
|
BT(:) )
|
|
! if(mpiid==0)then
|
|
! print*,'+++','iprofile=',iprofile
|
|
! endif
|
|
|
|
! print*,'quality=',quality
|
|
! print*,BT(1:5)
|
|
if(minval(BT)==0)then
|
|
print*,'BT=0'
|
|
BT(:) = -9999.9
|
|
endif
|
|
|
|
! print*,'BT=',BT(:)
|
|
! print*,'mpiid=',mpiid,'iprofile=',iprofile
|
|
! stop
|
|
|
|
write(mpiid+300,rec=iprofile-s_row+1)BT,quality
|
|
|
|
! if(BT(1)<100)then
|
|
! print*,'psfc=',psfc
|
|
! endif
|
|
deallocate(T)
|
|
deallocate(q)
|
|
deallocate(cloud)
|
|
deallocate(cfrac)
|
|
deallocate(O3)
|
|
deallocate(BT)
|
|
deallocate(quality)
|
|
deallocate(obs)
|
|
deallocate(emissivity)
|
|
enddo
|
|
!print*,'the number of exe is',mpiid
|
|
|
|
close(mpiid+300)
|
|
|
|
deallocate(pressure)
|
|
deallocate(pronum)
|
|
deallocate(CO2)
|
|
deallocate(N2O)
|
|
deallocate(CO)
|
|
deallocate(CH4)
|
|
deallocate(rawdata)
|
|
close(100)
|
|
call rttov_dealloc_coefs(errorstatus,coefs)
|
|
deallocate(chanprof)
|
|
deallocate(input_chan)
|
|
call mpi_finalize(ierror)
|
|
|
|
!print*,'finished'
|
|
end program simulate
|
|
|
|
|
|
|
|
|