Skip to content

Commit

Permalink
Changed a section of CADS code inside setuprad to a subdirectory.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wx20jjung
wx20jjung committed Nov 7, 2023
1 parent e4646fd commit 1d1715b
Show file tree
Hide file tree
Showing 2 changed files with 188 additions and 116 deletions.
2 changes: 1 addition & 1 deletion src/gsi/qcmod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2180,7 +2180,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs
real(r_kind),dimension(nsig,nchanl),intent(in ) :: ptau5,temp,wmix
real(r_kind),dimension(nsig), intent(in ) :: prsltmp,tvp
real(r_kind),dimension(nchanl), intent(inout) :: errf,varinv,varinv_use
real(r_kind), intent(in ) :: cluster_fraction(:)
real(r_kind),dimension(7), intent(in ) :: cluster_fraction
real(r_kind),dimension(2,7), intent(in ) :: cluster_bt
real(r_kind),dimension(2), intent(in ) :: chan_stdev, model_bt

Expand Down
302 changes: 187 additions & 115 deletions src/gsi/setuprad.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -443,18 +443,9 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,&
logical:: muse_ii

! variables added for CADS
integer(i_kind) :: itmp1_cads, itmp2_cads, nchanl_cads, maxinfo, dval_info, cads_info
integer(i_kind),allocatable,dimension(:) :: ich_cads
logical :: imager_spccoeff, imager_taucoeff
real(r_kind),allocatable,dimension(:) :: tsim_cads, emissivity_cads, chan_level_cads
real(r_kind),allocatable,dimension(:) :: ts_cads, emissivity_k_cads,data_s_cads
real(r_kind),allocatable,dimension(:,:) :: ptau5_cads, temp_cads, wmix_cads, jacobian_cads
real(r_kind),dimension(7) :: imager_cluster_fraction
real(r_kind),dimension(2,7) :: imager_cluster_bt
real(r_kind),dimension(2) :: imager_chan_stdev, imager_model_bt
character(len=80) :: spc_filename, tau_filename
character(len=20) :: isis_cads
character(len=10) :: obstype_cads
real(r_kind),dimension(7,nobs) :: imager_cluster_fraction
real(r_kind),dimension(2,7,nobs) :: imager_cluster_bt
real(r_kind),dimension(2,nobs) :: imager_chan_stdev, imager_model_bt

! Notations in use: for a single obs. or a single obs. type
! nchanl : a known channel count of a given type obs stream
Expand Down Expand Up @@ -616,108 +607,12 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,&
end if

call dtime_setup()
cads_info = 0
! If using CADS do some initial checks, setup arrayss and calculate imager BTs
! cads_info = 0
! If using CADS setup arrays and calculate imager BTs
if ((iasi_cads .and. iasi) .or. (cris_cads .and. cris)) then
cads_info = 23
itmp1_cads = len(trim(obstype))
itmp2_cads = len(trim(isis))

if ( iasi ) then
isis_cads = 'avhrr3'//isis(itmp1_cads+1:itmp2_cads)
obstype_cads = 'avhrr'
nchanl_cads = 3 !channels 3 - 5
elseif ( cris ) then
! isis_cads = 'viirs-m'//isis(itmp1+1:itmp2) When naming convention becomes standarized with CrIS
if ( isis == 'cris-fsr_npp' .or. isis == 'cris_npp' ) then
isis_cads = 'viirs-m_npp'
elseif ( isis == 'cris-fsr_n20' ) then
isis_cads = 'viirs-m_n20'
spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin'
inquire(file=trim(spc_filename), exist=imager_spccoeff)
if ( .not. imager_spccoeff ) isis_cads = 'viirs-m_j1'
elseif ( isis == 'cris-fsr_n21' ) then
isis_cads = 'viirs-m_n21'
spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin'
inquire(file=trim(spc_filename), exist=imager_spccoeff)
if ( .not. imager_spccoeff ) isis_cads = 'viirs-m_j2'
endif
obstype_cads = 'viirs-m'
nchanl_cads = 5 ! channels 12 - 16
endif

spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin'
inquire(file=trim(spc_filename), exist=imager_spccoeff)
tau_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.TauCoeff.bin'
inquire(file=trim(tau_filename), exist=imager_taucoeff)

! IF the RTM files exist allocage and setup various arrays for the RTM
if ( imager_spccoeff .and. imager_taucoeff) then

nchanl_cads = 0
do i=1,jpch_rad
if (trim(isis_cads) == nusis(i)) then
nchanl_cads = nchanl_cads +1
endif
end do

allocate( ich_cads(nchanl_cads) )
jc = 0
do i=1,jpch_rad
if (trim(isis_cads) == nusis(i)) then
jc = jc +1
ich_cads(jc) = i
endif
end do

call init_crtm(init_pass,-99,mype,nchanl_cads,nreal,isis_cads,obstype_cads,radmod)

! Initialize variables needed for the infrared cloud and aerosol detections
! softwre.
allocate(data_s_cads(nreal+nchanl_cads),tsim_cads(nchanl_cads),emissivity_cads(nchanl_cads), &
chan_level_cads(nchanl_cads),ptau5_cads(nsig,nchanl_cads),ts_cads(nchanl_cads),emissivity_k_cads(nchanl_cads), &
temp_cads(nsig,nchanl_cads),wmix_cads(nsig,nchanl_cads), jacobian_cads(nsigradjac,nchanl_cads))

do n = 1,nobs ! loop to derive imager BTs for CADS
! Extract analysis relative observation time.
dtime = data_s(itime,n)
call dtime_check(dtime, in_curbin, in_anybin)
if(.not.in_anybin) cycle

maxinfo = nreal - cads_info - dval_info - nstinfo
if ( sum(data_s(maxinfo+1:maxinfo+7,n)) > 0.90_r_kind ) then ! imager cluster information exists for this profile
data_s_cads = data_s(1:nreal+nchanl_cads,n)
call call_crtm(obstype_cads,dtime,data_s_cads,nchanl_cads,nreal,ich_cads, &
tvp,qvp,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, &
trop5,tzbgr,dtsavg,sfc_speed,tsim_cads,emissivity_cads,chan_level_cads, &
ptau5_cads,ts_cads,emissivity_k_cads,temp_cads,wmix_cads,jacobian_cads,error_status)

! Transfer imager data to arrays for qc_irsnd
imager_cluster_fraction = data_s(maxinfo+1:maxinfo+7,n)
imager_cluster_bt(1,:) = data_s(maxinfo+8:maxinfo+14,n)
imager_cluster_bt(2,:) = data_s(maxinfo+15:maxinfo+21,n)
imager_chan_stdev = data_s(maxinfo+22:maxinfo+23,n)
imager_model_bt(1:2) = tsim_cads(nchanl_cads-1:nchanl_cads)
else
! if the imager cluster information does not exist, set arrays to zero
imager_cluster_fraction(:) = zero
imager_cluster_bt(:,:) = zero
imager_chan_stdev(:) = zero
imager_model_bt(:) = zero
endif ! imager information exists
end do ! End loop to derive imager BTs

call destroy_crtm
deallocate(data_s_cads,tsim_cads,emissivity_cads, ich_cads,chan_level_cads,ptau5_cads,&
ts_cads,emissivity_k_cads, temp_cads,wmix_cads, jacobian_cads)
else

! if the RTM files do not exist, set arrays to zero so CADS will ignore the imager tests.
imager_cluster_fraction(:) = zero
imager_cluster_bt(:,:) = zero
imager_chan_stdev(:) = zero
imager_model_bt(:) = zero
endif ! RTM files exist
call cads_imager_calc(obstype,isis,nobs,nreal,nchanl,nsig,data_s,init_pass,mype,radmod, &
imager_cluster_fraction,imager_cluster_bt,imager_chan_stdev, imager_model_bt)
endif ! using cads

if ( mype == 0 .and. .not.l_may_be_passive) write(6,*)mype,'setuprad: passive obs',is,isis
Expand Down Expand Up @@ -980,13 +875,11 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,&

! Set relative weight value
val_obs=one
dval_info = 0
if(dval_use)then
ixx=nint(data_s(nreal-nstinfo,n))
if (ixx > 0 .and. super_val1(ixx) >= one) then
val_obs=data_s(nreal-nstinfo-1,n)/super_val1(ixx)
endif
dval_info = 2
endif

do jc=1,nchanl
Expand Down Expand Up @@ -1500,7 +1393,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,&
hirs,zsges,cenlat,cenlon,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, &
wavenumber,ptau5,prsltmp,tvp,temp,wmix,chan_level,emissivity_k,ts,tsim, &
id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole(n), &
imager_cluster_fraction, imager_cluster_bt, imager_chan_stdev,imager_model_bt )
imager_cluster_fraction(7,n), imager_cluster_bt(2,7,n), imager_chan_stdev(2,n),imager_model_bt(2,n))

! --------- MSU -------------------
! QC MSU data
Expand Down Expand Up @@ -2870,4 +2763,183 @@ subroutine final_binary_diag_
close(4)
end subroutine final_binary_diag_
end subroutine setuprad

subroutine cads_imager_calc(obstype,isis,nobs,nreal,nchanl,nsig,data_s,init_pass,mype,radmod, &
imager_cluster_fraction,imager_cluster_bt,imager_chan_stdev, imager_model_bt)

!$$$ subprogram documentation block
!
! subprogram: cads_imager_calc compute model equivalent to the imager channels used by CADS
! prgmmr: Jung
!
! abstract: accumulate the data necessary to derive the model equivalent brightness temperatures
! used by the cloud and aerosol detection software for the imager cloud tests.
!
! program history log:
!
!
!
! subroutines included:
!
!
! input argument list:
!
! obstype - type of tb observation
! isis - sensor/instrument/satellite id
! nobs - number of observations
! nreal - number of pieces of info (location, time, etc) per obs
! nchanl - number of channels per obs
! nsig - number of model layers
! data_s - array containing input data information for a specific sensor
! init_pass - state of "setup" processing
! mype - mpi task id
! radmod - structure for cloud and aerosol information
!
! output argument list:

! imager_cluster_fraction - CADS cluster fraction ( dimension 7)
! imager_cluster_bt - avreage brightness temperature of a cluster
! imager_chan stdev - brightness temperature standard deviation of the cluster
! imager_model_bt - model derived brightness temperature
!
!
!$$$ end documentation block

use kinds, only: i_kind, r_kind
use constants, only: zero
use radiance_mod, only: rad_obs_type
use radinfo, only: jpch_rad, nusis, crtm_coeffs_path, nsigradjac
use crtm_interface, only: init_crtm, call_crtm, destroy_crtm, itime
use obsmod, only: dval_use
use gsi_nstcouplermod, only: nstinfo

implicit none

logical, intent(in) :: init_pass
character(len=10), intent(in) :: obstype
character(len=20), intent(in) :: isis
integer(i_kind), intent(in) :: nobs, nreal, nchanl, nsig
integer(i_kind), intent(in) :: mype
type(rad_obs_type), intent(in) :: radmod
real(r_kind),dimension(nreal+nchanl,nobs),intent(in) :: data_s
real(r_kind),dimension(7,nobs), intent(out) :: imager_cluster_fraction
real(r_kind),dimension(2,7,nobs), intent(out) :: imager_cluster_bt
real(r_kind),dimension(2,nobs), intent(out) :: imager_chan_stdev, imager_model_bt

! local variables
integer(i_kind) :: jc, i, n
integer(i_kind) :: itmp1_cads, itmp2_cads, nchanl_cads, maxinfo, dval_info, cads_info, error_status
integer(i_kind),allocatable,dimension(:) :: ich_cads
logical :: imager_spccoeff, imager_taucoeff
real(r_kind) :: dtime, clw_guess, ciw_guess, rain_guess, snow_guess
real(r_kind) :: trop5, tzbgr, dtsavg, sfc_speed
real(r_kind),dimension(nsig) :: qvp, tvp, qs, prsltmp
real(r_kind),dimension(nsig+1) :: prsitmp
real(r_kind),allocatable,dimension(:) :: tsim_cads, emissivity_cads, chan_level_cads
real(r_kind),allocatable,dimension(:) :: ts_cads, emissivity_k_cads,data_s_cads
real(r_kind),allocatable,dimension(:,:) :: ptau5_cads, temp_cads, wmix_cads, jacobian_cads
character(len=80) :: spc_filename, tau_filename
character(len=20) :: isis_cads
character(len=10) :: obstype_cads

cads_info = 23
dval_info = 0
if (dval_use) dval_info = 2

itmp1_cads = len(trim(obstype))
itmp2_cads = len(trim(isis))

if ( obstype == 'iasi' ) then
isis_cads = 'avhrr3'//isis(itmp1_cads+1:itmp2_cads)
obstype_cads = 'avhrr'
! nchanl_cads = 3 !channels 3 - 5
elseif ( obstype == 'cris' .or. obstype == 'cris-fsr' ) then
! isis_cads = 'viirs-m'//isis(itmp1+1:itmp2) When naming convention becomes standarized with CrIS
if ( isis == 'cris-fsr_npp' .or. isis == 'cris_npp' ) then
isis_cads = 'viirs-m_npp'
elseif ( isis == 'cris-fsr_n20' ) then
isis_cads = 'viirs-m_n20'
spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin'
inquire(file=trim(spc_filename), exist=imager_spccoeff)
if ( .not. imager_spccoeff ) isis_cads = 'viirs-m_j1'
elseif ( isis == 'cris-fsr_n21' ) then
isis_cads = 'viirs-m_n21'
spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin'
inquire(file=trim(spc_filename), exist=imager_spccoeff)
if ( .not. imager_spccoeff ) isis_cads = 'viirs-m_j2'
endif
obstype_cads = 'viirs-m'
! nchanl_cads = 5 ! channels 12 - 16
endif

spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin'
inquire(file=trim(spc_filename), exist=imager_spccoeff)
tau_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.TauCoeff.bin'
inquire(file=trim(tau_filename), exist=imager_taucoeff)

! IF the RTM files exist allocate and setup various arrays for the RTM
if ( imager_spccoeff .and. imager_taucoeff) then
nchanl_cads = 0
do i=1,jpch_rad
if (trim(isis_cads) == nusis(i)) then
nchanl_cads = nchanl_cads +1
endif
end do

allocate( ich_cads(nchanl_cads) )
jc = 0
do i=1,jpch_rad
if (trim(isis_cads) == nusis(i)) then
jc = jc +1
ich_cads(jc) = i
endif
end do

call init_crtm(init_pass,-99,mype,nchanl_cads,nreal,isis_cads,obstype_cads,radmod)

! Initialize variables needed for the infrared cloud and aerosol detection software
allocate(data_s_cads(nreal+nchanl_cads),tsim_cads(nchanl_cads),emissivity_cads(nchanl_cads), &
chan_level_cads(nchanl_cads),ptau5_cads(nsig,nchanl_cads),ts_cads(nchanl_cads),emissivity_k_cads(nchanl_cads), &
temp_cads(nsig,nchanl_cads),wmix_cads(nsig,nchanl_cads), jacobian_cads(nsigradjac,nchanl_cads))

do n = 1,nobs ! loop to derive imager BTs for CADS
! Extract analysis relative observation time.
dtime = data_s(itime,n)
maxinfo = nreal - cads_info - dval_info - nstinfo
if ( sum(data_s(maxinfo+1:maxinfo+7,n)) > 0.90_r_kind ) then ! imager cluster information exists for this profile
data_s_cads = data_s(1:nreal+nchanl_cads,n)
call call_crtm(obstype_cads,dtime,data_s_cads,nchanl_cads,nreal,ich_cads, &
tvp,qvp,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, &
trop5,tzbgr,dtsavg,sfc_speed,tsim_cads,emissivity_cads,chan_level_cads, &
ptau5_cads,ts_cads,emissivity_k_cads,temp_cads,wmix_cads,jacobian_cads,error_status)

! Transfer imager data to arrays for qc_irsnd
imager_cluster_fraction(1:7,n) = data_s(maxinfo+1:maxinfo+7,n)
imager_cluster_bt(1,1:7,n) = data_s(maxinfo+8:maxinfo+14,n)
imager_cluster_bt(2,1:7,n) = data_s(maxinfo+15:maxinfo+21,n)
imager_chan_stdev(1:2,n) = data_s(maxinfo+22:maxinfo+23,n)
imager_model_bt(1:2,n) = tsim_cads(nchanl_cads-1:nchanl_cads)
else
! if the imager cluster information does not exist, set arrays to zero
imager_cluster_fraction(1:7,n) = zero
imager_cluster_bt(1:2,1:7,n) = zero
imager_chan_stdev(1:2,n) = zero
imager_model_bt(1:2,n) = zero
endif ! imager information exists
end do ! End loop to derive imager BTs

call destroy_crtm
deallocate(data_s_cads,tsim_cads,emissivity_cads, ich_cads,chan_level_cads,ptau5_cads,&
ts_cads,emissivity_k_cads, temp_cads,wmix_cads, jacobian_cads)

! if the RTM files do not exist, set arrays to zero so CADS will ignore the imager tests.
else
imager_cluster_fraction(1:7,1:nobs) = zero
imager_cluster_bt(1:2,1:7,1:nobs) = zero
imager_chan_stdev(1:2,1:nobs) = zero
imager_model_bt(1:2,1:nobs) = zero
endif ! RTM files exist

end subroutine cads_imager_calc

end module rad_setup

0 comments on commit 1d1715b

Please sign in to comment.