module m_prad
!$$$ subprogram documentation block
! . . . .
! subprogram: module m_prad
! prgmmr: j guo <jguo@nasa.gov>
! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3
! date: 2015-08-28
!
! abstract: a modular version of prad_bias and related data
!
! program history log:
! 2015-08-28 j guo - created this module on top of prad_bias();
! . completed with data components from obsmod;
! . changed code where this module needs to be used;
! . added this document block;
! . for earlier history log, see the history log section
! inside ::prad_bias() below.
! 2015-09-03 j guo - removed all older interface names.
! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting
!
! input argument list: see Fortran 90 style document below
!
! output argument list: see Fortran 90 style document below
!
! attributes:
! language: Fortran 90 and/or above
! machine:
!
!$$$ end subprogram documentation block
! module interface:
use m_obsNode , only: obsNode
use m_radNode , only: radNode
use m_radNode , only: radNode_typecast
use m_radNode , only: radNode_nextcast
use m_obsLList, only: obsLList
use m_obsLList, only: obsLList_headNode
use m_obsLList, only: obsLList_reset
implicit none
! Data objects:
public:: radheadm
! interfaces:
public:: prad_create !, create_passive_obsmod_vars
public:: prad_destroy !, destroyobs_passive
interface prad_create ; module procedure create_passive_obsmod_vars; end interface
interface prad_destroy; module procedure destroyobs_passive; end interface
public:: prad_updatePredx !, prad_bias
interface prad_updatePredx ; module procedure prad_bias; end interface
! Synopsis:
! - []_create: allocated in gsimod::gsimain_initialize(), through
! ::create_passive_obsmod_vars();
!
! - externally built, node-by-node, in setuprad();
!
! - []_updatePredx: used to update radinfo::predx, in glbsoi(), through
! ::prad_bias();
!
! - []_destroy: deallocated in gsimod::gsimain_finalize(), through
! ::create_passive_obsmod_vars().
! def radheadm - radiance linked list head for monitored radiance data
type(obsLList),dimension(:),pointer :: radheadm => null()
type(radNode),target,save:: radNode_mold
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
character(len=*),parameter :: myname='prad_bias'
contains
! ----------------------------------------------------------------------
subroutine create_passive_obsmod_vars
!$$$ subprogram documentation block
! . . . .
! subprogram: create_passive_obsmod_vars
! prgmmr: zhu org: np23 date: 2010-05-12
!
! abstract: allocate arrays to hold observation related information
!
! program history log:
! 2010-05-12 zhu
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$ end documentation block
use gsi_4dvar, only: nobs_bins
implicit none
ALLOCATE(radheadm(nobs_bins))
call lreset_(radheadm(:))
return
end subroutine create_passive_obsmod_vars
subroutine lreset_(llists)
use kinds, only: i_kind
implicit none
type(obsLList),dimension(:),intent(inout) :: llists
integer(kind=i_kind):: ib
do ib=lbound(llists,1),ubound(llists,1)
call obsLList_reset(llists(ib),mold=radNode_mold)
enddo
end subroutine lreset_
! ----------------------------------------------------------------------
subroutine destroyobs_passive
!$$$ subprogram documentation block
! . . . .
! subprogram: destroyobs_passive
! prgmmr: zhu org: np23 date: 2010-05-12
!
! abstract: deallocate arrays that hold observation information for
! use in outer and inner loops
!
! program history log:
! 2010-05-12 zhu
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$ end documentation block
implicit none
call lreset_(radheadm(:))
deallocate(radheadm)
return
end subroutine destroyobs_passive
! ----------------------------------------------------------------------
subroutine prad_bias()
!$$$ subprogram documentation block
! . . . .
! subprogram: prad_bias compute bias coefficients for passive radiances
! prgmmr: zhu org: np23 date: 2010-05-15
!
! abstract: compute bias coefficients for passive radiances
!
! program history log:
! 2010-05-15 zhu
! 2013-10-27 todling - move destroy to module where init reside (gsimod)
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind,r_quad
use mpimod, only: mype
use radinfo, only: npred,jpch_rad,iuse_rad,inew_rad ! intent(in)
use radinfo, only: predx ! intent(inout)
use gsi_4dvar, only: nobs_bins
use berror, only: varprd ! intent(in)
use mpl_allreducemod, only: mpl_allreduce
use timermod, only: timer_ini,timer_fnl
use constants, only : zero,one,zero_quad
implicit none
integer(i_kind),parameter :: nthreshold=100
real(r_kind),parameter :: atiny=1.0e-10_r_kind
integer(i_kind) i,n,j,ii,jj,jpassive,ibin,ic,mp,mm,kpred
integer(i_kind),dimension(jpch_rad) :: icp
integer(i_kind),dimension(npred) :: iorder
real(r_kind) varc
real(r_quad),dimension(npred,npred) :: Atmp
real(r_kind),allocatable,dimension(:) :: iobs
real(r_quad),allocatable,dimension(:,:,:) :: A
real(r_quad),allocatable,dimension(:,:) :: b
real(r_kind),allocatable,dimension(:,:) :: AA
real(r_kind),allocatable,dimension(:) :: be
class(obsNode),pointer:: obsptrm
type(radNode),pointer:: radptrm
! Initialize timer
call timer_ini('prad_bias')
! Allocate arrays
icp = 0
jpassive = 0
do j=1,jpch_rad
if (iuse_rad(j)==-1 .or. iuse_rad(j)==0) then
jpassive=jpassive+1
icp(j) = jpassive
end if
end do
if (mype==0) write(6,*) 'prad_bias: number of passive channels=', jpassive
if (jpassive<1) return
! Allocate arrays and initialize
allocate(A(npred,npred,jpassive),b(npred,jpassive))
allocate(iobs(jpassive))
do n=1,jpassive
iobs(n)=zero
do j=1,npred
b(j,n)=zero_quad
do i=1,npred
A(i,j,n)=zero_quad
end do
end do
end do
! Big loop for observations
do ibin=1,nobs_bins
!radptrm => obsLList_headNode(radheadm(ibin))
obsptrm => obsLList_headNode(radheadm(ibin))
radptrm => radNode_typecast(obsptrm)
do while (associated(radptrm))
if (radptrm%luse) then
! begin channel specific calculations
do n=1,radptrm%nchan
ic=radptrm%icx(n)
if (inew_rad(ic) .and. all(predx(:,ic)==zero)) cycle
mp=icp(ic)
! Assign arrays Ax=b
varc=radptrm%err2(n)*radptrm%raterr2(n)
if (sqrt(varc)<atiny) cycle
iobs(mp)=iobs(mp)+one
do i=1,npred
do j=1,npred
A(i,j,mp)=A(i,j,mp)+radptrm%pred(i,n)*radptrm%pred(j,n)*varc
end do
end do
do i=1,npred
b(i,mp)=b(i,mp)+radptrm%pred(i,n)*radptrm%res(n)*varc
end do
end do ! <n channel loop>
end if ! <luse>
! radptrm => radptrm%llpoint
radptrm => radNode_nextcast(radptrm)
end do
end do ! <ibin loop>
! Collect data from all processors
call mpl_allreduce(jpassive,rpvals=iobs)
call mpl_allreduce(npred,jpassive,b)
do n = 1,jpassive
if (iobs(n)<nthreshold) cycle
do i=1,jpch_rad
if (icp(i)==n) then
mm=i
exit
end if
end do
! Collect data from all processors for each channel
Atmp(:,:) = A(:,:,n)
call mpl_allreduce(npred,npred,Atmp)
! Solve linear system
iorder=0
kpred=0
do i=1,npred
use: do j=1,npred
if(Atmp(i,j) /= zero_quad)then
kpred = kpred+1
iorder(kpred) = i
exit use
end if
end do use
end do
if (kpred==0) cycle
allocate(AA(kpred,kpred),be(kpred))
do i = 1,kpred
ii=iorder(i)
be(i) = b(ii,n)
do j = 1,kpred
jj=iorder(j)
AA(i,j) = Atmp(ii,jj)
end do
end do
do i = 1,kpred
ii=iorder(i)
jj=(mm-1)*npred+ii
AA(i,i) = AA(i,i)+one/varprd(jj)
end do
if (all(abs(AA)<atiny)) cycle
if (all(abs(be)<atiny)) cycle
call linmm(AA,be,kpred,1,kpred,kpred)
! Update bias coefficients for passive channels
do i=1,kpred
ii=iorder(i)
predx(ii,mm)=predx(ii,mm)+be(i)
end do
deallocate(AA,be)
end do ! end of jpassive loop
deallocate(A,b,iobs)
! Finalize timer
call timer_fnl('prad_bias')
! End of routine
return
end subroutine prad_bias
end module m_prad