-
Andi Walther authoredAndi Walther authored
cx_rttov_bridge_mod.f90 19.18 KiB
module cx_rttov_bridge_mod
!use cx_pfaast_constants_mod, only: tstd, ostd, wstd,pstd
! rttov_const contains useful RTTOV constants
USE rttov_const, ONLY : &
errorstatus_success, &
errorstatus_fatal, &
platform_name, &
inst_name
! rttov_types contains definitions of all RTTOV data types
USE rttov_types, ONLY : &
rttov_options, &
rttov_coefs, &
rttov_profile, &
rttov_transmission, &
rttov_radiance, &
rttov_chanprof, &
rttov_emissivity, &
rttov_reflectance
USE parkind1, ONLY : jpim, jprb, jplm
USE rttov_unix_env, ONLY : rttov_exit
IMPLICIT NONE
#include "rttov_direct.interface"
#include "rttov_parallel_direct.interface"
#include "rttov_read_coefs.interface"
#include "rttov_dealloc_coefs.interface"
#include "rttov_alloc_direct.interface"
#include "rttov_user_options_checkinput.interface"
#include "rttov_print_opts.interface"
#include "rttov_print_profile.interface"
#include "rttov_skipcommentline.interface"
!--------------------------
!
INTEGER(KIND=jpim), PARAMETER :: iup = 20 ! unit for input profile file
INTEGER(KIND=jpim), PARAMETER :: ioout = 21 ! unit for output
! RTTOV variables/structures
!====================
TYPE(rttov_options) :: opts ! Options structure
TYPE(rttov_coefs) :: coefs ! Coefficients structure
TYPE(rttov_chanprof), POINTER :: chanprof(:) => NULL() ! Input channel/profile list
LOGICAL(KIND=jplm), POINTER :: calcemis(:) => NULL() ! Flag to indicate calculation of emissivity within RTTOV
TYPE(rttov_emissivity), POINTER :: emissivity(:) => NULL() ! Input/output surface emissivity
LOGICAL(KIND=jplm), POINTER :: calcrefl(:) => NULL() ! Flag to indicate calculation of BRDF within RTTOV
TYPE(rttov_reflectance), POINTER :: reflectance(:) => NULL() ! Input/output surface BRDF
TYPE(rttov_profile), POINTER :: profiles(:) => NULL() ! Input profiles
TYPE(rttov_transmission) :: transmission ! Output transmittances
TYPE(rttov_radiance) :: radiance ! Output radiances
INTEGER(KIND=jpim) :: errorstatus ! Return error status of RTTOV subroutine calls
INTEGER(KIND=jpim) :: alloc_status
CHARACTER(LEN=11) :: NameOfRoutine = 'test_rttovv'
! variables for input
!====================
CHARACTER(LEN=256) :: coef_filename,cld_coef_filename
INTEGER(KIND=jpim) :: nthreads
INTEGER(KIND=jpim) :: dosolar
INTEGER(KIND=jpim) :: nlevels
INTEGER(KIND=jpim) :: nprof INTEGER(KIND=jpim) :: nchannels
INTEGER(KIND=jpim) :: nchanprof
INTEGER(KIND=jpim), ALLOCATABLE :: channel_list(:)
REAL(KIND=jprb) :: trans_out(10)
! loop variables
INTEGER(KIND=jpim) :: j, jch
INTEGER(KIND=jpim) :: np, nch
INTEGER(KIND=jpim) :: ilev, nprint
INTEGER(KIND=jpim) :: iprof, joff
INTEGER :: ios
contains
function channel_map (sensor,ancil_data_path, chn,coef_filename,cld_coef_filename) result (list)
#define STRINGIFY(x) x
implicit none
character(len = *) :: sensor
character(len = * ) :: ancil_data_path
integer :: chn
integer :: list
integer, allocatable:: chn_list(:)
! integer :: mod_list(5) = [20,22,29,31,32]
integer :: kk
character (len = *) :: coef_filename
character (len = *) :: cld_coef_filename
character (len=:), allocatable :: path
character (len = 100) :: sensor_string
character (len=1) :: rttov_version_string
character (len = 1) :: metop_nr
integer :: avhrr_num
#ifndef RTTOVPATH
print*,'RTTOV PATH not set ..stop'
print*,__TIMESTAMP__
! print*,RTTOVPATH
stop
#endif
path = trim(ancil_data_path) // "static/rttov/"
rttov_version_string = '9'
! --
! the mapping translates one channel infot
! RTTOV channel number
!
!
allocate(chn_list(45))
chn_list = -1
! metop
if (sensor .eq. 'AVHRR-METOPA') metop_nr = '2'
if (sensor .eq. 'AVHRR-METOPB') metop_nr = '1'
if (sensor .eq. 'AVHRR-METOPC') metop_nr = '3'
select case(sensor)
case('MTSAT-1') list = chn
sensor_string = 'mtsat_1_imager'
case('MTSAT-2')
list = chn
sensor_string = 'mtsat_2_imager'
case( 'MODIS-AQUA')
!list = chn -19
! if ( chn .gt. 26) list = chn - 20
list = chn
sensor_string = 'eos_2_modis'
case( 'MODIS-TERRA')
!list = chn -19
! if ( chn .gt. 26) list = chn - 20
list = chn
sensor_string = 'eos_1_modis'
case ( 'VIIRS-SNPP')
chn_list = [6,9,3,4,10,14,15,1,2,-1,-1,-1,-1,-1,7,-1,-1,-1,-1,17,-1,18,-1,-1,-1,12,-1,-1,19,-1,20,22]
list = chn_list(chn)
sensor_string = 'jpss_0_viirs'
case ( 'VIIRS-N20')
chn_list = [6,9,3,4,10,14,15,1,2,-1,-1,-1,-1,-1,7,-1,-1,-1,-1,17,-1,18,-1,-1,-1,12,-1,-1,19,-1,20,22]
list = chn_list(chn)
sensor_string = 'noaa_20_viirs'
case ('AHI8')
chn_list(20:38) = [7,-1,-1,-1,-1,-1,-1,9,10,11,12,14,15,16,-1,-1,-1,8,13]
list = chn_list(chn)
sensor_string = 'himawari_8_ahi'
case ('AHI9')
chn_list(20:38) = [7,-1,-1,-1,-1,-1,-1,9,10,11,12,14,15,16,-1,-1,-1,8,13]
list = chn_list(chn)
sensor_string = 'himawari_9_ahi'
case ('GOES-16')
chn_list = [2,3,1,-1,-1,5,6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,7,-1,7,-1,-1,-1,4,9,10,11,12,14,15,16,-1,-1,-1,8,13,-1,-1,-1,-1,-1,-1,-1]
list = chn_list(chn)
sensor_string = 'goes_16_abi'
case ('GOES-17')
chn_list = [2,3,1,-1,-1,5,6,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,7,-1,7,-1,-1,-1,4,9,10,11,12,14,15,16,-1,-1,-1,8,13,-1,-1,-1,-1,-1,-1,-1]
list = chn_list(chn)
sensor_string = 'goes_17_abi'
case ( 'SEVIRI-MSG08')
!---stw chn_list(20:33) = [1,-1,-1,-1,-1,-1,-1,2,3,4,5,6,7,8]
chn_list(20:33) = [4,-1,-1,-1,-1,-1,-1,5,6,7,8,9,10,11]
list = chn_list(chn)
sensor_string = 'msg_1_seviri'
case ( 'SEVIRI-MSG09')
!---stw chn_list(20:33) = [1,-1,-1,-1,-1,-1,-1,2,3,4,5,6,7,8]
chn_list(20:33) = [4,-1,-1,-1,-1,-1,-1,5,6,7,8,9,10,11]
list = chn_list(chn)
sensor_string = 'msg_2_seviri'
case ( 'SEVIRI-MSG10')
!---stw chn_list(20:33) = [1,-1,-1,-1,-1,-1,-1,2,3,4,5,6,7,8]
chn_list(20:33) = [4,-1,-1,-1,-1,-1,-1,5,6,7,8,9,10,11]
list = chn_list(chn)
sensor_string = 'msg_3_seviri'
case ( 'SEVIRI-MSG11')
!---stw chn_list(20:33) = [1,-1,-1,-1,-1,-1,-1,2,3,4,5,6,7,8]
chn_list(20:33) = [4,-1,-1,-1,-1,-1,-1,5,6,7,8,9,10,11] list = chn_list(chn)
sensor_string = 'msg_4_seviri'
case('OLCI')
chn_list = [6,9,3,4,10,14,15,1,2,-1,-1,-1,-1,-1,7,-1,-1,-1,-1,7,-1,18,-1,-1,-1,12,-1,-1,19,-1,20,22]
list = chn
list = (chn-1)*5*3 + 1
sensor_string = 'sentinel3_1_olci'
case ('SLSTR2')
chn_list = [6,9,3,4,10,14,15,1,2,-1,-1,-1,-1,-1,7,-1,-1,-1,-1,7,-1,18,-1,-1,-1,12,-1,-1,19,-1,20,22]
list = chn
sensor_string = 'sentinel3_2_slstr'
case('FY3-D')
Chn_List(20:32) = [ 20,-1,-1,21,-1,-1,-1,-1, 22,23,-1,24,25 ]
list = chn_list(chn)
sensor_string = 'fy3_4_mersi2'
case default
if (index(sensor,'AVHRR') .gt. 0) then
chn_list(1) = 1
chn_list(2) = 2
select case(sensor(11:12))
! -- sensor name can be AVHRR-NOAA??
case ('15','16','17','18','19')
chn_list(6) = 3
chn_list(20) = 4
chn_list(31) = 5
chn_list(32) = 6
case default
chn_list(20) = 3
chn_list(31) = 4
chn_list(32) = 5
end select
if (index(sensor,'AVHRR-METOP') .gt. 0) then
chn_list(6) = 3
chn_list(20) = 4
chn_list(31) = 5
chn_list(32) = 6
end if
sensor_string = 'noaa_'//sensor(11:12)//'_avhrr'
if ( index(sensor, 'NOAA0' ) .gt. 0 ) then
sensor_string = 'noaa_'//sensor(12:12)//'_avhrr'
rttov_version_string = '8'
chn_list(20) = 1
chn_list(31) = 2
chn_list(32) = 3
end if
list = chn_list(chn)
if (index(sensor,'AVHRR-METOP') .gt. 0 ) then
sensor_string = 'metop_'//metop_nr//'_avhrr'
end if
if (index(sensor,'AVHRR-TIROSN') .gt. 0) then
sensor_string = 'noaa_5_avhrr'
rttov_version_string = '8'
chn_list(20) = 1
chn_list(31) = 2
chn_list(32) = 3
end if
else if (index(sensor,'HIRS') .gt. 0) then
! 20 23 24 25 27 28 30 31 32 33 34 35 36
chn_List(20:36) = [ 19,-1,-1,18,15,14,-1,12,11,-1,9, 8,10, 7, 6, 5, 4 ]
list = chn_list(chn)
sensor_string = 'noaa_'//sensor(10:11)//'_hirs-shifted'
rttov_version_string = '8'
if (index(sensor,'AVHRR-METOP') .gt. 0 ) then
sensor_string = 'metop_'//metop_nr//'_hirs-shifted'
end if
if (index(sensor,'AVHRR-TIROSN') .gt. 0) then
sensor_string = 'noaa_5_hirs-shifted'
end if
else if (index(sensor,'GOES-') .gt. 0) then
rttov_version_string = '8'
chn_list(1) = 1
chn_list(20) = 2
chn_list(27) = 3
chn_list(31) = 4
chn_list(32) = 5
chn_list(33) = 6
sensor_string = 'goes_'//sensor(6:7)//'_imager'
if ( sensor(6:6) .ne. '1' ) then
sensor_string = 'goes_'//sensor(6:6)//'_imager'
rttov_version_string = '7'
end if
list = chn_list(chn)
else
print*,trim(sensor)
stop ' missing sensor in channel mapping of RTTOV inform andi.walther@ssec.wisc.edu '
end if
end select
coef_filename = trim(path)//'/rtcoef_rttov12/rttov'//rttov_version_string//'pred54L/rtcoef_'//trim(sensor_string)//'.dat'
cld_coef_filename = trim(path)//'cldaer_ir/sccldcoef_'//trim(sensor_string)//'.dat'
if ( list .eq. 0) list = -1
!print*,'Sensor rttov mapping: ',sensor,chn,list
deallocate(chn_list)
end function channel_map
subroutine compute_transmission_rttov ( &
ancil_data_path &
& , pstd &
& ,temp &
& ,wvmr &
& ,ozmr &
& ,theta &
& ,sensor &
& ,kban_in &
& ,taut &
& , use_modis_channel_equivalent )
implicit none
character(len = * ) , intent(in) :: ancil_data_path
real, intent(in) :: pstd (:,:)
real, intent(in) :: temp (:,:)
real, intent(in) :: wvmr(:,:)
real, intent(in) :: ozmr(:,:) real, intent(in) :: theta (:)
character (len =* ), intent(in) :: sensor
integer, intent(in) :: kban_in
logical , optional , intent(in) :: use_modis_channel_equivalent
real, intent(out) :: taut (:,:)
integer :: lll
real, parameter :: Q_MIXRATIO_TO_PPMV = 1.60771704e+6
opts % rt_ir % addsolar = .FALSE.
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 % addclouds = .FALSE. ! Don't include cloud effects
opts % rt_ir % addaerosl = .FALSE. ! Don't include aerosol effects
opts % rt_ir % ozone_data = .TRUE. ! Set the relevant flag to .TRUE.
opts % rt_ir % co2_data = .FALSE. ! when supplying a profile of the
opts % rt_ir % n2o_data = .FALSE. ! given trace gas (ensure the
opts % rt_ir % ch4_data = .FALSE. ! coef file supports the gas)
opts % rt_ir % co_data = .FALSE. !
opts % rt_ir % so2_data = .FALSE. !
opts % rt_mw % clw_data = .FALSE. !
opts % rt_mw % clw_scheme = 0
opts % config % verbose = .FALSE. ! Enable printing of warnings
!lll = channel_map (sensor,kban_in , coef_filename)
lll = channel_map (sensor,ancil_data_path, kban_in , coef_filename,cld_coef_filename)
nprof = size(temp(1,:))
nlevels = size(temp(:,1))
nthreads = 1
nchannels = 1
allocate(channel_list(nchannels))
if ( lll .lt. 0 ) then
if ( allocated(channel_list)) deallocate(channel_list)
return
end if
channel_list = [lll]
! --------------------------------------------------------------------------
! 2. Read coefficients
! --------------------------------------------------------------------------
CALL rttov_read_coefs(errorstatus, coefs, opts, file_coef=coef_filename)
IF (errorstatus /= errorstatus_success) THEN
WRITE(*,*) 'fatal error reading coefficients'
CALL rttov_exit(errorstatus)
ENDIF
! Ensure input number of channels is not higher than number stored in coefficient file
IF (nchannels > coefs % coef % fmv_chn) THEN
nchannels = coefs % coef % fmv_chn
ENDIF
! Ensure the options and coefficients are consistent
CALL rttov_user_options_checkinput(errorstatus, opts, coefs)
IF (errorstatus /= errorstatus_success) THEN
WRITE(*,*) 'error in rttov options' CALL rttov_exit(errorstatus)
ENDIF
! --------------------------------------------------------------------------
! 3. Allocate RTTOV input and output structures
! --------------------------------------------------------------------------
! Determine the total number of radiances to simulate (nchanprof).
! In this example we simulate all specified channels for each profile, but
! in general one can simulate a different number of channels for each profile.
nchanprof = nchannels * nprof
! Allocate structures for rttov_direct
CALL rttov_alloc_direct( &
errorstatus, &
1_jpim, & ! 1 => allocate
nprof, &
nchanprof, &
nlevels, &
chanprof, &
opts, &
profiles, &
coefs, &
transmission, &
radiance, &
calcemis=calcemis, &
emissivity=emissivity, &
calcrefl=calcrefl, &
reflectance=reflectance, &
init=.TRUE._jplm)
IF (errorstatus /= errorstatus_success) THEN
WRITE(*,*) 'allocation error for rttov_direct structures'
CALL rttov_exit(errorstatus)
ENDIF
! --------------------------------------------------------------------------
! 4. Build the list of profile/channel indices in chanprof
! --------------------------------------------------------------------------
nch = 0_jpim
DO j = 1, nprof
DO jch = 1, nchannels
nch = nch + 1_jpim
chanprof(nch)%prof = j
chanprof(nch)%chan = channel_list(jch)
ENDDO
ENDDO
deallocate ( channel_list)
! --------------------------------------------------------------------------
! 5. Read profile data
! --------------------------------------------------------------------------
profiles(:) % gas_units = 2 !
! Gas units (must be same for all profiles)
! 0 => ppmv over dry air
! 1 => kg/kg over moist air
! 2 => ppmv over moist air
!
! Loop over all profiles and read data for each one
DO iprof = 1, nprof
profiles(iprof) % t(:) = temp(:,iprof)
profiles(iprof) % p(:) = pstd(:,iprof)
profiles(iprof) % o3(:) = max(ozmr(:,iprof),0.1001E-10)
profiles(iprof) % q(:) = wvmr(:, iprof) * Q_MIXRATIO_TO_PPMV /1000.
profiles(iprof) % s2m % p =pstd(96,iprof)
profiles(iprof) % s2m % t = temp(96,iprof)
profiles(iprof) % s2m % q = wvmr(96,iprof) * Q_MIXRATIO_TO_PPMV /1000.
profiles(iprof) % s2m % u = 0.
profiles(iprof) % s2m % v = 0.
profiles(iprof) % s2m % wfetc = 1000000.
profiles(iprof) % skin % t = temp(96,iprof) + 5.
profiles(iprof) % zenangle = theta(iprof)
profiles(iprof) % azangle = 0.
profiles(iprof) % latitude = 45.
profiles(iprof) % longitude = 19.
! profiles(iprof) % zenangle = theta
! profiles(iprof) % zenangle = theta
ENDDO
CLOSE(iup)
! --------------------------------------------------------------------------
! 6. Specify surface emissivity and reflectance
! --------------------------------------------------------------------------
! In this example we have no values for input emissivities
emissivity(:) % emis_in = 0._jprb
! Calculate emissivity within RTTOV where the input emissivity value is
! zero or less (all channels in this case)
calcemis(:) = (emissivity(:) % emis_in <= 0._jprb)
! In this example we have no values for input reflectances
reflectance(:) % refl_in = 0._jprb
! Calculate BRDF within RTTOV where the input BRDF value is zero or less
! (all channels in this case)
calcrefl(:) = (reflectance(:) % refl_in <= 0._jprb)
! Use default cloud top BRDF for simple cloud in VIS/NIR channels
reflectance(:) % refl_cloud_top = 0._jprb
! --------------------------------------------------------------------------
! 7. Call RTTOV forward model
! --------------------------------------------------------------------------
CALL rttov_direct( &
errorstatus, &! out error flag
chanprof, &! in channel and profile index structure
opts, &! in options structure
profiles, &! in profile array
coefs, &! in coefficients structure
transmission, &! inout computed transmittances
radiance, &! inout computed radiances
calcemis = calcemis, &! in flag for internal emissivity calcs
emissivity = emissivity, &! inout input/output emissivities per channel
calcrefl = calcrefl, &! in flag for internal BRDF calcs
reflectance = reflectance) ! inout input/output BRDFs per channel
IF (errorstatus /= errorstatus_success) THEN
WRITE (*,*) 'rttov_direct error'
CALL rttov_exit(errorstatus)
ENDIF
taut = transmission % tau_levels
! Deallocate structures for rttov_direct
CALL rttov_alloc_direct( &
errorstatus, &
0_jpim, & ! 0 => deallocate
nprof, &
nchanprof, &
nlevels, &
chanprof, &
opts, &
profiles, &
coefs, &
transmission, &
radiance, &
calcemis=calcemis, &
emissivity=emissivity, &
calcrefl=calcrefl, &
reflectance=reflectance)
IF (errorstatus /= errorstatus_success) THEN
WRITE(*,*) 'deallocation error for rttov_direct structures'
CALL rttov_exit(errorstatus)
ENDIF
CALL rttov_dealloc_coefs(errorstatus, coefs)
IF (errorstatus /= errorstatus_success) THEN
WRITE(*,*) 'coefs deallocation error'
ENDIF
end subroutine compute_transmission_rttov
end module cx_rttov_bridge_mod