From e9aa185e579ca5adf1eea2d334b927c7af0085f2 Mon Sep 17 00:00:00 2001 From: "michael.lueken" Date: Fri, 10 Jul 2020 17:58:15 +0000 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#13. Continuing to clear through coding standard issues in the master. Finished through src/gsi/m_find.f90. --- src/gsi/m_colvkNode.F90 | 366 ----- src/gsi/m_colvknode.F90 | 366 +++++ ...{m_cvgridLookup.F90 => m_cvgridlookup.F90} | 290 ++-- src/gsi/m_dbzNode.F90 | 264 --- src/gsi/m_dbznode.F90 | 264 +++ src/gsi/m_dgeevx.F90 | 21 +- src/gsi/m_distance.f90 | 12 +- src/gsi/m_dtime.F90 | 22 +- src/gsi/m_dwNode.F90 | 256 --- src/gsi/m_dwnode.F90 | 256 +++ src/gsi/m_extOzone.F90 | 1452 ----------------- src/gsi/m_extozone.F90 | 1452 +++++++++++++++++ src/gsi/m_find.f90 | 20 +- 13 files changed, 2521 insertions(+), 2520 deletions(-) delete mode 100644 src/gsi/m_colvkNode.F90 create mode 100644 src/gsi/m_colvknode.F90 rename src/gsi/{m_cvgridLookup.F90 => m_cvgridlookup.F90} (60%) delete mode 100644 src/gsi/m_dbzNode.F90 create mode 100644 src/gsi/m_dbznode.F90 delete mode 100644 src/gsi/m_dwNode.F90 create mode 100644 src/gsi/m_dwnode.F90 delete mode 100644 src/gsi/m_extOzone.F90 create mode 100644 src/gsi/m_extozone.F90 diff --git a/src/gsi/m_colvkNode.F90 b/src/gsi/m_colvkNode.F90 deleted file mode 100644 index 3db8c2f73e..0000000000 --- a/src/gsi/m_colvkNode.F90 +++ /dev/null @@ -1,366 +0,0 @@ -module m_colvkNode -!$$$ subprogram documentation block -! . . . . -! subprogram: module m_colvkNode -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 -! date: 2016-05-18 -! -! abstract: class-module of obs-type colvkNode (sbuv co) -! -! program history log: -! 2016-05-18 j guo - added this document block for the initial polymorphic -! implementation. -! -! 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_obsdiagNode, only: obs_diag,aofp_obs_diag => fptr_obsdiagNode - use m_obsdiagNode, only: obs_diags - use kinds , only: i_kind,r_kind - use mpeu_util, only: assert_,die,perr,warn,tell - use m_obsNode, only: obsNode - implicit none - private - - public:: colvkNode - - type,extends(obsNode):: colvkNode - !type(colvk_ob_type),pointer :: llpoint => NULL() - type(aofp_obs_diag), dimension(:), pointer :: diags => NULL() - real(r_kind),dimension(:),pointer :: res => NULL() - ! co residual - real(r_kind),dimension(:),pointer :: err2 => NULL() - ! co error squared - real(r_kind),dimension(:),pointer :: raterr2 => NULL() - ! square of ratio of final obs error - ! to original obs error - !real(r_kind) :: time ! observation time in sec - real(r_kind),dimension( :),pointer :: wk => NULL() - real(r_kind),dimension(:,:),pointer :: wij => NULL() - ! horizontal interpolation weights - real(r_kind),dimension(:),pointer :: prs => NULL() - ! pressure levels - integer(i_kind),dimension(:),pointer :: ipos => NULL() - real(r_kind),dimension(:,:),pointer :: ak => NULL() - ! MOPITT vertical averaging kernel - real(r_kind),dimension(:),pointer :: ap => NULL() - ! MOPITT a priori - !real(r_kind),dimension(:),pointer :: wkk1 => NULL() - !real(r_kind),dimension(:),pointer :: wkk2 => NULL() - ! vertical intropolation weights for MOPITT - - integer(i_kind) :: nlco =0 ! number of levels for this profile - integer(i_kind) :: ij(4)=0 ! horizontal locations - !logical :: luse ! flag indicating if ob is used in pen. - - !integer(i_kind) :: idv,iob ! device id and obs index for sorting - !real (r_kind) :: elat, elon ! earth lat-lon for redistribution - !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution - contains - procedure,nopass:: mytype - procedure:: setHop => obsNode_setHop_ - procedure:: xread => obsNode_xread_ - procedure:: xwrite => obsNode_xwrite_ - procedure:: isvalid => obsNode_isvalid_ - procedure:: getTLDdp => getTLDdp_ - - ! procedure, nopass:: headerRead => obsHeader_read_ - ! procedure, nopass:: headerWrite => obsHeader_write_ - ! procedure:: init => obsNode_init_ - procedure:: clean => obsNode_clean_ - - end type colvkNode - - public:: colvkNode_typecast - public:: colvkNode_nextcast - interface colvkNode_typecast; module procedure typecast_ ; end interface - interface colvkNode_nextcast; module procedure nextcast_ ; end interface - - public:: colvkNode_appendto - interface colvkNode_appendto; module procedure appendto_ ; end interface - - character(len=*),parameter:: MYNAME="m_colvkNode" - -#include "myassert.H" -#include "mytrace.H" -contains -function typecast_(aNode) result(ptr_) -!-- cast a class(obsNode) to a type(colvkNode) - use m_obsNode, only: obsNode - implicit none - type(colvkNode),pointer:: ptr_ - class(obsNode ),pointer,intent(in):: aNode - ptr_ => null() - if(.not.associated(aNode)) return - ! logically, typecast of a null-reference is a null pointer. - select type(aNode) - type is(colvkNode) - ptr_ => aNode - end select -return -end function typecast_ - -function nextcast_(aNode) result(ptr_) -!-- cast an obsNode_next(obsNode) to a type(colvkNode) - use m_obsNode, only: obsNode,obsNode_next - implicit none - type(colvkNode),pointer:: ptr_ - class(obsNode ),target ,intent(in):: aNode - - class(obsNode),pointer:: inode_ - inode_ => obsNode_next(aNode) - ptr_ => typecast_(inode_) -return -end function nextcast_ - -subroutine appendto_(aNode,oll) - use m_obsNode , only: obsNode - use m_obsLList, only: obsLList,obsLList_appendNode - implicit none - type(colvkNode),pointer,intent(in):: aNode - type(obsLList),intent(inout):: oLL - - class(obsNode),pointer:: inode_ - inode_ => aNode - call obsLList_appendNode(oLL,inode_) - inode_ => null() -end subroutine appendto_ - -! obsNode implementations - -function mytype() - implicit none - character(len=:),allocatable:: mytype - mytype="[colvkNode]" -end function mytype - -subroutine obsNode_clean_(aNode) - implicit none - class(colvkNode),intent(inout):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_clean_' -_ENTRY_(myname_) -!_TRACEV_(myname_,'%mytype() =',aNode%mytype()) - if(associated(aNode%res )) deallocate(aNode%res ) - if(associated(aNode%diags )) deallocate(aNode%diags ) - if(associated(aNode%err2 )) deallocate(aNode%err2 ) - if(associated(aNode%raterr2)) deallocate(aNode%raterr2) - if(associated(aNode%prs )) deallocate(aNode%prs ) - if(associated(aNode%ipos )) deallocate(aNode%ipos ) - if(associated(aNode%ak )) deallocate(aNode%ak ) - if(associated(aNode%ap )) deallocate(aNode%ap ) - if(associated(aNode%wk )) deallocate(aNode%wk ) - if(associated(aNode%wij )) deallocate(aNode%wij ) - -_EXIT_(myname_) -return -end subroutine obsNode_clean_ - -subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) - use gridmod, only: nsig - use m_obsdiagNode, only: obsdiagLookup_locate - implicit none - class(colvkNode) , intent(inout):: aNode - integer(i_kind) , intent(in ):: iunit - integer(i_kind) , intent( out):: istat - type(obs_diags) , intent(in ):: diagLookup - logical,optional, intent(in ):: skip - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xread_' - integer(i_kind),allocatable,dimension(:):: ich - integer(i_kind):: k,nlco,nlevp - logical:: skip_ -_ENTRY_(myname_) - skip_=.false. - if(present(skip)) skip_=skip - - istat=0 - if(skip_) then - read(iunit,iostat=istat) - if (istat/=0) then - call perr(myname_,'skipping read(%nlco), istat =',istat) - _EXIT_(myname_) - return - end if - read(iunit,iostat=istat) - if(istat/=0) then - call perr(myname_,'skipping read(%(res,err2,...)), istat =',istat) - _EXIT_(myname_) - return - endif - - else - read(iunit,iostat=istat) aNode%nlco - if (istat/=0) then - call perr(myname_,'read(%nlco), istat =',istat) - _EXIT_(myname_) - return - end if - - call die(myname_,'unfinished implementation here') - ! The size of %wij(8,:) seems in a mess. It is nlev, nlevs, or nsig? - ! See setupco() for additional details. - - ! re-allocation is needed, because a node may have been read in - ! but excluded later. - if(associated(aNode%res )) deallocate(aNode%res ) - if(associated(aNode%diags )) deallocate(aNode%diags ) - if(associated(aNode%err2 )) deallocate(aNode%err2 ) - if(associated(aNode%raterr2)) deallocate(aNode%raterr2) - if(associated(aNode%prs )) deallocate(aNode%prs ) - if(associated(aNode%ipos )) deallocate(aNode%ipos ) - if(associated(aNode%ak )) deallocate(aNode%ak ) - if(associated(aNode%ap )) deallocate(aNode%ap ) - if(associated(aNode%wk )) deallocate(aNode%wk ) - if(associated(aNode%wij )) deallocate(aNode%wij ) - - nlco = aNode%nlco - nlevp= max(nlco,1) - - allocate( aNode%res (nlco) , & - aNode%diags (nlco) , & - aNode%err2 (nlco) , & - aNode%raterr2(nlco) , & - aNode%prs (nlevp), & - aNode%ipos (nlco) , & - aNode%ak(nlco,nlco) , & - aNode%ap (nlco) , & - aNode%wk (nsig) , & - aNode%wij (8,nsig) ) - - allocate( ich (nlco) ) - - read(iunit,iostat=istat) ich , & !(nlco) - aNode%res , & !(nlco) - aNode%err2 , & !(nlco) - aNode%raterr2, & !(nlco) - aNode%prs , & !(nlevp) - aNode%ipos , & !(nlco) - aNode%ak , & !(nlco,nlco) - aNode%ap , & !(nlco) - aNode%wk , & !( nsig) - aNode%wij , & !(8,nsig) - aNode%ij !(4) - if (istat/=0) then - call perr(myname_,'read(ich,%(...)), istat =',istat) - _EXIT_(myname_) - return - end if - - do k=1,nlco - ! Not sure if ich=k or ich=ich(k) for now. More verification is needed. - aNode%diags(k)%ptr => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,ich(k)) - if(.not.associated(aNode%diags(k)%ptr)) then - call perr(myname_,'obsdiagLookup_locate(k), k =',aNode%idv) - call perr(myname_,' %idv =',aNode%idv) - call perr(myname_,' %iob =',aNode%iob) - call perr(myname_,' ich =',ich(k)) - call die(myname_) - istat=-1 - _EXIT_(myname_) - return - endif - enddo - deallocate(ich) - endif -_EXIT_(myname_) -return -end subroutine obsNode_xread_ - -subroutine obsNode_xwrite_(aNode,junit,jstat) - implicit none - class(colvkNode),intent(in):: aNode - integer(i_kind),intent(in ):: junit - integer(i_kind),intent( out):: jstat - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xwrite_' - integer(i_kind):: k -_ENTRY_(myname_) - - jstat=0 - write(junit,iostat=jstat) aNode%nlco - if (jstat/=0) then - call perr(myname_,'write(%nlco), jstat =',jstat) - _EXIT_(myname_) - return - end if - - write(junit,iostat=jstat) (/ (aNode%diags(k)%ptr%ich, k=1,aNode%nlco) /), & !(nlco) - aNode%res , & !(nlco) - aNode%err2 , & !(nlco) - aNode%raterr2, & !(nlco) - aNode%prs , & !(nlevp) - aNode%ipos , & !(nlco) - aNode%ak , & !(nlco,nlco) - aNode%ap , & !(nlco) - aNode%wk , & !( nsig) - aNode%wij , & !(8,nsig) - aNode%ij !(4) - if (jstat/=0) then - call perr(myname_,'write(%(res,err2,...)), istat =',jstat) - _EXIT_(myname_) - return - end if -_EXIT_(myname_) -return -end subroutine obsNode_xwrite_ - -subroutine obsNode_setHop_(aNode) - !use m_obsNode, only: obstype_getHop => obsNode_getHop - use m_cvgridLookup, only: cvgridLookup_getiw - use gridmod, only: nsig - implicit none - class(colvkNode),intent(inout):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' - real(r_kind),dimension(4):: zwij - integer(i_kind):: k -_ENTRY_(myname_) - call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%ij,zwij) - do k=1,nsig - aNode%wij(1:4,k) = zwij(1:4)*aNode%wk(k) - aNode%wij(5:8,k) = zwij(1:4)*(1._r_kind-aNode%wk(k)) - enddo -_EXIT_(myname_) -return -end subroutine obsNode_setHop_ - -function obsNode_isvalid_(aNode) result(isvalid_) - implicit none - logical:: isvalid_ - class(colvkNode),intent(in):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_' - integer(i_kind):: k -_ENTRY_(myname_) - isvalid_ = all( (/ (associated(aNode%diags(k)%ptr), k=1,aNode%nlco) /) ) -_EXIT_(myname_) -return -end function obsNode_isvalid_ - -pure subroutine getTLDdp_(aNode,jiter,tlddp,nob) - use kinds, only: r_kind - implicit none - class(colvkNode), intent(in):: aNode - integer(kind=i_kind),intent(in):: jiter - real(kind=r_kind),intent(inout):: tlddp - integer(kind=i_kind),optional,intent(inout):: nob - - integer(kind=i_kind):: k - do k=1,aNode%nlco - tlddp = tlddp + aNode%diags(k)%ptr%tldepart(jiter)*aNode%diags(k)%ptr%tldepart(jiter) - enddo - if(present(nob)) nob=nob+aNode%nlco -return -end subroutine getTLDdp_ - -end module m_colvkNode diff --git a/src/gsi/m_colvknode.F90 b/src/gsi/m_colvknode.F90 new file mode 100644 index 0000000000..ca64b66769 --- /dev/null +++ b/src/gsi/m_colvknode.F90 @@ -0,0 +1,366 @@ +module m_colvknode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_colvknode +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2016-05-18 +! +! abstract: class-module of obs-type colvknode (sbuv co) +! +! program history log: +! 2016-05-18 j guo - added this document block for the initial polymorphic +! implementation. +! +! 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_obsdiagnode, only: obs_diag,aofp_obs_diag => fptr_obsdiagnode + use m_obsdiagnode, only: obs_diags + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsnode, only: obsnode + implicit none + private + + public:: colvknode + + type,extends(obsnode):: colvknode + !type(colvk_ob_type),pointer :: llpoint => null() + type(aofp_obs_diag), dimension(:), pointer :: diags => null() + real(r_kind),dimension(:),pointer :: res => null() + ! co residual + real(r_kind),dimension(:),pointer :: err2 => null() + ! co error squared + real(r_kind),dimension(:),pointer :: raterr2 => null() + ! square of ratio of final obs error + ! to original obs error + !real(r_kind) :: time ! observation time in sec + real(r_kind),dimension( :),pointer :: wk => null() + real(r_kind),dimension(:,:),pointer :: wij => null() + ! horizontal interpolation weights + real(r_kind),dimension(:),pointer :: prs => null() + ! pressure levels + integer(i_kind),dimension(:),pointer :: ipos => null() + real(r_kind),dimension(:,:),pointer :: ak => null() + ! mopitt vertical averaging kernel + real(r_kind),dimension(:),pointer :: ap => null() + ! mopitt a priori + !real(r_kind),dimension(:),pointer :: wkk1 => null() + !real(r_kind),dimension(:),pointer :: wkk2 => null() + ! vertical intropolation weights for mopitt + + integer(i_kind) :: nlco =0 ! number of levels for this profile + integer(i_kind) :: ij(4)=0 ! horizontal locations + !logical :: luse ! flag indicating if ob is used in pen. + + !integer(i_kind) :: idv,iob ! device id and obs index for sorting + !real (r_kind) :: elat, elon ! earth lat-lon for redistribution + !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution + contains + procedure,nopass:: mytype + procedure:: sethop => obsnode_sethop_ + procedure:: xread => obsnode_xread_ + procedure:: xwrite => obsnode_xwrite_ + procedure:: isvalid => obsnode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerread => obsheader_read_ + ! procedure, nopass:: headerwrite => obsheader_write_ + ! procedure:: init => obsnode_init_ + procedure:: clean => obsnode_clean_ + + end type colvknode + + public:: colvknode_typecast + public:: colvknode_nextcast + interface colvknode_typecast; module procedure typecast_ ; end interface + interface colvknode_nextcast; module procedure nextcast_ ; end interface + + public:: colvknode_appendto + interface colvknode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: myname="m_colvknode" + +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(anode) result(ptr_) +!-- cast a class(obsnode) to a type(colvknode) + use m_obsnode, only: obsnode + implicit none + type(colvknode),pointer:: ptr_ + class(obsnode ),pointer,intent(in):: anode + ptr_ => null() + if(.not.associated(anode)) return + ! logically, typecast of a null-reference is a null pointer. + select type(anode) + type is(colvknode) + ptr_ => anode + end select + return +end function typecast_ + +function nextcast_(anode) result(ptr_) +!-- cast an obsnode_next(obsnode) to a type(colvknode) + use m_obsnode, only: obsnode,obsnode_next + implicit none + type(colvknode),pointer:: ptr_ + class(obsnode ),target ,intent(in):: anode + + class(obsnode),pointer:: inode_ + inode_ => obsnode_next(anode) + ptr_ => typecast_(inode_) + return +end function nextcast_ + +subroutine appendto_(anode,oll) + use m_obsnode , only: obsnode + use m_obsllist, only: obsllist,obsllist_appendnode + implicit none + type(colvknode),pointer,intent(in):: anode + type(obsllist),intent(inout):: oll + + class(obsnode),pointer:: inode_ + inode_ => anode + call obsllist_appendnode(oll,inode_) + inode_ => null() +end subroutine appendto_ + +! obsnode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[colvknode]" +end function mytype + +subroutine obsnode_clean_(anode) + implicit none + class(colvknode),intent(inout):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_clean_' +_ENTRY_(myname_) +!_TRACEV_(myname_,'%mytype() =',anode%mytype()) + if(associated(anode%res )) deallocate(anode%res ) + if(associated(anode%diags )) deallocate(anode%diags ) + if(associated(anode%err2 )) deallocate(anode%err2 ) + if(associated(anode%raterr2)) deallocate(anode%raterr2) + if(associated(anode%prs )) deallocate(anode%prs ) + if(associated(anode%ipos )) deallocate(anode%ipos ) + if(associated(anode%ak )) deallocate(anode%ak ) + if(associated(anode%ap )) deallocate(anode%ap ) + if(associated(anode%wk )) deallocate(anode%wk ) + if(associated(anode%wij )) deallocate(anode%wij ) + +_EXIT_(myname_) + return +end subroutine obsnode_clean_ + +subroutine obsnode_xread_(anode,iunit,istat,diaglookup,skip) + use gridmod, only: nsig + use m_obsdiagnode, only: obsdiaglookup_locate + implicit none + class(colvknode) , intent(inout):: anode + integer(i_kind) , intent(in ):: iunit + integer(i_kind) , intent( out):: istat + type(obs_diags) , intent(in ):: diaglookup + logical,optional, intent(in ):: skip + + character(len=*),parameter:: myname_=myname//'::obsnode_xread_' + integer(i_kind),allocatable,dimension(:):: ich + integer(i_kind):: k,nlco,nlevp + logical:: skip_ +_ENTRY_(myname_) + skip_=.false. + if(present(skip)) skip_=skip + + istat=0 + if(skip_) then + read(iunit,iostat=istat) + if (istat/=0) then + call perr(myname_,'skipping read(%nlco), istat =',istat) + _EXIT_(myname_) + return + end if + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(res,err2,...)), istat =',istat) + _EXIT_(myname_) + return + endif + + else + read(iunit,iostat=istat) anode%nlco + if (istat/=0) then + call perr(myname_,'read(%nlco), istat =',istat) + _EXIT_(myname_) + return + end if + + call die(myname_,'unfinished implementation here') + ! The size of %wij(8,:) seems in a mess. It is nlev, nlevs, or nsig? + ! See setupco() for additional details. + + ! re-allocation is needed, because a node may have been read in + ! but excluded later. + if(associated(anode%res )) deallocate(anode%res ) + if(associated(anode%diags )) deallocate(anode%diags ) + if(associated(anode%err2 )) deallocate(anode%err2 ) + if(associated(anode%raterr2)) deallocate(anode%raterr2) + if(associated(anode%prs )) deallocate(anode%prs ) + if(associated(anode%ipos )) deallocate(anode%ipos ) + if(associated(anode%ak )) deallocate(anode%ak ) + if(associated(anode%ap )) deallocate(anode%ap ) + if(associated(anode%wk )) deallocate(anode%wk ) + if(associated(anode%wij )) deallocate(anode%wij ) + + nlco = anode%nlco + nlevp= max(nlco,1) + + allocate( anode%res (nlco) , & + anode%diags (nlco) , & + anode%err2 (nlco) , & + anode%raterr2(nlco) , & + anode%prs (nlevp), & + anode%ipos (nlco) , & + anode%ak(nlco,nlco) , & + anode%ap (nlco) , & + anode%wk (nsig) , & + anode%wij (8,nsig) ) + + allocate( ich (nlco) ) + + read(iunit,iostat=istat) ich , & !(nlco) + anode%res , & !(nlco) + anode%err2 , & !(nlco) + anode%raterr2, & !(nlco) + anode%prs , & !(nlevp) + anode%ipos , & !(nlco) + anode%ak , & !(nlco,nlco) + anode%ap , & !(nlco) + anode%wk , & !( nsig) + anode%wij , & !(8,nsig) + anode%ij !(4) + if (istat/=0) then + call perr(myname_,'read(ich,%(...)), istat =',istat) + _EXIT_(myname_) + return + end if + + do k=1,nlco + ! Not sure if ich=k or ich=ich(k) for now. More verification is needed. + anode%diags(k)%ptr => obsdiaglookup_locate(diaglookup,anode%idv,anode%iob,ich(k)) + if(.not.associated(anode%diags(k)%ptr)) then + call perr(myname_,'obsdiaglookup_locate(k), k =',anode%idv) + call perr(myname_,' %idv =',anode%idv) + call perr(myname_,' %iob =',anode%iob) + call perr(myname_,' ich =',ich(k)) + call die(myname_) + istat=-1 + _EXIT_(myname_) + return + endif + enddo + deallocate(ich) + endif +_EXIT_(myname_) + return +end subroutine obsnode_xread_ + +subroutine obsnode_xwrite_(anode,junit,jstat) + implicit none + class(colvknode),intent(in):: anode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=myname//'::obsnode_xwrite_' + integer(i_kind):: k +_ENTRY_(myname_) + + jstat=0 + write(junit,iostat=jstat) anode%nlco + if (jstat/=0) then + call perr(myname_,'write(%nlco), jstat =',jstat) + _EXIT_(myname_) + return + end if + + write(junit,iostat=jstat) (/ (anode%diags(k)%ptr%ich, k=1,anode%nlco) /), & !(nlco) + anode%res , & !(nlco) + anode%err2 , & !(nlco) + anode%raterr2, & !(nlco) + anode%prs , & !(nlevp) + anode%ipos , & !(nlco) + anode%ak , & !(nlco,nlco) + anode%ap , & !(nlco) + anode%wk , & !( nsig) + anode%wij , & !(8,nsig) + anode%ij !(4) + if (jstat/=0) then + call perr(myname_,'write(%(res,err2,...)), istat =',jstat) + _EXIT_(myname_) + return + end if +_EXIT_(myname_) + return +end subroutine obsnode_xwrite_ + +subroutine obsnode_sethop_(anode) + !use m_obsnode, only: obstype_gethop => obsnode_gethop + use m_cvgridlookup, only: cvgridlookup_getiw + use gridmod, only: nsig + implicit none + class(colvknode),intent(inout):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_sethop_' + real(r_kind),dimension(4):: zwij + integer(i_kind):: k +_ENTRY_(myname_) + call cvgridlookup_getiw(anode%elat,anode%elon,anode%ij,zwij) + do k=1,nsig + anode%wij(1:4,k) = zwij(1:4)*anode%wk(k) + anode%wij(5:8,k) = zwij(1:4)*(1._r_kind-anode%wk(k)) + enddo +_EXIT_(myname_) + return +end subroutine obsnode_sethop_ + +function obsnode_isvalid_(anode) result(isvalid_) + implicit none + logical:: isvalid_ + class(colvknode),intent(in):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_isvalid_' + integer(i_kind):: k +_ENTRY_(myname_) + isvalid_ = all( (/ (associated(anode%diags(k)%ptr), k=1,anode%nlco) /) ) +_EXIT_(myname_) + return +end function obsnode_isvalid_ + +pure subroutine gettlddp_(anode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(colvknode), intent(in):: anode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + integer(kind=i_kind):: k + do k=1,anode%nlco + tlddp = tlddp + anode%diags(k)%ptr%tldepart(jiter)*anode%diags(k)%ptr%tldepart(jiter) + enddo + if(present(nob)) nob=nob+anode%nlco + return +end subroutine gettlddp_ + +end module m_colvknode diff --git a/src/gsi/m_cvgridLookup.F90 b/src/gsi/m_cvgridlookup.F90 similarity index 60% rename from src/gsi/m_cvgridLookup.F90 rename to src/gsi/m_cvgridlookup.F90 index 89f8989c02..7acb5927e4 100644 --- a/src/gsi/m_cvgridLookup.F90 +++ b/src/gsi/m_cvgridlookup.F90 @@ -1,8 +1,8 @@ -module m_cvgridLookup +module m_cvgridlookup !$$$ subprogram documentation block ! . . . . -! subprogram: module m_cvgridLookup +! subprogram: module m_cvgridlookup ! prgmmr: j guo ! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 ! date: 2015-08-18 @@ -10,8 +10,8 @@ module m_cvgridLookup ! abstract: cv (control-vector) grid utilities for obs. partitioning ! ! program history log: -! 2015-08-18 j guo - initial implementation, from m_obsNode.F90 -! 2016-06-22 j guo - added []_sdget() to support m_latlonRange +! 2015-08-18 j guo - initial implementation, from m_obsnode.F90 +! 2016-06-22 j guo - added []_sdget() to support m_latlonrange ! . refined []_islocal() algorithm and related module ! variables, ilon_lbound, ilon_ubound, etc., to support ! regional/non-periodic subdomain grid types. @@ -34,25 +34,25 @@ module m_cvgridLookup use mpeu_util, only: stdout implicit none private ! except - public :: cvgridLookup_sdget ! get subdomain parameters - public :: cvgridLookup_islocal ! a lat-lon pair is "local" on a PE. - public :: cvgridLookup_isluse ! a lat-lon pair is "luse" on a PE - public :: cvgridLookup_getiw ! get interpolation operator (index,weight) - - interface cvgridLookup_sdget ; module procedure sdget_ ; end interface - interface cvgridLookup_islocal; module procedure islocal_; end interface - ! ... = []_local(elat,elon,myPE) - interface cvgridLookup_isluse ; module procedure isluse_ ; end interface - ! ... = []_luse(elat,elon,myPE) - interface cvgridLookup_getiw; module procedure & - get_Dij_ , & ! call []_getiw(elat,elon,ij(1:4),wij(1:4) [, debugging arguments]) - get_Dijk_, & ! call []_getiw(elat,elon,dlev,ijk(1:8),wijk(1:8) [, debugging arguments]) - get_Hij_ , & ! call []_getiw(elat,elon,ij(1:4),wij(1:4)) - get_Hijk_ ! call []_getiw(elat,elon,dlev,ijk(1:8),wijk(1:8)) + public :: cvgridlookup_sdget ! get subdomain parameters + public :: cvgridlookup_islocal ! a lat-lon pair is "local" on a pe. + public :: cvgridlookup_isluse ! a lat-lon pair is "luse" on a pe + public :: cvgridlookup_getiw ! get interpolation operator (index,weight) + + interface cvgridlookup_sdget ; module procedure sdget_ ; end interface + interface cvgridlookup_islocal; module procedure islocal_; end interface + ! ... = []_local(elat,elon,mype) + interface cvgridlookup_isluse ; module procedure isluse_ ; end interface + ! ... = []_luse(elat,elon,mype) + interface cvgridlookup_getiw; module procedure & + get_dij_ , & ! call []_getiw(elat,elon,ij(1:4),wij(1:4) [, debugging arguments]) + get_dijk_, & ! call []_getiw(elat,elon,dlev,ijk(1:8),wijk(1:8) [, debugging arguments]) + get_hij_ , & ! call []_getiw(elat,elon,ij(1:4),wij(1:4)) + get_hijk_ ! call []_getiw(elat,elon,dlev,ijk(1:8),wijk(1:8)) end interface !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - character(len=*),parameter :: myname='m_cvgridLookup' + character(len=*),parameter :: myname='m_cvgridlookup' logical,save:: isubdom_configured_=.false. integer(i_kind),allocatable,dimension(:),save:: ilatbox_lbound, ilatbox_ubound @@ -62,10 +62,10 @@ module m_cvgridLookup real(r_kind),save:: sdlat_lbnd_, sdlon_lbnd_ ! sub-domain lower-left corner in degree lat-lon real(r_kind),save:: sdlat_ubnd_, sdlon_ubnd_ ! sub-domain upper-right corner in degree lat-lon - real(r_kind),parameter:: DEG_IN_RAD = 4._r_kind*atan(1._r_kind)/180._r_kind - real(r_kind),parameter:: RAD_IN_DEG = 1._r_kind/DEG_IN_RAD + real(r_kind),parameter:: deg_in_rad = 4._r_kind*atan(1._r_kind)/180._r_kind + real(r_kind),parameter:: rad_in_deg = 1._r_kind/deg_in_rad - integer(i_kind),parameter:: iHalo = 1 + integer(i_kind),parameter:: ihalo = 1 #include "mytrace.H" #include "myassert.H" @@ -75,8 +75,8 @@ subroutine isubdom_config_() use gridmod, only: nlon, rlons, jlon1, jstart !, periodic_s use gridmod, only: nlat, rlats, ilat1, istart use gridmod, only: regional - use mpimod , only: nPEs => npe ! number of PEs in GSI_comm_world - use mpimod , only: myPE ! my rank in GSI_comm_world + use mpimod , only: npes => npe ! number of pes in gsi_comm_world + use mpimod , only: mype ! my rank in gsi_comm_world implicit none character(len=*),parameter:: myname_=myname//"::isubdom_config_" @@ -86,22 +86,22 @@ subroutine isubdom_config_() !integer(i_kind):: mlon_lbound, mlon_ubound ! lower/upper lon. boundaries of the local subdomain real(r_kind):: r,rad_to_deg - rad_to_deg(r) = r*RAD_IN_DEG + rad_to_deg(r) = r*rad_in_deg ! The lower and the upper grid index boundaries of the subdomains of - ! iHalo==1, are defined in the same way as they are in obs_para(), - ! except that the array index are confined in 0:nPEs-1, instead of 1:nPEs. + ! ihalo==1, are defined in the same way as they are in obs_para(), + ! except that the array index are confined in 0:npes-1, instead of 1:npes. - allocate(ilatbox_lbound(0:nPEs-1),ilatbox_ubound(0:nPEs-1)) - allocate(ilonbox_lbound(0:nPEs-1),ilonbox_ubound(0:nPEs-1)) + allocate(ilatbox_lbound(0:npes-1),ilatbox_ubound(0:npes-1)) + allocate(ilonbox_lbound(0:npes-1),ilonbox_ubound(0:npes-1)) ! It was assumed in this isubdom_ module, the halo grid is +/-1 grid - ! line. The current implementation let a single iHalo be a configurable + ! line. The current implementation let a single ihalo be a configurable ! parameter for both latitude and longitude grids. However, no test has - ! been done for iHalo>1, which would need some changes to other part of - ! the GSI to run. + ! been done for ihalo>1, which would need some changes to other part of + ! the gsi to run. - do ipe=0,nPEs-1 + do ipe=0,npes-1 ! However, the periodicity variable periodic_s(:) are not used here, ! since its true definition is unknown. For example, they are undefined @@ -109,49 +109,49 @@ subroutine isubdom_config_() ! definitions for some regional/non-periodic longitude grids (expected ! to be .false.?) are also suspecious. - ip1=ipe+1 + ip1=ipe+1 ! latitude grid box index range of the local subdomain - ilatbox_lbound(ipe) = istart(ip1)-iHalo ! lower grid box boundary - ilatbox_ubound(ipe) = istart(ip1)+ilat1(ip1)-1 - ilatbox_ubound(ipe) = ilatbox_ubound(ipe)-1 +iHalo ! upper grid box boundary without halo + ilatbox_lbound(ipe) = istart(ip1)-ihalo ! lower grid box boundary + ilatbox_ubound(ipe) = istart(ip1)+ilat1(ip1)-1 + ilatbox_ubound(ipe) = ilatbox_ubound(ipe)-1 +ihalo ! upper grid box boundary without halo ! longitude grid box index range of the local subdomain - ilonbox_lbound(ipe) = jstart(ip1)-iHalo ! lower grid box boundary - ilonbox_ubound(ipe) = jstart(ip1)+jlon1(ip1)-1 - ilonbox_ubound(ipe) = ilonbox_ubound(ipe)-1 +iHalo ! upper grid box boundary without halo + ilonbox_lbound(ipe) = jstart(ip1)-ihalo ! lower grid box boundary + ilonbox_ubound(ipe) = jstart(ip1)+jlon1(ip1)-1 + ilonbox_ubound(ipe) = ilonbox_ubound(ipe)-1 +ihalo ! upper grid box boundary without halo enddo - ip1=myPE+1 + ip1=mype+1 sdlat_ref_ = rad_to_deg(rlats(istart(ip1))) sdlon_ref_ = rad_to_deg(rlons(jstart(ip1))) ! Map local grid box range into local grid point range ! First, for the local latitude grid, the lower-bound and the upper-bound - ilat= ilatbox_lbound(myPE) ! the lower grid point is the lower grid box + ilat= ilatbox_lbound(mype) ! the lower grid point is the lower grid box ilat= min(max(1,ilat),nlat) ! map it to a valid grid index sdlat_lbnd_= rad_to_deg(rlats(ilat)) ! get the grid value - ilat= ilatbox_ubound(myPE)+1 ! the upper grid point is the outer edge of the upper grid box + ilat= ilatbox_ubound(mype)+1 ! the upper grid point is the outer edge of the upper grid box ilat= min(max(1,ilat),nlat) ! map it to a valid grid index sdlat_ubnd_= rad_to_deg(rlats(ilat)) ! get the grid value ! Then, for the local longitude grid, the lower-bound and the upper-bound - ilon=ilonbox_lbound(myPE) ! the lower grid point is the lower grid box + ilon=ilonbox_lbound(mype) ! the lower grid point is the lower grid box if(regional) then - ilon=min(max(1,ilon),nlon) ! map it to a valid non-periodic grid index + ilon=min(max(1,ilon),nlon) ! map it to a valid non-periodic grid index else - if(ilon<1) ilon= nlon+ilon ! map it to a valid periodic grid index + if(ilon<1) ilon= nlon+ilon ! map it to a valid periodic grid index endif sdlon_lbnd_= rad_to_deg(rlons(ilon)) ! get the grid value - ilon=ilonbox_ubound(myPE)+1 ! the upper grid point is the outer edge of the upper grid box + ilon=ilonbox_ubound(mype)+1 ! the upper grid point is the outer edge of the upper grid box if(regional) then - ilon=min(max(1,ilon),nlon) ! map it to a valid non-periodic grid index + ilon=min(max(1,ilon),nlon) ! map it to a valid non-periodic grid index else - if(ilon>nlon) ilon= ilon-nlon ! map it to a valid periodic grid index + if(ilon>nlon) ilon= ilon-nlon ! map it to a valid periodic grid index endif sdlon_ubnd_= rad_to_deg(rlons(ilon)) ! get the grid value @@ -219,14 +219,14 @@ subroutine sdget_( sdlat_ref, sdlat_lbnd, sdlat_ubnd, & end subroutine sdget_ !----- external routine: determine if (elat,elon) is local on the cv-grid -function islocal_(elat,elon,iPE) +function islocal_(elat,elon,ipe) use gridmod, only: nlat use gridmod, only: nlon !,periodic_s use gridmod, only: regional implicit none logical :: islocal_ real (r_kind), intent(in):: elat,elon - integer(i_kind), intent(in):: iPE + integer(i_kind), intent(in):: ipe integer(kind=i_kind):: ilat,ilon _ENTRY_('isluse_') @@ -234,8 +234,8 @@ function islocal_(elat,elon,iPE) call isubdom_index_(elat,elon,ilat,ilon) ASSERT( allocated(ilonbox_lbound )) - ASSERT(iPE>=lbound(ilonbox_lbound,1)) - ASSERT(iPE<=ubound(ilonbox_lbound,1)) + ASSERT(ipe>=lbound(ilonbox_lbound,1)) + ASSERT(ipe<=ubound(ilonbox_lbound,1)) ! The conditions defining "local" were original kept the same as they ! were in obs_para(). However, the way it uses periodic_s(:) was @@ -245,22 +245,22 @@ function islocal_(elat,elon,iPE) ! configuration variable _regional_ instead. ilat=min(max(1,ilat),nlat) ! ilat==nlat should not happen. See grdcrd() through isubdom_index_() - islocal_ = ilat>=ilatbox_lbound(iPE) .and. ilat<=ilatbox_ubound(iPE) + islocal_ = ilat>=ilatbox_lbound(ipe) .and. ilat<=ilatbox_ubound(ipe) if(islocal_) then - ilon=min(max(0,ilon),nlon) ! ilon==0 should not happen. See grdcrd() through isubdom_index_() - ! ilon==nlon should happend only for periodic grid - islocal_ = ilon>=ilonbox_lbound(iPE) .and. ilon<=ilonbox_ubound(iPE) - if(.not.islocal_) then - if(.not.regional) then ! periodicity applies - islocal_ = ilon-nlon >= ilonbox_lbound(iPE) .or. & ! left to grid box 1 - ilon+nlon <= ilonbox_ubound(iPE) ! right to grid box nlon-1 - endif - endif - - !-- This is the original algorithm as it is defined in obs_para() --- - !islocal_ = ( ilon>=ilon_west(iPE) .and. ilon<=ilon_east(iPE) ) .or. & - ! ( ilon==0 .and. ilon_east(iPE)>=nlon ) .or. & - ! ( ilon==nlon .and. ilon_west(iPE)<=1 ) .or. ilon_periodic(iPE) + ilon=min(max(0,ilon),nlon) ! ilon==0 should not happen. See grdcrd() through isubdom_index_() + ! ilon==nlon should happend only for periodic grid + islocal_ = ilon>=ilonbox_lbound(ipe) .and. ilon<=ilonbox_ubound(ipe) + if(.not.islocal_) then + if(.not.regional) then ! periodicity applies + islocal_ = ilon-nlon >= ilonbox_lbound(ipe) .or. & ! left to grid box 1 + ilon+nlon <= ilonbox_ubound(ipe) ! right to grid box nlon-1 + endif + endif + + !-- This is the original algorithm as it is defined in obs_para() --- + !islocal_ = ( ilon>=ilon_west(ipe) .and. ilon<=ilon_east(ipe) ) .or. & + ! ( ilon==0 .and. ilon_east(ipe)>=nlon ) .or. & + ! ( ilon==nlon .and. ilon_west(ipe)<=1 ) .or. ilon_periodic(ipe) endif _EXIT_('isluse_') end function islocal_ @@ -271,85 +271,85 @@ function isluse_(elat,elon,mype) result(luse) implicit none logical:: luse real (r_kind), intent(in ):: elat,elon - integer(i_kind), intent(in ):: myPE + integer(i_kind), intent(in ):: mype call setluse_(luse,elat,elon,mype) end function isluse_ subroutine setluse_(luse,elat,elon,mype) implicit none logical , intent(out):: luse real (r_kind), intent(in ):: elat,elon - integer(i_kind), intent(in ):: myPE + integer(i_kind), intent(in ):: mype integer(kind=i_kind):: ipe _ENTRY_('setluse_') ! Similar to what it was in obs_para() (and almost the same lately), for ! a given location (elat,elon), luse is set to .true., if and only if, - ! islocal(myPE) .and. .not.any(islocal_(myPE-1:0:-1)) + ! islocal(mype) .and. .not.any(islocal_(mype-1:0:-1)) luse=islocal_(elat,elon,mype) - ! It would not be here, if it were not even islocal() on the current PE. + ! It would not be here, if it were not even islocal() on the current pe. ! So does this assertion. ASSERT(luse) do ipe=mype-1,0,-1 - ! On any PE of lower rank (iPE = mype-1, mype-2, ..., 0), if this ob + ! On any pe of lower rank (ipe = mype-1, mype-2, ..., 0), if this ob ! islocal(), then it is .not. luse. - if(islocal_(elat,elon,iPE)) then - luse=.false. - exit - endif + if(islocal_(elat,elon,ipe)) then + luse=.false. + exit + endif enddo _EXIT_('setluse_') end subroutine setluse_ -subroutine get_Hij_(elat,elon,ij,wij) +subroutine get_hij_(elat,elon,ij,wij) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use gridmod, only: get_ij use constants, only: deg2rad - use mpimod, only: myPE + use mpimod, only: mype implicit none real(r_kind),intent(in):: elat,elon integer(i_kind),dimension(:),intent(inout):: ij real (r_kind),dimension(:),intent(inout):: wij real (r_kind):: dlat,dlon -_ENTRY_('get_Hij_') +_ENTRY_('get_hij_') dlat=elat*deg2rad dlon=elon*deg2rad call grdcrd1(dlat,rlats,nlat,1) ! returns dlat:=. call grdcrd1(dlon,rlons,nlon,1) ! returns dlon:=. - call get_ij(myPE+1,dlat,dlon,ij(1:4),wij(1:4)) -_EXIT_('get_Hij_') -end subroutine get_Hij_ + call get_ij(mype+1,dlat,dlon,ij(1:4),wij(1:4)) +_EXIT_('get_hij_') +end subroutine get_hij_ -subroutine get_Hijk_(elat,elon,dlev,ijk,wijk) +subroutine get_hijk_(elat,elon,dlev,ijk,wijk) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use gridmod, only: get_ijk use constants, only: deg2rad - use mpimod, only: myPE + use mpimod, only: mype implicit none real (r_kind),intent(in):: elat,elon,dlev integer(i_kind),dimension(:),intent(inout):: ijk real (r_kind),dimension(:),intent(inout):: wijk real (r_kind):: dlat,dlon -_ENTRY_('get_Hijk_') +_ENTRY_('get_hijk_') dlat=elat*deg2rad dlon=elon*deg2rad call grdcrd1(dlat,rlats,nlat,1) ! returns dlat:=. call grdcrd1(dlon,rlons,nlon,1) ! returns dlon:=. - call get_ijk(myPE+1,dlat,dlon,dlev,ijk(1:8),wijk(1:8)) -_EXIT_('get_Hijk_') -end subroutine get_Hijk_ + call get_ijk(mype+1,dlat,dlon,dlev,ijk(1:8),wijk(1:8)) +_EXIT_('get_hijk_') +end subroutine get_hijk_ -subroutine get_Dij_(elat,elon,dlat,dlon,ij,wij,jtype,asis) +subroutine get_dij_(elat,elon,dlat,dlon,ij,wij,jtype,asis) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use gridmod, only: get_ij use constants, only: deg2rad - use mpimod, only: myPE + use mpimod, only: mype implicit none real(r_kind),intent(in):: elat,elon real(r_kind),intent(in):: dlat,dlon @@ -364,12 +364,12 @@ subroutine get_Dij_(elat,elon,dlat,dlon,ij,wij,jtype,asis) ! integer(i_kind):: ilat,ilon ! real (r_kind):: wlat,wlon logical:: asis_ -_ENTRY_('get_Dij_') +_ENTRY_('get_dij_') asis_=.false. if(present(asis)) asis_=asis if(asis_) then - _EXIT_('get_Dij_') - return + _EXIT_('get_dij_') + return endif dlat_=elat*deg2rad @@ -379,42 +379,42 @@ subroutine get_Dij_(elat,elon,dlat,dlon,ij,wij,jtype,asis) ij_(1:4)= ij(1:4) wij_(1:4)=wij_(1:4) - call get_ij(myPE+1,dlat_,dlon_,ij(1:4),wij(1:4)) -! call get_ij(myPE+1,dlat_,dlon_,ij(1:4),wij(1:4), & + call get_ij(mype+1,dlat_,dlon_,ij(1:4),wij(1:4)) +! call get_ij(mype+1,dlat_,dlon_,ij(1:4),wij(1:4), & ! jjlat=ilat,jjlon=ilon,wwlat=wlat,wwlon=wlon,verbose=.true.) if(.not.all(ij_(1:4)==ij(1:4))) then - call tell('get_Dij_','all( ij_(1:4)== ij(1:4)) =',all( ij_(1:4)== ij(1:4))) - call tell('get_Dij_','all(wij_(1:4)==wij(1:4)) =',all(wij_(1:4)==wij(1:4))) - call tell('get_Dij_',' jtype =',jtype) - call tell('get_Dij_',' elat =',elat ) - call tell('get_Dij_',' dlat =',dlat_) - call tell('get_Dij_',' dlat0=',dlat ) - call tell('get_Dij_',' dlat-dlat0=',dlat_-dlat) - call tell('get_Dij_',' elon =',elon ) - call tell('get_Dij_',' dlon =',dlon_) - call tell('get_Dij_',' dlon0=',dlon ) - call tell('get_Dij_',' dlon-dlon0=',dlon_-dlon) -! call tell('get_Dij_',' ilat =',ilat) -! call tell('get_Dij_',' ilon =',ilon) -! call tell('get_Dij_',' wlat =',wlat) -! call tell('get_Dij_',' wlon =',wlon) - write(stdout,'(a,8i8 )') 'get_Dij(): ij (1:4) =', ij (1:4) - write(stdout,'(a,8i8 )') 'get_Dij(): ij~(1:4) =', ij_(1:4) - write(stdout,'(a,8f8.4)') 'get_Dij(): wij (1:4) =',wij (1:4) - write(stdout,'(a,8f8.4)') 'get_Dij(): wij~(1:4) =',wij_(1:4) - call die('get_Dij_') + call tell('get_dij_','all( ij_(1:4)== ij(1:4)) =',all( ij_(1:4)== ij(1:4))) + call tell('get_dij_','all(wij_(1:4)==wij(1:4)) =',all(wij_(1:4)==wij(1:4))) + call tell('get_dij_',' jtype =',jtype) + call tell('get_dij_',' elat =',elat ) + call tell('get_dij_',' dlat =',dlat_) + call tell('get_dij_',' dlat0=',dlat ) + call tell('get_dij_',' dlat-dlat0=',dlat_-dlat) + call tell('get_dij_',' elon =',elon ) + call tell('get_dij_',' dlon =',dlon_) + call tell('get_dij_',' dlon0=',dlon ) + call tell('get_dij_',' dlon-dlon0=',dlon_-dlon) +! call tell('get_dij_',' ilat =',ilat) +! call tell('get_dij_',' ilon =',ilon) +! call tell('get_dij_',' wlat =',wlat) +! call tell('get_dij_',' wlon =',wlon) + write(stdout,'(a,8i8 )') 'get_dij(): ij (1:4) =', ij (1:4) + write(stdout,'(a,8i8 )') 'get_dij(): ij~(1:4) =', ij_(1:4) + write(stdout,'(a,8f8.4)') 'get_dij(): wij (1:4) =',wij (1:4) + write(stdout,'(a,8f8.4)') 'get_dij(): wij~(1:4) =',wij_(1:4) + call die('get_dij_') endif ! ASSERT(all( ij_(1:4)== ij(1:4))) ! ASSERT(all(wij_(1:4)==wij(1:4))) -_EXIT_('get_Dij_') -end subroutine get_Dij_ +_EXIT_('get_dij_') +end subroutine get_dij_ -subroutine get_Dijk_(elat,elon,dlat,dlon,dlev,ijk,wijk,jtype,asis) +subroutine get_dijk_(elat,elon,dlat,dlon,dlev,ijk,wijk,jtype,asis) use gridmod, only: nlat,rlats use gridmod, only: nlon,rlons use gridmod, only: get_ijk use constants, only: deg2rad - use mpimod, only: myPE + use mpimod, only: mype implicit none real (r_kind),intent(in):: elat,elon,dlev real (r_kind),intent(in):: dlat,dlon @@ -427,12 +427,12 @@ subroutine get_Dijk_(elat,elon,dlat,dlon,dlev,ijk,wijk,jtype,asis) real (r_kind),dimension(8):: wijk_ real (r_kind):: dlat_,dlon_ logical:: asis_ -_ENTRY_('get_Dijk_') +_ENTRY_('get_dijk_') asis_=.false. if(present(asis)) asis_=asis if(asis_) then - _EXIT_('get_Dijk_') - return + _EXIT_('get_dijk_') + return endif dlat_=elat*deg2rad @@ -441,30 +441,30 @@ subroutine get_Dijk_(elat,elon,dlat,dlon,dlev,ijk,wijk,jtype,asis) call grdcrd1(dlon_,rlons,nlon,1) ! returns dlon:=. ijk_(1:8)= ijk(1:8) wijk_(1:8)=wijk(1:8) - call get_ijk(myPE+1,dlat_,dlon_,dlev,ijk(1:8),wijk(1:8)) + call get_ijk(mype+1,dlat_,dlon_,dlev,ijk(1:8),wijk(1:8)) if(.not.all(ijk_(1:8)==ijk(1:8))) then - call tell('get_Dijk_','all( ijk_(1:8)== ijk(1:8)) =',all( ijk_(1:8)== ijk(1:8))) - call tell('get_Dijk_','all(wijk_(1:8)==wijk(1:8)) =',all(wijk_(1:8)==wijk(1:8))) - call tell('get_Dijk_',' jtype =',jtype) - call tell('get_DijK_',' elat =',elat ) - call tell('get_Dijk_',' dlat =',dlat_) - call tell('get_Dijk_',' dlat0=',dlat ) - call tell('get_Dijk_',' dlat-dlat0=',dlat_-dlat) - call tell('get_DijK_',' elon =',elon ) - call tell('get_Dijk_',' dlon =',dlon_) - call tell('get_Dijk_',' dlon0=',dlon ) - call tell('get_Dijk_',' dlon-dlon0=',dlon_-dlon) - call tell('get_Dijk_',' dlev =',dlev ) - write(stdout,'(a,8i8 )') 'get_Dijk(): ijk (1:8) =', ijk (1:8) - write(stdout,'(a,8i8 )') 'get_Dijk(): ijk~(1:8) =', ijk_(1:8) - write(stdout,'(a,8f8.4)') 'get_Dijk(): wijk (1:8) =',wijk (1:8) - write(stdout,'(a,8f8.4)') 'get_Dijk(): wijk~(1:8) =',wijk_(1:8) - call die('get_Dijk_') + call tell('get_dijk_','all( ijk_(1:8)== ijk(1:8)) =',all( ijk_(1:8)== ijk(1:8))) + call tell('get_dijk_','all(wijk_(1:8)==wijk(1:8)) =',all(wijk_(1:8)==wijk(1:8))) + call tell('get_dijk_',' jtype =',jtype) + call tell('get_dijK_',' elat =',elat ) + call tell('get_dijk_',' dlat =',dlat_) + call tell('get_dijk_',' dlat0=',dlat ) + call tell('get_dijk_',' dlat-dlat0=',dlat_-dlat) + call tell('get_dijK_',' elon =',elon ) + call tell('get_dijk_',' dlon =',dlon_) + call tell('get_dijk_',' dlon0=',dlon ) + call tell('get_dijk_',' dlon-dlon0=',dlon_-dlon) + call tell('get_dijk_',' dlev =',dlev ) + write(stdout,'(a,8i8 )') 'get_dijk(): ijk (1:8) =', ijk (1:8) + write(stdout,'(a,8i8 )') 'get_dijk(): ijk~(1:8) =', ijk_(1:8) + write(stdout,'(a,8f8.4)') 'get_dijk(): wijk (1:8) =',wijk (1:8) + write(stdout,'(a,8f8.4)') 'get_dijk(): wijk~(1:8) =',wijk_(1:8) + call die('get_dijk_') endif ! ASSERT(all( ijk_(1:8)/= 0)) ! ASSERT(all( ijk_(1:8)== ijk(1:8))) ! ASSERT(all(wijk_(1:8)==wijk(1:8))) -_EXIT_('get_Dijk_') -end subroutine get_Dijk_ +_EXIT_('get_dijk_') +end subroutine get_dijk_ -end module m_cvgridLookup +end module m_cvgridlookup diff --git a/src/gsi/m_dbzNode.F90 b/src/gsi/m_dbzNode.F90 deleted file mode 100644 index 6a80f0df8d..0000000000 --- a/src/gsi/m_dbzNode.F90 +++ /dev/null @@ -1,264 +0,0 @@ -module m_dbzNode -!$$$ subprogram documentation block -! . . . . -! subprogram: module m_dbzNode -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 -! date: 2016-05-18 -! -! abstract: class-module of obs-type dbzNode (radar reflectivity) -! -! program history log: -! 2016-05-18 j guo - added this document block for the initial polymorphic -! implementation. -! 2017-05-12 Y. Wang and X. Wang - module for defining reflectivity observation, -! POC: xuguang.wang@ou.edu -! -! 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_obsdiagNode, only: obs_diag - use m_obsdiagNode, only: obs_diags - use kinds , only: i_kind,r_kind - use mpeu_util, only: assert_,die,perr,warn,tell - use m_obsNode, only: obsNode - implicit none - private - - public:: dbzNode - - type,extends(obsNode):: dbzNode -! type(dbz_ob_type),pointer :: llpoint => NULL() - type(obs_diag), pointer :: diags => NULL() - real(r_kind) :: res ! Radar reflectivity residual (ob-ges) - real(r_kind) :: err2 ! radar reflectivity error squared - real(r_kind) :: raterr2 ! square of ratio of final obs error - ! to original obs error - real(r_kind) :: jqr ! for TL and ADJ !modified: - real(r_kind) :: jqs ! for TL and ADJ - !real(r_kind) :: jqi ! for TL and ADJ - real(r_kind) :: jqg ! for TL and ADJ - !real(r_kind) :: jnr ! for TL and ADJ - !real(r_kind) :: jni ! for TL and ADJ - real(r_kind) :: jqli ! for TL and ADJ - !real(r_kind) :: time ! observation time in sec - real(r_kind) :: b ! variational quality control parameter - real(r_kind) :: pg ! variational quality control parameter - real(r_kind) :: wij(8) ! horizontal interpolation weights - integer(i_kind) :: ij(8) ! horizontal locations -! logical :: luse ! flag indicating if ob is used in pen. - -! integer(i_kind) :: idv,iob ! device id and obs index for sorting - - real (r_kind) :: dlev ! reference to the vertical grid - contains - procedure,nopass:: mytype - procedure:: setHop => obsNode_setHop_ - procedure:: xread => obsNode_xread_ - procedure:: xwrite => obsNode_xwrite_ - procedure:: isvalid => obsNode_isvalid_ - procedure:: gettlddp => gettlddp_ - - ! procedure, nopass:: headerRead => obsHeader_read_ - ! procedure, nopass:: headerWrite => obsHeader_write_ - ! procedure:: init => obsNode_init_ - ! procedure:: clean => obsNode_clean_ - end type dbzNode - - public:: dbzNode_typecast - public:: dbzNode_nextcast - interface dbzNode_typecast; module procedure typecast_ ; end interface - interface dbzNode_nextcast; module procedure nextcast_ ; end interface - - public:: dbzNode_appendto - interface dbzNode_appendto; module procedure appendto_ ; end interface - - character(len=*),parameter:: MYNAME="m_dbzNode" - -!#define CHECKSUM_VERBOSE -!#define DEBUG_TRACE -#include "myassert.H" -#include "mytrace.H" -contains -function typecast_(aNode) result(ptr_) -!-- cast a class(obsNode) to a type(dbzNode) - use m_obsNode, only: obsNode - implicit none - type (dbzNode),pointer:: ptr_ - class(obsNode),pointer,intent(in):: aNode - ptr_ => null() - if(.not.associated(aNode)) return - ! logically, typecast of a null-reference is a null pointer - select type(aNode) - type is(dbzNode) - ptr_ => aNode - end select -return -end function typecast_ - -function nextcast_(aNode) result(ptr_) -!-- cast an obsNode_next(obsNode) to a type(dbzNode) - use m_obsNode, only: obsNode,obsNode_next - implicit none - type (dbzNode),pointer:: ptr_ - class(obsNode),target ,intent(in):: aNode - - class(obsNode),pointer:: inode_ - inode_ => obsNode_next(aNode) - ptr_ => typecast_(inode_) -return -end function nextcast_ - -subroutine appendto_(aNode,oll) -!-- append aNode to linked-list oLL - use m_obsNode , only: obsNode - use m_obsLList, only: obsLList,obsLList_appendNode - implicit none - type(dbzNode),pointer,intent(in):: aNode - type(obsLList),intent(inout):: oLL - - class(obsNode),pointer:: inode_ - inode_ => aNode - call obsLList_appendNode(oLL,inode_) - inode_ => null() -end subroutine appendto_ - -! obsNode implementations - -function mytype() - implicit none - character(len=:),allocatable:: mytype - mytype="[dbzNode]" -end function mytype - -subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) - use m_obsdiagNode, only: obsdiagLookup_locate - implicit none - class(dbzNode),intent(inout):: aNode - integer(i_kind),intent(in ):: iunit - integer(i_kind),intent( out):: istat - type(obs_diags),intent(in ):: diagLookup - logical,optional,intent(in ):: skip - - character(len=*),parameter:: myname_=MYNAME//'.obsNode_xread_' - logical:: skip_ -_ENTRY_(myname_) - skip_=.false. - if(present(skip)) skip_=skip - - istat=0 - if(skip_) then - read(iunit,iostat=istat) - if(istat/=0) then - call perr(myname_,'skipping read(%(res,err2,...)), iostat =',istat) - _EXIT_(myname_) - return - endif - - else - read(iunit,iostat=istat) aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%b , & - aNode%pg , & - aNode%jqr , & - aNode%jqs , & - aNode%jqg , & - aNode%jqli , & - aNode%dlev , & - aNode%wij , & - aNode%ij - if (istat/=0) then - call perr(myname_,'read(%(res,err2,...)), iostat =',istat) - _EXIT_(myname_) - return - end if - - aNode%diags => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,1_i_kind) - if(.not.associated(aNode%diags)) then - call perr(myname_,'obsdiagLookup_locate(), %idv =',aNode%idv) - call perr(myname_,' %iob =',aNode%iob) - call die(myname_) - endif - endif -_EXIT_(myname_) -return -end subroutine obsNode_xread_ - -subroutine obsNode_xwrite_(aNode,junit,jstat) - implicit none - class(dbzNode),intent(in):: aNode - integer(i_kind),intent(in ):: junit - integer(i_kind),intent( out):: jstat - - character(len=*),parameter:: myname_=MYNAME//'.obsNode_xwrite_' -_ENTRY_(myname_) - - jstat=0 - write(junit,iostat=jstat) aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%b , & - aNode%pg , & - aNode%jqr , & - aNode%jqs , & - aNode%jqg , & - aNode%jqli , & - aNode%dlev , & - aNode%wij , & - aNode%ij - if (jstat/=0) then - call perr(myname_,'write(%(res,err2,...)), iostat =',jstat) - _EXIT_(myname_) - return - end if -_EXIT_(myname_) -return -end subroutine obsNode_xwrite_ - -subroutine obsNode_setHop_(aNode) - use m_cvgridLookup, only: cvgridLookup_getiw - implicit none - class(dbzNode),intent(inout):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_' -_ENTRY_(myname_) - call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%dlev,aNode%ij,aNode%wij) -_EXIT_(myname_) -return -end subroutine obsNode_setHop_ - -function obsNode_isvalid_(aNode) result(isvalid_) - implicit none - logical:: isvalid_ - class(dbzNode),intent(in):: aNode - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_' -_ENTRY_(myname_) - isvalid_=associated(aNode%diags) -_EXIT_(myname_) -return -end function obsNode_isvalid_ - -pure subroutine gettlddp_(aNode,jiter,tlddp,nob) - use kinds, only: r_kind - implicit none - class(dbzNode), intent(in):: aNode - integer(kind=i_kind),intent(in):: jiter - real(kind=r_kind),intent(inout):: tlddp - integer(kind=i_kind),optional,intent(inout):: nob - - tlddp = tlddp + aNode%diags%tldepart(jiter)*aNode%diags%tldepart(jiter) - if(present(nob)) nob=nob+1 -return -end subroutine gettlddp_ - -end module m_dbzNode diff --git a/src/gsi/m_dbznode.F90 b/src/gsi/m_dbznode.F90 new file mode 100644 index 0000000000..59971ce2b5 --- /dev/null +++ b/src/gsi/m_dbznode.F90 @@ -0,0 +1,264 @@ +module m_dbznode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_dbznode +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2016-05-18 +! +! abstract: class-module of obs-type dbznode (radar reflectivity) +! +! program history log: +! 2016-05-18 j guo - added this document block for the initial polymorphic +! implementation. +! 2017-05-12 Y. Wang and X. Wang - module for defining reflectivity observation, +! POC: xuguang.wang@ou.edu +! +! 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_obsdiagnode, only: obs_diag + use m_obsdiagnode, only: obs_diags + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsnode, only: obsnode + implicit none + private + + public:: dbznode + + type,extends(obsnode):: dbznode +! type(dbz_ob_type),pointer :: llpoint => null() + type(obs_diag), pointer :: diags => null() + real(r_kind) :: res ! Radar reflectivity residual (ob-ges) + real(r_kind) :: err2 ! radar reflectivity error squared + real(r_kind) :: raterr2 ! square of ratio of final obs error + ! to original obs error + real(r_kind) :: jqr ! for TL and ADJ !modified: + real(r_kind) :: jqs ! for TL and ADJ + !real(r_kind) :: jqi ! for TL and ADJ + real(r_kind) :: jqg ! for TL and ADJ + !real(r_kind) :: jnr ! for TL and ADJ + !real(r_kind) :: jni ! for TL and ADJ + real(r_kind) :: jqli ! for TL and ADJ + !real(r_kind) :: time ! observation time in sec + real(r_kind) :: b ! variational quality control parameter + real(r_kind) :: pg ! variational quality control parameter + real(r_kind) :: wij(8) ! horizontal interpolation weights + integer(i_kind) :: ij(8) ! horizontal locations +! logical :: luse ! flag indicating if ob is used in pen. + +! integer(i_kind) :: idv,iob ! device id and obs index for sorting + + real (r_kind) :: dlev ! reference to the vertical grid + contains + procedure,nopass:: mytype + procedure:: sethop => obsnode_sethop_ + procedure:: xread => obsnode_xread_ + procedure:: xwrite => obsnode_xwrite_ + procedure:: isvalid => obsnode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerread => obsheader_read_ + ! procedure, nopass:: headerwrite => obsheader_write_ + ! procedure:: init => obsnode_init_ + ! procedure:: clean => obsnode_clean_ + end type dbznode + + public:: dbznode_typecast + public:: dbznode_nextcast + interface dbznode_typecast; module procedure typecast_ ; end interface + interface dbznode_nextcast; module procedure nextcast_ ; end interface + + public:: dbznode_appendto + interface dbznode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: myname="m_dbznode" + +!#define CHECKSUM_VERBOSE +!#define DEBUG_TRACE +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(anode) result(ptr_) +!-- cast a class(obsnode) to a type(dbznode) + use m_obsnode, only: obsnode + implicit none + type (dbznode),pointer:: ptr_ + class(obsnode),pointer,intent(in):: anode + ptr_ => null() + if(.not.associated(anode)) return + ! logically, typecast of a null-reference is a null pointer + select type(anode) + type is(dbznode) + ptr_ => anode + end select + return +end function typecast_ + +function nextcast_(anode) result(ptr_) +!-- cast an obsnode_next(obsnode) to a type(dbznode) + use m_obsnode, only: obsnode,obsnode_next + implicit none + type (dbznode),pointer:: ptr_ + class(obsnode),target ,intent(in):: anode + + class(obsnode),pointer:: inode_ + inode_ => obsnode_next(anode) + ptr_ => typecast_(inode_) + return +end function nextcast_ + +subroutine appendto_(anode,oll) +!-- append anode to linked-list oll + use m_obsnode , only: obsnode + use m_obsllist, only: obsllist,obsllist_appendnode + implicit none + type(dbznode),pointer,intent(in):: anode + type(obsllist),intent(inout):: oll + + class(obsnode),pointer:: inode_ + inode_ => anode + call obsllist_appendnode(oll,inode_) + inode_ => null() +end subroutine appendto_ + +! obsnode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[dbznode]" +end function mytype + +subroutine obsnode_xread_(anode,iunit,istat,diaglookup,skip) + use m_obsdiagnode, only: obsdiaglookup_locate + implicit none + class(dbznode),intent(inout):: anode + integer(i_kind),intent(in ):: iunit + integer(i_kind),intent( out):: istat + type(obs_diags),intent(in ):: diaglookup + logical,optional,intent(in ):: skip + + character(len=*),parameter:: myname_=myname//'.obsnode_xread_' + logical:: skip_ +_ENTRY_(myname_) + skip_=.false. + if(present(skip)) skip_=skip + + istat=0 + if(skip_) then + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(res,err2,...)), iostat =',istat) + _EXIT_(myname_) + return + endif + + else + read(iunit,iostat=istat) anode%res , & + anode%err2 , & + anode%raterr2, & + anode%b , & + anode%pg , & + anode%jqr , & + anode%jqs , & + anode%jqg , & + anode%jqli , & + anode%dlev , & + anode%wij , & + anode%ij + if (istat/=0) then + call perr(myname_,'read(%(res,err2,...)), iostat =',istat) + _EXIT_(myname_) + return + end if + + anode%diags => obsdiaglookup_locate(diaglookup,anode%idv,anode%iob,1) + if(.not.associated(anode%diags)) then + call perr(myname_,'obsdiaglookup_locate(), %idv =',anode%idv) + call perr(myname_,' %iob =',anode%iob) + call die(myname_) + endif + endif +_EXIT_(myname_) + return +end subroutine obsnode_xread_ + +subroutine obsnode_xwrite_(anode,junit,jstat) + implicit none + class(dbznode),intent(in):: anode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=myname//'.obsnode_xwrite_' +_ENTRY_(myname_) + + jstat=0 + write(junit,iostat=jstat) anode%res , & + anode%err2 , & + anode%raterr2, & + anode%b , & + anode%pg , & + anode%jqr , & + anode%jqs , & + anode%jqg , & + anode%jqli , & + anode%dlev , & + anode%wij , & + anode%ij + if (jstat/=0) then + call perr(myname_,'write(%(res,err2,...)), iostat =',jstat) + _EXIT_(myname_) + return + end if +_EXIT_(myname_) + return +end subroutine obsnode_xwrite_ + +subroutine obsnode_sethop_(anode) + use m_cvgridlookup, only: cvgridlookup_getiw + implicit none + class(dbznode),intent(inout):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_sethop_' +_ENTRY_(myname_) + call cvgridlookup_getiw(anode%elat,anode%elon,anode%dlev,anode%ij,anode%wij) +_EXIT_(myname_) + return +end subroutine obsnode_sethop_ + +function obsnode_isvalid_(anode) result(isvalid_) + implicit none + logical:: isvalid_ + class(dbznode),intent(in):: anode + + character(len=*),parameter:: myname_=myname//'::obsnode_isvalid_' +_ENTRY_(myname_) + isvalid_=associated(anode%diags) +_EXIT_(myname_) + return +end function obsnode_isvalid_ + +pure subroutine gettlddp_(anode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(dbznode), intent(in):: anode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + tlddp = tlddp + anode%diags%tldepart(jiter)*anode%diags%tldepart(jiter) + if(present(nob)) nob=nob+1 + return +end subroutine gettlddp_ + +end module m_dbznode diff --git a/src/gsi/m_dgeevx.F90 b/src/gsi/m_dgeevx.F90 index bf5b378c1e..82a944f35d 100644 --- a/src/gsi/m_dgeevx.F90 +++ b/src/gsi/m_dgeevx.F90 @@ -6,7 +6,7 @@ module m_dgeevx ! org: NASA/GSFC, Global Modeling and Assimilation Office, 900.3 ! date: 2010-03-24 ! -! abstract: an alternative interface of LAPACK DGEEV() +! abstract: an alternative interface of lapack dgeev() ! ! program history log: ! 2010-03-24 j guo - added this document block @@ -27,7 +27,7 @@ module m_dgeevx ! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS ! !BOP ------------------------------------------------------------------- ! -! !MODULE: m_dgeevx - an alternative interface of LAPACK DGEEV() +! !MODULE: m_dgeevx - an alternative interface of lapack dgeev() ! ! !DESCRIPTION: ! @@ -36,7 +36,7 @@ module m_dgeevx implicit none private ! except - public :: dgeevx ! alternative DGEEV() + public :: dgeevx ! alternative dgeev() ! !REVISION HISTORY: ! 29May08 - Jing Guo @@ -89,7 +89,7 @@ subroutine dgeevx(nsig,qmat,ldqmat,www,wwwd,zzz,ldzzz,zzzd,ldzzzd,info,mype) integer(i_kind) :: naux real(r_kind) :: aux,factor - ! Use IBM ESSL routine dgeev() + ! Use ibm essl routine dgeev() ! Below is a copy of the original code in mod_vtrans.f90_1.2. !<<<< @@ -163,10 +163,10 @@ subroutine dgeevx(nsig,qmat,ldqmat,www,wwwd,zzz,ldzzz,zzzd,ldzzzd,info,mype) info=0 #else - ! Use standard LAPACK routine dgeev() + ! Use standard lapack routine dgeev() !---------------------------------------- -! Modification for using GSI dgeev -RY +! Modification for using gsi dgeev !---------------------------------------- character*1,parameter:: jobvr='V' character*1,parameter:: jobvl='V' @@ -182,8 +182,7 @@ subroutine dgeevx(nsig,qmat,ldqmat,www,wwwd,zzz,ldzzz,zzzd,ldzzzd,info,mype) lwork=size(work) !------------------------------------ -! ryang: -! use SGI scslib dgeev subroutine +! use sgi scslib dgeev subroutine do j=1,nsig do i=1,nsig @@ -197,14 +196,14 @@ subroutine dgeevx(nsig,qmat,ldqmat,www,wwwd,zzz,ldzzz,zzzd,ldzzzd,info,mype) if (mype ==0) write (6,*) myname, ' after CALL DGEEV', 'status: info = ', info -! use Dave's array names +! use dave's array names do j=1,nsig if (wi(j) /= zero) write (6,*) & 'wrong eigen computation: create_vtrans' www(1,j)=wr(j) www(2,j)=wi(j) !----------------------------------------------------------- -!ADD explanation:!!! +!add explanation:!!! ! ! vl is zzzd !the following line is not generic, only hold when wi (j)=0.0 @@ -217,7 +216,7 @@ subroutine dgeevx(nsig,qmat,ldqmat,www,wwwd,zzz,ldzzz,zzzd,ldzzzd,info,mype) zzzd(2,k,j)=zero enddo enddo -! back to Dave's code +! back to dave's code !sort from largest to smallest eigenvalues ! sort from largest to smallest eigenvalue do j=1,nsig-1 diff --git a/src/gsi/m_distance.f90 b/src/gsi/m_distance.f90 index 17b1581a29..d3745a9370 100644 --- a/src/gsi/m_distance.f90 +++ b/src/gsi/m_distance.f90 @@ -9,8 +9,8 @@ module m_distance ! assimilation ! ! program history log: -! 1996-10-01 Joiner/Karki - initial coding from NASA/GMAO -! 2012-02-15 eliu - reformat to use in GSI +! 1996-10-01 Joiner/Karki - initial coding from nasa/gmao +! 2012-02-15 eliu - reformat to use in gsi ! ! subroutines included: ! @@ -22,8 +22,9 @@ module m_distance ! !$$$ end documentation block - use constants, only: deg2rad,half,two + use constants, only: deg2rad,zero,half,two use kinds, only: i_kind, r_kind + implicit none interface distance module procedure distance1 @@ -33,7 +34,7 @@ module m_distance contains function distance1( alat, alon, blat, blon ) result (dist) - implicit NONE + implicit none !INPUT PARAMETERS: @@ -53,7 +54,7 @@ function distance1( alat, alon, blat, blon ) result (dist) !Initialize !---------- - dist=0. + dist=zero !Compute distances !================= @@ -90,6 +91,7 @@ function distance1( alat, alon, blat, blon ) result (dist) end function distance1 function distance2( alat, alon, blat, blon ) result (dist) + implicit none real(r_kind), intent(in) :: alat real(r_kind), intent(in) :: alon diff --git a/src/gsi/m_dtime.F90 b/src/gsi/m_dtime.F90 index 367f932f8c..6342bc0f8a 100644 --- a/src/gsi/m_dtime.F90 +++ b/src/gsi/m_dtime.F90 @@ -21,7 +21,7 @@ module m_dtime ! !$$$ end subprogram documentation block -#define ZERODIFFTEST +#define zerodifftest ! module interface: @@ -128,7 +128,7 @@ subroutine dtime_check(dtime, in_curbin,in_anybin) at=at+(dtime-at)/nt in_curbin = (dtime>hrdifsig(1) .and. dtime<=hrdifsig(nfldsig)) -#ifdef ZERODIFFTEST +#ifdef zerodifftest if(hrdifsig(1)==hrdifsig_all(1)) in_curbin = in_curbin .or. dtime<=hrdifsig(1) if(hrdifsig(nfldsig)==hrdifsig_all(nfldsig_all)) in_curbin = in_curbin .or. dtime>hrdifsig(nfldsig) in_anybin = .true. @@ -140,22 +140,22 @@ subroutine dtime_check(dtime, in_curbin,in_anybin) if(in_curbin) then - nm=nm+1 - am=am+(dtime-am)/nm - return + nm=nm+1 + am=am+(dtime-am)/nm + return endif if(in_anybin) return if(dtime <= hrdifsig_all(1)) then - nl=nl+1 - al=al+(dtime-al)/nl - return + nl=nl+1 + al=al+(dtime-al)/nl + return endif if(dtime > hrdifsig_all(nfldsig_all)) then - nr=nr+1 - ar=ar+(dtime-ar)/nr - return + nr=nr+1 + ar=ar+(dtime-ar)/nr + return endif end subroutine dtime_check diff --git a/src/gsi/m_dwNode.F90 b/src/gsi/m_dwNode.F90 deleted file mode 100644 index 3e211e078f..0000000000 --- a/src/gsi/m_dwNode.F90 +++ /dev/null @@ -1,256 +0,0 @@ -module m_dwNode -!$$$ subprogram documentation block -! . . . . -! subprogram: module m_dwNode -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 -! date: 2016-05-18 -! -! abstract: class-module of obs-type dwNode (Doppler wind) -! -! program history log: -! 2016-05-18 j guo - added this document block for the initial polymorphic -! implementation. -! -! 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_obsdiagNode, only: obs_diag - use m_obsdiagNode, only: obs_diags - use kinds , only: i_kind,r_kind - use mpeu_util, only: assert_,die,perr,warn,tell - use m_obsNode, only: obsNode - implicit none - private - - public:: dwNode - - type,extends(obsNode):: dwNode - !type(dw_ob_type),pointer :: llpoint => NULL() - type(obs_diag), pointer :: diags => NULL() - real(r_kind) :: res =0._r_kind ! doppler wind residual - real(r_kind) :: err2 =0._r_kind ! radial wind error squared - real(r_kind) :: raterr2=0._r_kind ! square of ratio of final obs error - ! to original obs error - !real(r_kind) :: time ! observation time in sec - real(r_kind) :: b =0._r_kind ! variational quality control parameter - real(r_kind) :: pg =0._r_kind ! variational quality control parameter - real(r_kind) :: cosazm =0._r_kind ! v factor - real(r_kind) :: sinazm =0._r_kind ! u factor - real(r_kind) :: wij(8) =0._r_kind ! horizontal interpolation weights - integer(i_kind) :: ij(8) =0_i_kind ! horizontal locations - !logical :: luse ! flag indicating if ob is used in pen. - - !integer(i_kind) :: idv,iob ! device id and obs index for sorting - !real (r_kind) :: elat, elon ! earth lat-lon for redistribution - !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution - real (r_kind) :: dlev ! reference to the vertical grid - real (r_kind) :: factw ! factor of 10m wind - contains - procedure,nopass:: mytype - procedure:: setHop => obsNode_setHop_ - procedure:: xread => obsNode_xread_ - procedure:: xwrite => obsNode_xwrite_ - procedure:: isvalid => obsNode_isvalid_ - procedure:: gettlddp => gettlddp_ - - ! procedure, nopass:: headerRead => obsHeader_read_ - ! procedure, nopass:: headerWrite => obsHeader_write_ - ! procedure:: init => obsNode_init_ - ! procedure:: clean => obsNode_clean_ - end type dwNode - - public:: dwNode_typecast - public:: dwNode_nextcast - interface dwNode_typecast; module procedure typecast_ ; end interface - interface dwNode_nextcast; module procedure nextcast_ ; end interface - - public:: dwNode_appendto - interface dwNode_appendto; module procedure appendto_ ; end interface - - character(len=*),parameter:: MYNAME="m_dwNode" - -#include "myassert.H" -#include "mytrace.H" -contains -function typecast_(aNode) result(ptr_) -!-- cast a class(obsNode) to a type(dwNode) - use m_obsNode, only: obsNode - implicit none - type (dwNode ),pointer:: ptr_ - class(obsNode),pointer,intent(in):: aNode - ptr_ => null() - if(.not.associated(aNode)) return - ! logically, typecast of a null-reference is a null pointer. - select type(aNode) - type is(dwNode) - ptr_ => aNode - end select -return -end function typecast_ - -function nextcast_(aNode) result(ptr_) -!-- cast an obsNode_next(obsNode) to a type(dwNode) - use m_obsNode, only: obsNode,obsNode_next - implicit none - type (dwNode ),pointer:: ptr_ - class(obsNode),target ,intent(in):: aNode - - class(obsNode),pointer:: inode_ - inode_ => obsNode_next(aNode) - ptr_ => typecast_(inode_) -return -end function nextcast_ - -subroutine appendto_(aNode,oll) -!-- append aNode to linked-list oLL - use m_obsNode , only: obsNode - use m_obsLList, only: obsLList,obsLList_appendNode - implicit none - type(dwNode),pointer,intent(in):: aNode - type(obsLList),intent(inout):: oLL - - class(obsNode),pointer:: inode_ - inode_ => aNode - call obsLList_appendNode(oLL,inode_) - inode_ => null() -end subroutine appendto_ - -! obsNode implementations - -function mytype() - implicit none - character(len=:),allocatable:: mytype - mytype="[dwNode]" -end function mytype - -subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip) - use m_obsdiagNode, only: obsdiagLookup_locate - implicit none - class(dwNode) , intent(inout):: aNode - integer(i_kind) , intent(in ):: iunit - integer(i_kind) , intent( out):: istat - type(obs_diags) , intent(in ):: diagLookup - logical,optional, intent(in ):: skip - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xread_' - logical:: skip_ -_ENTRY_(myname_) - skip_=.false. - if(present(skip)) skip_=skip - - istat=0 - if(skip_) then - read(iunit,iostat=istat) - if(istat/=0) then - call perr(myname_,'skipping read(%(res,...)), istat =',istat) - _EXIT_(myname_) - return - endif - - else - read(iunit,iostat=istat) aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%b , & - aNode%pg , & - aNode%cosazm , & - aNode%sinazm , & - aNode%dlev , & - aNode%factw , & - aNode%wij , & - aNode%ij - if (istat/=0) then - call perr(myname_,'read(%(res,err2,...)), istat =',istat) - _EXIT_(myname_) - return - end if - - aNode%diags => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,1_i_kind) - if(.not.associated(aNode%diags)) then - call perr(myname_,'obsdiagLookup_locate(), %idv =',aNode%idv) - call perr(myname_,' %iob =',aNode%iob) - call die(myname_) - endif - endif -_EXIT_(myname_) -return -end subroutine obsNode_xread_ - -subroutine obsNode_xwrite_(aNode,junit,jstat) - implicit none - class(dwNode),intent(in):: aNode - integer(i_kind),intent(in ):: junit - integer(i_kind),intent( out):: jstat - - character(len=*),parameter:: myname_=MYNAME//'::obsNode_xwrite_' -_ENTRY_(myname_) - - jstat=0 - write(junit,iostat=jstat) aNode%res , & - aNode%err2 , & - aNode%raterr2, & - aNode%b , & - aNode%pg , & - aNode%cosazm , & - aNode%sinazm , & - aNode%dlev , & - aNode%factw , & - aNode%wij , & - aNode%ij - if (jstat/=0) then - call perr(myname_,'write(%(res,err2,...)), jstat =',jstat) - _EXIT_(myname_) - return - end if -_EXIT_(myname_) -return -end subroutine obsNode_xwrite_ - -subroutine obsNode_setHop_(aNode) - use m_cvgridLookup, only: cvgridLookup_getiw - implicit none - class(dwNode),intent(inout):: aNode - - character(len=*),parameter:: myname_=myname//"::obsNode_setHop_" -_ENTRY_(myname_) - call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%dlev,aNode%ij,aNode%wij) - aNode%wij(1:8) = aNode%wij(1:8)*aNode%factw -_EXIT_(myname_) -return -end subroutine obsNode_setHop_ - -function obsNode_isvalid_(aNode) result(isvalid_) - implicit none - logical:: isvalid_ - class(dwNode),intent(in):: aNode - - character(len=*),parameter:: myname_=myname//"::obsNode_isvalid_" -_ENTRY_(myname_) - isvalid_=associated(aNode%diags) -_EXIT_(myname_) -return -end function obsNode_isvalid_ - -pure subroutine gettlddp_(aNode,jiter,tlddp,nob) - use kinds, only: r_kind - implicit none - class(dwNode), intent(in):: aNode - integer(kind=i_kind),intent(in):: jiter - real(kind=r_kind),intent(inout):: tlddp - integer(kind=i_kind),optional,intent(inout):: nob - - tlddp = tlddp + aNode%diags%tldepart(jiter)*aNode%diags%tldepart(jiter) - if(present(nob)) nob=nob+1 -return -end subroutine gettlddp_ - -end module m_dwNode diff --git a/src/gsi/m_dwnode.F90 b/src/gsi/m_dwnode.F90 new file mode 100644 index 0000000000..56aea79f5d --- /dev/null +++ b/src/gsi/m_dwnode.F90 @@ -0,0 +1,256 @@ +module m_dwnode +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_dwnode +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.3 +! date: 2016-05-18 +! +! abstract: class-module of obs-type dwnode (doppler wind) +! +! program history log: +! 2016-05-18 j guo - added this document block for the initial polymorphic +! implementation. +! +! 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_obsdiagnode, only: obs_diag + use m_obsdiagnode, only: obs_diags + use kinds , only: i_kind,r_kind + use mpeu_util, only: assert_,die,perr,warn,tell + use m_obsnode, only: obsnode + implicit none + private + + public:: dwnode + + type,extends(obsnode):: dwnode + !type(dw_ob_type),pointer :: llpoint => null() + type(obs_diag), pointer :: diags => null() + real(r_kind) :: res =0._r_kind ! doppler wind residual + real(r_kind) :: err2 =0._r_kind ! radial wind error squared + real(r_kind) :: raterr2=0._r_kind ! square of ratio of final obs error + ! to original obs error + !real(r_kind) :: time ! observation time in sec + real(r_kind) :: b =0._r_kind ! variational quality control parameter + real(r_kind) :: pg =0._r_kind ! variational quality control parameter + real(r_kind) :: cosazm =0._r_kind ! v factor + real(r_kind) :: sinazm =0._r_kind ! u factor + real(r_kind) :: wij(8) =0._r_kind ! horizontal interpolation weights + integer(i_kind) :: ij(8) =0 ! horizontal locations + !logical :: luse ! flag indicating if ob is used in pen. + + !integer(i_kind) :: idv,iob ! device id and obs index for sorting + !real (r_kind) :: elat, elon ! earth lat-lon for redistribution + !real (r_kind) :: dlat, dlon ! earth lat-lon for redistribution + real (r_kind) :: dlev ! reference to the vertical grid + real (r_kind) :: factw ! factor of 10m wind + contains + procedure,nopass:: mytype + procedure:: sethop => obsnode_sethop_ + procedure:: xread => obsnode_xread_ + procedure:: xwrite => obsnode_xwrite_ + procedure:: isvalid => obsnode_isvalid_ + procedure:: gettlddp => gettlddp_ + + ! procedure, nopass:: headerread => obsheader_read_ + ! procedure, nopass:: headerwrite => obsheader_write_ + ! procedure:: init => obsnode_init_ + ! procedure:: clean => obsnode_clean_ + end type dwnode + + public:: dwnode_typecast + public:: dwnode_nextcast + interface dwnode_typecast; module procedure typecast_ ; end interface + interface dwnode_nextcast; module procedure nextcast_ ; end interface + + public:: dwnode_appendto + interface dwnode_appendto; module procedure appendto_ ; end interface + + character(len=*),parameter:: myname="m_dwnode" + +#include "myassert.H" +#include "mytrace.H" +contains +function typecast_(anode) result(ptr_) +!-- cast a class(obsnode) to a type(dwnode) + use m_obsnode, only: obsnode + implicit none + type (dwnode ),pointer:: ptr_ + class(obsnode),pointer,intent(in):: anode + ptr_ => null() + if(.not.associated(anode)) return + ! logically, typecast of a null-reference is a null pointer. + select type(anode) + type is(dwnode) + ptr_ => anode + end select + return +end function typecast_ + +function nextcast_(anode) result(ptr_) +!-- cast an obsnode_next(obsnode) to a type(dwnode) + use m_obsnode, only: obsnode,obsnode_next + implicit none + type (dwnode ),pointer:: ptr_ + class(obsnode),target ,intent(in):: anode + + class(obsnode),pointer:: inode_ + inode_ => obsnode_next(anode) + ptr_ => typecast_(inode_) + return +end function nextcast_ + +subroutine appendto_(anode,oll) +!-- append anode to linked-list oll + use m_obsnode , only: obsnode + use m_obsllist, only: obsllist,obsllist_appendnode + implicit none + type(dwnode),pointer,intent(in):: anode + type(obsllist),intent(inout):: oll + + class(obsnode),pointer:: inode_ + inode_ => anode + call obsllist_appendnode(oll,inode_) + inode_ => null() +end subroutine appendto_ + +! obsnode implementations + +function mytype() + implicit none + character(len=:),allocatable:: mytype + mytype="[dwnode]" +end function mytype + +subroutine obsnode_xread_(anode,iunit,istat,diaglookup,skip) + use m_obsdiagnode, only: obsdiaglookup_locate + implicit none + class(dwnode) , intent(inout):: anode + integer(i_kind) , intent(in ):: iunit + integer(i_kind) , intent( out):: istat + type(obs_diags) , intent(in ):: diaglookup + logical,optional, intent(in ):: skip + + character(len=*),parameter:: myname_=myname//'::obsnode_xread_' + logical:: skip_ +_ENTRY_(myname_) + skip_=.false. + if(present(skip)) skip_=skip + + istat=0 + if(skip_) then + read(iunit,iostat=istat) + if(istat/=0) then + call perr(myname_,'skipping read(%(res,...)), istat =',istat) + _EXIT_(myname_) + return + endif + + else + read(iunit,iostat=istat) anode%res , & + anode%err2 , & + anode%raterr2, & + anode%b , & + anode%pg , & + anode%cosazm , & + anode%sinazm , & + anode%dlev , & + anode%factw , & + anode%wij , & + anode%ij + if (istat/=0) then + call perr(myname_,'read(%(res,err2,...)), istat =',istat) + _EXIT_(myname_) + return + end if + + anode%diags => obsdiaglookup_locate(diaglookup,anode%idv,anode%iob,1) + if(.not.associated(anode%diags)) then + call perr(myname_,'obsdiaglookup_locate(), %idv =',anode%idv) + call perr(myname_,' %iob =',anode%iob) + call die(myname_) + endif + endif +_EXIT_(myname_) + return +end subroutine obsnode_xread_ + +subroutine obsnode_xwrite_(anode,junit,jstat) + implicit none + class(dwnode),intent(in):: anode + integer(i_kind),intent(in ):: junit + integer(i_kind),intent( out):: jstat + + character(len=*),parameter:: myname_=myname//'::obsnode_xwrite_' +_ENTRY_(myname_) + + jstat=0 + write(junit,iostat=jstat) anode%res , & + anode%err2 , & + anode%raterr2, & + anode%b , & + anode%pg , & + anode%cosazm , & + anode%sinazm , & + anode%dlev , & + anode%factw , & + anode%wij , & + anode%ij + if (jstat/=0) then + call perr(myname_,'write(%(res,err2,...)), jstat =',jstat) + _EXIT_(myname_) + return + end if +_EXIT_(myname_) + return +end subroutine obsnode_xwrite_ + +subroutine obsnode_sethop_(anode) + use m_cvgridlookup, only: cvgridlookup_getiw + implicit none + class(dwnode),intent(inout):: anode + + character(len=*),parameter:: myname_=myname//"::obsnode_sethop_" +_ENTRY_(myname_) + call cvgridlookup_getiw(anode%elat,anode%elon,anode%dlev,anode%ij,anode%wij) + anode%wij(1:8) = anode%wij(1:8)*anode%factw +_EXIT_(myname_) + return +end subroutine obsnode_sethop_ + +function obsnode_isvalid_(anode) result(isvalid_) + implicit none + logical:: isvalid_ + class(dwnode),intent(in):: anode + + character(len=*),parameter:: myname_=myname//"::obsnode_isvalid_" +_ENTRY_(myname_) + isvalid_=associated(anode%diags) +_EXIT_(myname_) + return +end function obsnode_isvalid_ + +pure subroutine gettlddp_(anode,jiter,tlddp,nob) + use kinds, only: r_kind + implicit none + class(dwnode), intent(in):: anode + integer(kind=i_kind),intent(in):: jiter + real(kind=r_kind),intent(inout):: tlddp + integer(kind=i_kind),optional,intent(inout):: nob + + tlddp = tlddp + anode%diags%tldepart(jiter)*anode%diags%tldepart(jiter) + if(present(nob)) nob=nob+1 + return +end subroutine gettlddp_ + +end module m_dwnode diff --git a/src/gsi/m_extOzone.F90 b/src/gsi/m_extOzone.F90 deleted file mode 100644 index 31cdedbebc..0000000000 --- a/src/gsi/m_extOzone.F90 +++ /dev/null @@ -1,1452 +0,0 @@ -module m_extOzone -!#define NO_EXTRA_ -!$$$ subprogram documentation block -! . . . . -! subprogram: module m_extOzone -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 900.3 -! date: 2013-09-27 -! -! abstract: a module for reading extra ozone observation data -! -! program history log: -! 2012-09-10 wargan - add OMI text option for OMI with efficiency factors -! 2012-12-06 jin - read MLS o3lev data from bufr files. -! keep's Kris's ascii reading and rename the type to 'o3levtext'. -! 2012-12-14 jin/Wargan - fixed mflg error under "o3levtext" handling block. -! 2013-01-18 Wargan - read omieffnc from a NetCDF file, and removed omieff text. -! 2013-02-05 guo - STOP statements were replaced with die() to signal an _abort_. -! - dfile_format() is used to handle formats of "omieff" and "o3lev". -! - code is slightly cleaned up for possible restructuring -! with alternative format handling. -! 2013-09-27 guo - added this document block. -! 2014-02-03 guo - restructured from GMAO extended read_ozone() routine -! as a separated module for extended ozone obstypes, in -! particular, if not in BUFR format. -! - removed "text" option of "o3lev", for it is not been -! used anymore. -! - Moved history log messages here from read_ozone. -! 2015-09-17 Thomas - add l4densvar and thin4d to data selection procedure -! 2015-10-01 guo - consolidate use of ob location (in deg) -! 2016-09-19 guo - moved function dfile_format() here from obsmod.F90. -! 2018-05-21 j.jin - added time-thinning. Moved the checking of thin4d into satthin.F90. -! 2018-05-25 wargan - added OMPS and hooks for LIMS, UARS MLS, MIPAS -! -! 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 kinds, only: i_kind,r_kind,r_double - use obsmod, only: time_window_max - implicit none - private ! except - - public:: is_extOzone - public:: extOzone_read - -!! public:: extOzone_setupoz -!! public:: extOzone_setupozlev -!! public:: extOzone_setupoztot - -!! public:: extOzone_intoz -!! public:: extOzone_intozlev -!! public:: extOzone_intoztot - - interface is_extOzone; module procedure is_extOzone_; end interface - interface extOzone_read; module procedure read_; end interface - - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - character(len=*) , parameter:: myname='m_extOzone' - - integer(kind=i_kind), parameter:: MAXOBS_ = 1000000 - - real (kind=r_kind), parameter:: RMISS = -9999.9_r_kind - real (kind=r_kind), parameter:: badoz = 10000.0_r_kind - real (kind=r_kind), parameter:: r6 = 6.0_r_kind - real (kind=r_kind), parameter:: r360 = 360.0_r_kind - real (kind=r_kind) :: ptime,timeinflat,crit0 - integer(kind=i_kind) :: ithin_time,n_tbin,it_mesh - -contains -function is_extOzone_(dfile,dtype,dplat,class) - - use mpeu_util, only: die,perr - implicit none - logical:: is_extozone_ ! this is the function return variable - character(len=*),intent(in):: dfile ! observation input filename (see gsiparm.anl:&OBSINPUT/) - character(len=*),intent(in):: dtype ! observation type (see gsiparm.anl:&OBSINPUT/) - character(len=*),intent(in):: dplat ! platform (see gsiparm.anl:&OBSINPUT/) - character(len=*),optional,intent(in):: class ! specify either "level" or "layer" for sub-class - -! Use Case 1: (in read_obs()) -! -! elseif(is_extOzone(dfile,obstype,dplat)) then -! call read_ozone(obstype,dplat,dsis,dfile,...) -! -! -! Use Case 2: (in setuprhsall()) -! -! elseif(is_extOzone(dfile,obstype,dplat,class='level')) then -! call setupozlev(obstype,dplat,dsis,...) -! elseif(is_extOzone(dfile,obstype,dplat,class='layer')) then -! call setupozlay(obstype,dplat,dsis,...) -! elseif(is_extOzone(dfile,obstype,dplat,class='total')) then -! call setupozlay(obstype,dplat,dsis,...) -! -! Use Case 3: (where intoz() or stpoz() are called) -! -! elseif(is_extOzone(dfile,obstype,dplat)) then -! call intoz(obstype,dplat,dsis,...) - - character(len=*), parameter:: myname_=myname//'::is_extOzone_' - - integer(kind=i_kind),parameter:: iANY = 0 - integer(kind=i_kind),parameter:: iUNKNOWN=-1 - - integer(kind=i_kind),parameter:: iLEVEL = 1 - integer(kind=i_kind),parameter:: iLAYER = 2 - integer(kind=i_kind),parameter:: iTOTAL = 3 - - integer(kind=i_kind),parameter:: iTEXT = 1 - integer(kind=i_kind),parameter:: iBUFR = 2 - integer(kind=i_kind),parameter:: iNC = 3 - - integer(kind=i_kind):: class_,ifile_ - - ifile_=iUNKNOWN - select case(dfile_format(dfile)) - case("text") - ifile_=iTEXT - case("bufr") - ifile_=iBUFR - case("nc") - ifile_=iNC - end select - - class_=iANY - if(present(class)) then - select case(class) - case('level') - class_=iLEVEL - case('layer') - class_=iLAYER - case('total') - class_=iTOTAL - case default - class_=iUNKNOWN - call perr(myname_,'unknown ozone class, class =',class) - call perr(myname_,' dfile =',dfile) - call perr(myname_,' dtype =',dtype) - call perr(myname_,' dplat =',dplat) - call die(myname_) - end select - endif - - is_extOzone_= .false. -#ifndef NO_EXTRA_ - select case(class_) - case(iANY) - is_extOzone_= & - ifile_==iBUFR .and. dtype == 'o3lev' .or. & - ifile_==iNC .and. dtype == 'mls55' .or. & - ifile_==iNC .and. dtype == 'ompslpvis' .or. & - ifile_==iNC .and. dtype == 'ompslpuv' .or. & - ifile_==iNC .and. dtype == 'ompslp' .or. & - ifile_==iNC .and. dtype == 'lims' .or. & - ifile_==iNC .and. dtype == 'uarsmls' .or. & - ifile_==iNC .and. dtype == 'mipas' .or. & - ifile_==iNC .and. dtype == 'omieff' .or. & - ifile_==iNC .and. dtype == 'tomseff' - - case(iLEVEL) - is_extOzone_= & - ifile_==iBUFR .and. dtype == 'o3lev' .or. & - ifile_==iNC .and. dtype == 'mls55' .or. & - ifile_==iNC .and. dtype == 'ompslpvis' .or. & - ifile_==iNC .and. dtype == 'ompslpuv' .or. & - ifile_==iNC .and. dtype == 'ompslp' .or. & - ifile_==iNC .and. dtype == 'lims' .or. & - ifile_==iNC .and. dtype == 'uarsmls' .or. & - ifile_==iNC .and. dtype == 'mipas' - - case(iLAYER) - is_extOzone_= .false. - - case(iTOTAL) - is_extOzone_= & - ifile_==iNC .and. dtype == 'omieff' .or. & - ifile_==iNC .and. dtype == 'tomseff' - - case default - is_extOzone_= .false. - end select -#endif - -end function is_extOzone_ - -! ---------------------------------------------------------------------- -function dfile_format(dfile) result(dform) -!$$$ subprogram documentation block -! . . . . -! subprogram: function dfile_format -! prgmmr: j guo -! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.1 -! date: 2013-02-04 -! -! abstract: - check filename suffix to guess its format -! -! program history log: -! 2013-02-04 j guo - added this document block -! -! 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 - -! function interface: - - implicit none - - character(len=len('unknown')):: dform ! a 2-4 byte code for a format guess, - ! from a list of filename suffixes, 'bufr', 'text' (also for 'txt', - ! 'tcvitle', or 'vitl'), 'nc', or return a default value 'unknown'. One - ! can extend the list to other suffixes, such as 'hdf', 'hdf4', 'hdf5', - ! etc., if they are needed in the future. - character(len=*),intent(in):: dfile ! a given filename - -!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - character(len=*),parameter :: myname_=myname//'::dfile_format' - integer(i_kind):: i,l - - dform='unknown' - l=len_trim(dfile) - - i=max(0,l-6)+1 ! 6 byte code? - select case(dfile(i:l)) - case ('tcvitl') - dform='text' - end select - if(dform/='unknown') return - - i=max(0,l-4)+1! 4 byte code? - select case(dfile(i:l)) - case ('bufr','avhr') - dform='bufr' - case ('text','vitl') - dform='text' - end select - if(dform/='unknown') return - - i=max(0,l-3)+1 ! 3 byte code? - select case(dfile(i:l)) - case ('txt') ! a short - dform='text' - end select - if(dform/='unknown') return - - i=max(0,l-2)+1 ! 2 byte code? - select case(dfile(i:l)) - case ('nc') - dform='nc' - end select - if(dform/='unknown') return - - dform='bufr' -return -end function dfile_format - -subroutine read_(dfile,dtype,dplat,dsis, & ! intent(in), keys for type managing - nread,npuse,nouse, & ! intent(out), beside other implicit output variables - jsatid,gstime,lunout,twind,ithin,rmesh,nobs) ! intent(in), beside other implicit input variables - - use constants, only: zero - use mpeu_util, only: die,perr,tell - use mpimod, only: npe - use satthin, only: satthin_time_info => radthin_time_info -! use mpeu_util, only: mprefix,stdout -! nobs - array of observations on each subdomain for each processor - - implicit none - character(len=*), intent(in):: dfile ! obs_input filename - character(len=*), intent(in):: dtype ! observation type (see gsiparm.anl:&OBSINPUT/) - character(len=*), intent(in):: dplat ! platform (see gsiparm.anl:&OBSINPUT/) - character(len=*), intent(in):: dsis ! sensor/instrument/satellite tag (see gsiparm.anl:&OBSINPUT/ and ozinfo) - -! Use Case: (in read_ozone()) -! -! elseif(is_extOzone(dfile,dtype,dplat)) then -! call extOzone_read(dtype,dplat,dsis,dfile,...) -! - - character(len=*), parameter:: myname_=myname//'::read_' - - integer(kind=i_kind), intent(out):: nread ! number of obs record reads in this call - integer(kind=i_kind),dimension(npe), intent(inout):: nobs ! number of obs record reads in this call - integer(kind=i_kind), intent(out):: npuse ! nnmber of "preofiles" retained in this call - integer(kind=i_kind), intent(out):: nouse ! nnmber of obs retained in this call - - ! jsatid,gstime,lunout,twind,ithin,rmesh - character(len=*) , intent(in ):: jsatid ! platform ID (verification) - real (kind=r_kind), intent(in ):: gstime ! analysis time (minute) from reference date - integer(kind=i_kind), intent(in ):: lunout ! logical unit to send obs output. - real (kind=r_kind), intent(in ):: twind ! input group time window (hour) - - integer(kind=i_kind), intent(in ):: ithin ! flag to thin data - real (kind=r_kind), intent(in ):: rmesh ! thining mesh size (km) - - real(kind=r_kind),pointer,dimension(:,:):: p_out - - integer(kind=i_kind):: nreal,nchan,ilat,ilon - integer(kind=i_kind):: maxobs - integer(kind=i_kind):: i,k - - nread=-1 - npuse=-1 - nouse=-1 - - if(.not.is_extOzone(dfile,dtype,dplat)) then - call perr(myname_,'unexpected use, dfile =',dfile) - call perr(myname_,' dtype =',dtype) - call perr(myname_,' dplat =',dplat) - call die (myname_) - endif - call satthin_time_info(dtype, dplat, dsis, ptime, ithin_time) - if( ptime > 0.0_r_kind) then - n_tbin=nint(2*time_window_max/ptime) - else - n_tbin=1 - endif - - select case(dtype) - case('omieff','tomseff') ! layer-ozone or total-ozone types - select case(dfile_format(dfile)) - case('nc') - call oztot_ncInquire_(nreal,nchan,ilat,ilon, & - ithin,rmesh,maxobs) - - allocate(p_out(nreal+nchan,maxobs)) - p_out(:,:)=RMISS - - call oztot_ncRead_(dfile,dtype,dsis, p_out,nread,npuse,nouse, & - gstime, twind, ithin) - - ! Skip all "thinned" data records, and reset values of - ! nread, npuse, and nouse, as they are required by - ! upper level read_obs(). - - k=0 - do i=1,maxobs - if(p_out(1,i)>zero) then - k=k+1 - if(i>k) p_out(:,k)=p_out(:,i) - endif - enddo - nouse=k - npuse=k - end select - - case('o3lev') ! level-ozone types - select case(dfile_format(dfile)) -! case('text') -! call ozlev_textInquire(dfile,dtype,dplat, & -! nreal,nchan,ilat,ilon, maxobs) -! -! allocate(p_out(nreal+nchan,maxobs)) -! p_out(:,:)=RMISS -! -! call ozlev_textRead_(dfile,dtype,dplat,dsis, p_out,nread,npuse,nouse, & -! gstime,twind, nreal,nchan,ilat,ilon) -! - case('bufr') - call ozlev_bufrInquire_(nreal,nchan,ilat,ilon,maxobs) - - allocate(p_out(nreal+nchan,maxobs)) - p_out(:,:)=RMISS - - call ozlev_bufrRead_(dfile,dtype,dsis, p_out,nread,npuse,nouse, & - jsatid, gstime,twind) - end select - - case('mls55','ompslpvis','ompslpuv','ompslp','lims','uarsmls','mipas') - select case(dfile_format(dfile)) - case('nc') - call ozlev_ncInquire_( nreal,nchan,ilat,ilon,maxobs) - - allocate(p_out(nreal+nchan,maxobs)) - p_out(:,:)=RMISS - - call ozlev_ncRead_(dfile,dtype, p_out,nread,npuse,nouse, gstime,twind) - - end select - - end select - - if(nouse<0 .or. .not.associated(p_out)) then - call perr(myname_,'can not process, dtype =',trim(dtype)) - call perr(myname_,' dplat =',trim(dplat)) - call perr(myname_,' dsis =',trim(dsis)) - call perr(myname_,' dfile =',trim(dfile)) - call perr(myname_,' associated(p_out) =',associated(p_out)) - if(associated(p_out)) then - call perr(myname_,' actual size(p_out,1) =',size(p_out,1)) - call perr(myname_,' actual size(p_out,2) =',size(p_out,2)) - call perr(myname_,' nread =',nread) - call perr(myname_,' npuse =',npuse) - call perr(myname_,' nouse =',nouse) - endif - call die(myname_) - endif - - if(nreal+nchan/=size(p_out,1)) then - call perr(myname_,'unexpected size(p_out,1), nreal+nchan =',nreal+nchan) - call perr(myname_,' nreal =',nreal) - call perr(myname_,' nchan =',nchan) - call perr(myname_,' actual size(p_out,1) =',size(p_out,1)) - call perr(myname_,' actual size(p_out,2) =',size(p_out,2)) - call die(myname_) - endif - -! Output candidate observations from _dfile_ to _obsfile_ (pre-opened in lunout). - -! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & -! ' obstype,nread,npuse,nouse,no/npuse = ',dtype,nread,npuse,nouse,real(nouse)/npuse - - ! While nreal+nchan defines the leading dimension of p_out(:,:), it is - ! obviously a bad idea that n_out has been missing from the header for - ! its second dimension. - - call count_obs(npuse,nreal+nchan,ilat,ilon,p_out,nobs) - write(lunout) dtype, dsis, nreal, nchan, ilat, ilon - write(lunout) p_out(:,1:npuse) - - deallocate(p_out) - -end subroutine read_ - -subroutine oztot_ncInquire_( nreal,nchan,ilat,ilon, ithin,rmesh,maxrec) - use satthin, only: satthin_itxmax => itxmax - use satthin, only: satthin_makegrids => makegrids - implicit none - - integer(kind=i_kind), intent(out):: nreal ! number of real parameters per record - integer(kind=i_kind), intent(out):: nchan ! number of channels or levels per record - integer(kind=i_kind), intent(out):: ilat ! index to latitude in nreal parameters. - integer(kind=i_kind), intent(out):: ilon ! index to longitude in nreal parameters. - - integer(kind=i_kind), intent(in ):: ithin ! flag to thin data - real (kind=r_kind), intent(in ):: rmesh ! size (km) of the thinning mesh - - integer(kind=i_kind), intent(out):: maxrec ! extimated input record count - - character(len=*), parameter:: myname_=myname//'::oztot_ncInquire_' - -! Configure the record buffer for this obs class - nreal=36 - nchan=1 - ilat=4 - ilon=3 - -! Make thinning grids, and to define the total record size. - call satthin_makegrids(rmesh,ithin,n_tbin=n_tbin) - maxrec=satthin_itxmax - -end subroutine oztot_ncInquire_ - -!.................................................................................. -subroutine oztot_ncread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, & - gstime,twind, ithin) -!.................................................................................. - - use netcdf, only: nf90_open - use netcdf, only: nf90_nowrite - use netcdf, only: nf90_noerr - use netcdf, only: nf90_inq_dimid - use netcdf, only: nf90_inquire_dimension - use netcdf, only: nf90_inq_varid - use netcdf, only: nf90_get_var - use netcdf, only: nf90_close - - use satthin, only: satthin_makegrids => makegrids - use satthin, only: satthin_tdiff2crit => tdiff2crit - use satthin, only: satthin_map2tgrid => map2tgrid - use satthin, only: satthin_finalcheck => finalcheck - use satthin, only: satthin_destroygrids => destroygrids - - use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons - use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar - use obsmod, only: nloz_omi - - use constants, only: deg2rad,zero,r60inv -! use mpeu_util, only: mprefix,stdout - - implicit none - character(len=*), intent(in):: dfile ! obs_input filename - character(len=*), intent(in):: dtype ! observation type (see gsiparm.anl:&OBSINPUT/) - character(len=*), intent(in):: dsis ! sensor/instrument/satellite tag (see gsiparm.anl:&OBSINPUT/ and ozinfo) - - real (kind=r_kind), dimension(:,:), intent(out):: ozout - integer(kind=i_kind), intent(out):: nmrecs ! read count in ozout - integer(kind=i_kind), intent(out):: ndata ! "good" profile count in ozout - integer(kind=i_kind), intent(out):: nodata ! "good" data count in ozout - - real (kind=r_kind), intent(in):: gstime ! analysis time (minute) from reference date - real (kind=r_kind), intent(in):: twind ! input group time window (hour) - - - integer(kind=i_kind), intent(in ):: ithin ! flag to thin data - - character(len=*), parameter:: myname_=myname//'::oztot_ncRead_' - -! parameters for output bookkeeping - integer(kind=i_kind):: i, irec - integer(kind=i_kind):: itx,itt - -! parameters for NetCDF arrays - integer(kind=i_kind), allocatable :: iya(:),ima(:),idda(:),ihha(:),imina(:),fovna(:) - integer(kind=i_kind), allocatable :: toqfa(:),alqfa(:) - real (kind=r_kind), allocatable :: rseca(:),slatsa(:),slonsa(:),totoza(:),szaa(:) - real (kind=r_kind), allocatable :: aprioria(:,:),efficiencya(:,:) - integer(kind=i_kind) nrecDimId,nomilevsDimID,lonVarID,latVarID,yyVarID,mmVarID - integer(kind=i_kind) ddVarID,hhVarID,minVarID,ssVarID,fovnVarID,toqfVarID,alqfVarID - integer(kind=i_kind) szaVarID,totozVarID,aprioriVarID,efficiencyVarID - integer(kind=i_kind) ier, ncid, nomilevs - integer(kind=i_kind):: iy,im,idd,ihh,imin,nmind - integer(kind=i_kind):: toqf,alqf,fovn - real (kind=r_kind):: rsec,slats,slons - real (kind=r_kind),dimension(nloz_omi):: apriori, efficiency - integer(kind=i_kind):: binary(17) - - real (kind=r_kind):: dlon,dlon_earth,dlon_earth_deg - real (kind=r_kind):: dlat,dlat_earth,dlat_earth_deg - real (kind=r_kind):: tdiff,sstime,t4dv,crit1,dist1,rsat - integer(kind=i_kind):: idate5(5) - - real (kind=r_double):: totoz, sza - - logical:: outside, removeScans, iuse - integer(kind=i_kind):: maxobs - - - removeScans = .true. ! remove OMI scan numbers >= 25? - maxobs=size(ozout,2) - rsat=999._r_kind - -! Using OMI/TOMS with efficience factors in NetCDF format - nodata = 0 - ndata = 0 - - call check(nf90_open(trim(dfile),nf90_nowrite,ncid),stat=ier) - - ! ignore if the file is not actually present. - if(ier/=nf90_noerr) then - call satthin_destroygrids() - return - end if - - ! Get dimensions from OMI input file - call check(nf90_inq_dimid(ncid, "nrec", nrecDimId),stat=ier) - - ! ignore if the file header is empty - if(ier/=nf90_noerr) then - call check(nf90_close(ncid),stat=ier) - call satthin_destroygrids() - return - endif - - ! Get dimensions from OMI/TOMS input file -!!!! nmrecs=0 - call check(nf90_inquire_dimension(ncid, nrecDimId, len = nmrecs),stat=ier) - - ! ignore if the file header is empty - if(ier/=nf90_noerr .or. nmrecs==0) then - call check(nf90_close(ncid),stat=ier) - call satthin_destroygrids() - return - endif - - ! Continue the input - call check(nf90_inq_dimid(ncid, "nlevs", nomilevsDimId)) - call check(nf90_inquire_dimension(ncid, nomilevsDimId, len = nomilevs)) - - ! We have dimensions so we can allocate arrays - allocate(iya(nmrecs),ima(nmrecs),idda(nmrecs),ihha(nmrecs),imina(nmrecs), & - rseca(nmrecs),fovna(nmrecs),slatsa(nmrecs),slonsa(nmrecs),totoza(nmrecs), & - toqfa(nmrecs),alqfa(nmrecs),szaa(nmrecs)) - allocate(aprioria(nomilevs,nmrecs),efficiencya(nomilevs,nmrecs)) - - ! Read variables and store them in these arrays - call check(nf90_inq_varid(ncid, "lon", lonVarId)) - call check(nf90_get_var(ncid, lonVarId, slonsa)) - - call check(nf90_inq_varid(ncid, "lat", latVarId)) - call check(nf90_get_var(ncid, latVarId, slatsa)) - - call check(nf90_inq_varid(ncid, "yy", yyVarId)) - call check(nf90_get_var(ncid, yyVarId, iya)) - - call check(nf90_inq_varid(ncid, "mm", mmVarId)) - call check(nf90_get_var(ncid, mmVarId, ima)) - - call check(nf90_inq_varid(ncid, "dd", ddVarId)) - call check(nf90_get_var(ncid, ddVarId, idda)) - - call check(nf90_inq_varid(ncid, "hh", hhVarId)) - call check(nf90_get_var(ncid, hhVarId, ihha)) - - call check(nf90_inq_varid(ncid, "min", minVarId)) - call check(nf90_get_var(ncid, minVarId, imina)) - - call check(nf90_inq_varid(ncid, "ss", ssVarId)) - call check(nf90_get_var(ncid, ssVarId, rseca)) - - call check(nf90_inq_varid(ncid, "fov", fovnVarId)) - call check(nf90_get_var(ncid, fovnVarId, fovna)) - - call check(nf90_inq_varid(ncid, "sza", szaVarId)) - call check(nf90_get_var(ncid, szaVarId, szaa)) - - call check(nf90_inq_varid(ncid, "toqf", toqfVarId)) - call check(nf90_get_var(ncid, toqfVarId, toqfa)) - - call check(nf90_inq_varid(ncid, "alqf", alqfVarId)) - call check(nf90_get_var(ncid, alqfVarId, alqfa)) - - call check(nf90_inq_varid(ncid, "toz", totozVarId)) - call check(nf90_get_var(ncid, totozVarId, totoza)) - - call check(nf90_inq_varid(ncid, "apriori", aprioriVarId)) - call check(nf90_get_var(ncid, aprioriVarId, aprioria)) - - call check(nf90_inq_varid(ncid, "efficiency", efficiencyVarId)) - call check(nf90_get_var(ncid, efficiencyVarId, efficiencya)) - - ! close the data file - call check(nf90_close(ncid)) - - ! now screen the data and put them into the right places - recloop: do irec = 1, nmrecs - iy = iya(irec) - im = ima(irec) - idd = idda(irec) - ihh = ihha(irec) - imin = imina(irec) - rsec = rseca(irec) - fovn = fovna(irec) - slats = slatsa(irec) - slons = slonsa(irec) - totoz = totoza(irec) - toqf = toqfa(irec) - alqf = alqfa(irec) - sza = szaa(irec) - do i = 1, nomilevs - apriori(i) = aprioria(i, irec) - ! Reported efficiencies from layer 4 up (counting from the surface - ! so from 1 - 8 here) are incorrect (too large) for high SZA - ! Setting them all to 1.0 per PK Bhartia''s recommendation - if (i .le. 8) then - efficiency(i) = 1._r_kind - else - efficiency(i) = efficiencya(i, irec) - endif - end do - -!!!! nmrecs=nmrecs+1 - ! Convert observation location to radians - if(abs(slats)>90._r_kind .or. abs(slons)>r360) cycle recloop - if(slons< zero) slons=slons+r360 - if(slons==r360) slons=zero - dlat_earth_deg = slats - dlon_earth_deg = slons - dlat_earth = slats * deg2rad - dlon_earth = slons * deg2rad - - if(regional)then - call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) - if(outside) cycle recloop - else - dlat = dlat_earth - dlon = dlon_earth - call grdcrd1(dlat,rlats,nlat,1) - call grdcrd1(dlon,rlons,nlon,1) - endif - -! convert observation time to relative time - idate5(1) = iy !year - idate5(2) = im !month - idate5(3) = idd !day - idate5(4) = ihh !hour - idate5(5) = imin !minute - call w3fs21(idate5,nmind) - - t4dv=real((nmind-iwinbgn),r_kind)*r60inv - sstime=real(nmind,r_kind) - tdiff=(sstime-gstime)*r60inv - if (l4dvar.or.l4densvar) then - if (t4dvwinlen) cycle recloop - else - if(abs(tdiff) > twind) cycle recloop - end if - - if (totoz > badoz ) cycle recloop - -! Apply data screening based on quality flags -! Bit 10 (from the left) in TOQF represents row anomaly. All 17 bits in toqf is -! supposed to converted into array elements of binary(:), either for "tomseff" or -! "omieff". - binary(:)=0 - select case(dtype) - case('omieff') - - if (toqf .ne. 0 .and. toqf .ne. 1) cycle recloop - -! Remove obs at high solar zenith angles - if (sza > 84.0_r_kind) cycle recloop - -! remove the bad scan position data: fovn beyond 25 - if (removeScans) then - if (fovn >=25_i_kind) cycle recloop - endif - if (fovn <=2_i_kind .or. fovn >=58_i_kind) cycle recloop - -! remove the data in which the C-pair algorithm ((331 and 360 nm) is used. - if (alqf == 3_i_kind .or. alqf == 13_i_kind) cycle recloop - - case('tomseff') - ! The meaning of quality flags for TOMS version 8 is similar to that - ! for SBUV: - ! 0 - good data, 1 - good data with SZA > 84 deg - if (toqf /= 0) cycle recloop - - case default - end select - -! thin OMI/TOMS data - - crit0 = 0.01_r_kind - timeinflat=r6 - call satthin_tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh) - if (ithin /= 0) then - call satthin_map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,dsis,it_mesh=it_mesh) - if(.not. iuse)cycle recloop - call satthin_finalcheck(dist1,crit1,itx,iuse) - if(.not. iuse)cycle recloop - ndata=ndata+1 - nodata=ndata - else - ndata=ndata+1 - nodata=ndata - itx= ndata - end if - -!! ASSERT_(size(ozout,2)>=itx) - - if(itx <= maxobs) then - ozout(1,itx)=rsat - ozout(2,itx)=t4dv - ozout(3,itx)=dlon ! grid relative longitude - ozout(4,itx)=dlat ! grid relative latitude - ozout(5,itx)=dlon_earth_deg ! earth relative longitude (degrees) - ozout(6,itx)=dlat_earth_deg ! earth relative latitude (degrees) - ozout(7,itx)=real(toqf) ! - total ozone quality code (not used) - ozout(8,itx)=real(sza) ! solar zenith angle - ozout(9,itx)=binary(10) ! row anomaly flag, is actually fixed to 0 - ozout(10,itx)=0. ! - cloud amount (not used) - ozout(11,itx)=0. ! - vzan (not used) - ozout(12,itx)=0. ! - aerosol index (not used) - ozout(13,itx)=0. ! - ascending/descending (not used) - ozout(14,itx)=real(fovn) ! scan position - ! "(not used)" flags above imply that they - ! are not used in setupozlay(). - -! Added apriori and efficiency profiles - ozout(15:25,itx)=apriori - ozout(26:36,itx)=efficiency - ozout(37,itx)=totoz - endif - -!! ASSERT_(size(ozout,1)==37) - - end do recloop - - deallocate(iya,ima,idda,ihha,imina, & - rseca,fovna,slatsa,slonsa,totoza, & - toqfa,alqfa,szaa,aprioria,efficiencya) -! end -! End of loop over observations -! End of OMI block with efficiency factors in NetCDF format - - call satthin_destroygrids() - - nodata=min(ndata,maxobs) -! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & -! ' obstype,nmrecs,ndata,nodata,no/ndata = ',dtype,nmrecs,ndata,nodata,real(nodata)/ndata - return - -end subroutine oztot_ncread_ -!.................................................................................. - - -subroutine ozlev_ncInquire_( nreal,nchan,ilat,ilon, maxrec) - implicit none - - integer(kind=i_kind), intent(out):: nreal ! number of real parameters per record - integer(kind=i_kind), intent(out):: nchan ! number of channels or levels per record - integer(kind=i_kind), intent(out):: ilat ! index to latitude in nreal parameters. - integer(kind=i_kind), intent(out):: ilon ! index to longitude in nreal parameters. - - integer(kind=i_kind), intent(out):: maxrec ! extimated input record count - - character(len=*), parameter:: myname_=myname//'::ozlev_ncInquire_' - -! Configure the record, they are not (dfile,dtype,dplat) dependent in this case. - nreal = 12 - nchan = 1 ! There are 'levs' levels but each is treated - ! as a separate observation so that nchanl = 1 - ilat=4 - ilon=3 - - maxrec = MAXOBS_ -end subroutine ozlev_ncInquire_ - -!.................................................................................. -subroutine ozlev_ncread_(dfile,dtype,ozout,nmrecs,ndata,nodata, gstime,twind) -!.................................................................................. - use netcdf, only: nf90_open - use netcdf, only: nf90_nowrite - use netcdf, only: nf90_noerr - use netcdf, only: nf90_inq_dimid - use netcdf, only: nf90_inquire_dimension - use netcdf, only: nf90_inq_varid - use netcdf, only: nf90_get_var - use netcdf, only: nf90_close - - use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons - use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar - - use constants, only: deg2rad,zero,one_tenth,r60inv - use ozinfo, only: jpch_oz,nusis_oz,iuse_oz - use mpeu_util, only: perr,die -! use mpeu_util, only: mprefix,stdout - - implicit none - character(len=*), intent(in):: dfile ! obs_input filename - character(len=*), intent(in):: dtype ! obs_input dtype - - real (kind=r_kind), dimension(:,:), intent(out):: ozout - integer(kind=i_kind), intent(out):: nmrecs ! count of records read - integer(kind=i_kind), intent(out):: ndata ! count of processed - integer(kind=i_kind), intent(out):: nodata ! count of retained - - real (kind=r_kind), intent(in):: gstime ! analysis time (minute) from reference date - real (kind=r_kind), intent(in):: twind ! input group time window (hour) - - - character(len=*), parameter:: myname_=myname//'::ozlev_ncRead_' - - integer(kind=i_kind):: ier, iprof, nprofs, maxobs - integer(kind=i_kind):: i, ilev, ikx, ncid, k0 - integer(kind=i_kind),allocatable,dimension(:):: ipos - - integer(kind=i_kind):: nrecDimId,lonVarID,latVarID,yyVarID,mmVarID - integer(kind=i_kind):: ddVarID,hhVarID,minVarID,ssVarID - integer(kind=i_kind):: pressVarID - integer(kind=i_kind):: convVarID, errVarID, ozoneVarID - integer(kind=i_kind):: levsDimID,levs - - integer(kind=i_kind), allocatable :: iya(:),ima(:),idda(:),ihha(:),imina(:),iseca(:) - real (kind=r_kind), allocatable :: slatsa(:),slonsa(:) - real (kind=r_kind), allocatable :: press(:), ozone(:,:), qual(:), press2d(:,:) - real (kind=r_kind), allocatable :: conv(:), err(:,:) - - integer(kind=i_kind):: nmind - - real (kind=r_kind):: slons0,slats0 - real (kind=r_kind):: ppmv, pres, pob, obserr, usage - - real (kind=r_kind):: dlon,dlon_earth,dlon_earth_deg - real (kind=r_kind):: dlat,dlat_earth,dlat_earth_deg - real (kind=r_kind):: tdiff,sstime,t4dv,rsat - integer(kind=i_kind):: idate5(5) - - logical:: outside - logical:: first - - maxobs=size(ozout,2) - rsat=999._r_kind - ndata = 0 - nmrecs=0 - nodata=-1 -!.................................................................................. - ! ---------------MLS and other limb ozone data in NetCDF format ---------- - ! Open file and read dimensions - call check(nf90_open(trim(dfile),nf90_nowrite,ncid),stat=ier) - - ! ignore if the file is not actually present. - if(ier/=nf90_noerr) return - - ! Get dimensions from OMI input file - call check(nf90_inq_dimid(ncid, "nprofiles", nrecDimId),stat=ier) - - ! ignore if the file header is empty - if(ier/=nf90_noerr) then - call check(nf90_close(ncid),stat=ier) - return - endif - - ! Get dimensions from the input file: # of profiles and # of levels - nprofs=0 - call check(nf90_inquire_dimension(ncid, nrecDimId, len = nprofs),stat=ier) - - ! ignore if the file header is empty - if(ier/=nf90_noerr) then - call check(nf90_close(ncid),stat=ier) - return - endif - - if(nprofs==0) then - nodata=0 - call check(nf90_close(ncid),stat=ier) - return - endif - - ! Continue the input - call check(nf90_inq_dimid(ncid, "nlevs", levsDimId)) - call check(nf90_inquire_dimension(ncid, levsDimId, len = levs)) - - ! NOTE: Make sure that 'ozinfo' has the same number of levels - ! for NRT it is 55 - allocate(ipos(levs)) - ipos=999 - - ! Process limb data in NetDCF format - ikx = 0 - first=.false. - do i=1,jpch_oz - if( (.not. first) .and. index(nusis_oz(i), trim(dtype))/=0) then - k0=i - first=.true. - end if - if(first.and.index(nusis_oz(i),trim(dtype))/=0) then - ikx=ikx+1 - ipos(ikx)=k0+ikx-1 - end if - end do - - if (ikx/=levs) call die(myname_//': inconsistent levs for '//dtype) - - nmrecs=0 - ! Allocate space and read data - - allocate(iya(nprofs),ima(nprofs),idda(nprofs),ihha(nprofs),imina(nprofs), & - iseca(nprofs),slatsa(nprofs),slonsa(nprofs), ozone(levs,nprofs), & - err(levs,nprofs),qual(nprofs), conv(nprofs)) - if (index(dtype, 'ompslp') == 0) then - allocate(press(levs)) - else - allocate(press2d(levs,nprofs)) - endif - - ! Read variables and store them in these arrays - call check(nf90_inq_varid(ncid, "lon", lonVarId)) - call check(nf90_get_var(ncid, lonVarId, slonsa)) - - call check(nf90_inq_varid(ncid, "lat", latVarId)) - call check(nf90_get_var(ncid, latVarId, slatsa)) - - call check(nf90_inq_varid(ncid, "yy", yyVarId)) - call check(nf90_get_var(ncid, yyVarId, iya)) - - call check(nf90_inq_varid(ncid, "mm", mmVarId)) - call check(nf90_get_var(ncid, mmVarId, ima)) - - call check(nf90_inq_varid(ncid, "dd", ddVarId)) - call check(nf90_get_var(ncid, ddVarId, idda)) - - call check(nf90_inq_varid(ncid, "hh", hhVarId)) - call check(nf90_get_var(ncid, hhVarId, ihha)) - - call check(nf90_inq_varid(ncid, "min", minVarId)) - call check(nf90_get_var(ncid, minVarId, imina)) - - call check(nf90_inq_varid(ncid, "ss", ssVarId)) - call check(nf90_get_var(ncid, ssVarId, iseca)) - - if (index(dtype, 'ompslp') == 0) then - call check(nf90_inq_varid(ncid, "press", pressVarId)) - call check(nf90_get_var(ncid, pressVarId, press)) - else - call check(nf90_inq_varid(ncid, "press", pressVarId)) - call check(nf90_get_var(ncid, pressVarId, press2d)) - endif - - if (index(dtype, 'ompslp') == 0) then - call check(nf90_inq_varid(ncid, "conv", convVarId)) - call check(nf90_get_var(ncid, convVarId, conv)) - end if - - ! call check(nf90_inq_varid(ncid, "qual", qualVarId)) - ! call check(nf90_get_var(ncid, qualVarId, qual)) - - call check(nf90_inq_varid(ncid, "oberr", errVarId)) - call check(nf90_get_var(ncid, errVarId, err)) - - call check(nf90_inq_varid(ncid, "ozone", ozoneVarId)) - call check(nf90_get_var(ncid, ozoneVarId, ozone)) - - ! close the data file - call check(nf90_close(ncid)) - - ! 'Unpack' the data - nmrecs = 0 - nodata = 0 - do iprof = 1, nprofs - do ilev = 1, levs - ! note that most of the quality control is done at the - ! pre-processing stage -!! if (press(ilev) > 262.0_r_kind .or. press(ilev) < 0.1_r_kind ) cycle ! undefined - if (ozone(ilev, iprof) < -900.0_r_kind) cycle ! undefined - if (err(ilev, iprof) < -900.0_r_kind) cycle ! undefined - if (iuse_oz(ipos(ilev)) < 0) then - usage = 100._r_kind - else - usage = zero - endif - nmrecs=nmrecs+1 - - ! convert observation location to radians - slons0=slonsa(iprof) - slats0=slatsa(iprof) - if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle - if(slons0< zero) slons0=slons0+r360 - if(slons0==r360) slons0=zero - dlat_earth_deg = slats0 - dlon_earth_deg = slons0 - dlat_earth = slats0 * deg2rad - dlon_earth = slons0 * deg2rad - - if(regional)then - call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) - if(outside) cycle - else - dlat = dlat_earth - dlon = dlon_earth - call grdcrd1(dlat,rlats,nlat,1) - call grdcrd1(dlon,rlons,nlon,1) - endif - - idate5(1) = iya(iprof) !year - idate5(2) = ima(iprof) !month - idate5(3) = idda(iprof) !day - idate5(4) = ihha(iprof) !hour - idate5(5) = imina(iprof) !minute - call w3fs21(idate5,nmind) - t4dv=real((nmind-iwinbgn),r_kind)*r60inv - if (l4dvar.or.l4densvar) then - if (t4dvwinlen) then - write(6,*)'read_ozone: ', dtype,' obs time idate5=',idate5,', t4dv=',& - t4dv,' is outside time window, sstime=',sstime*r60inv - cycle - end if - else - sstime=real(nmind,r_kind) - tdiff=(sstime-gstime)*r60inv - if(abs(tdiff) > twind)then - write(6,*)'read_ozone: ',dtype,' obs time idate5=',idate5,', tdiff=',& - tdiff,' is outside time window=',twind - cycle - end if - end if - - obserr = err(ilev, iprof) - ppmv = ozone(ilev, iprof) - if (index(dtype, 'ompslp') == 0) then - pres = press(ilev) - else - pres = press2d(ilev, iprof) - end if - pob = log(pres * one_tenth) - ndata = ndata+1 - if(ndata<=maxobs) then - nodata = nodata + 1 - ozout(1,ndata)=rsat - ozout(2,ndata)=t4dv - ozout(3,ndata)=dlon ! grid relative longitude - ozout(4,ndata)=dlat ! grid relative latitude - ozout(5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) - ozout(6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) - - ozout(7,ndata)=rmiss ! used to be solar zenith angle - ozout(8,ndata)=usage - ozout(9,ndata)=pob ! pressure - ozout(10,ndata)=obserr ! ozone mixing ratio precision in ppmv - ozout(11,ndata)=float(ipos(ilev)) ! pointer of obs level index in ozinfo.txt - ozout(12,ndata)=levs ! # of vertical levels - ozout(13,ndata)=ppmv ! ozone mixing ratio in ppmv - endif - - end do - end do - - - deallocate(iya,ima,idda,ihha,imina,iseca,slatsa,slonsa, ozone, & - err,qual, conv) - if (index(dtype, 'ompslp') == 0) then - deallocate(press) - else - deallocate(press2d) - end if - deallocate(ipos) - - ! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & - ! ' obstype,nmrecs,ndata,nodata,no/ndata = ',dtype,nmrecs,ndata,nodata,real(nodata)/ndata - - if (ndata > maxobs) then - call perr(myname_,'Number of limb obs reached maxobs = ', maxobs) - call perr(myname_,' ndata = ', ndata) - call perr(myname_,' nodata = ', nodata) - call die(myname_) - endif - - !---------------END MLS NRT NetCDF--------------------------- -end subroutine ozlev_ncread_ - -subroutine ozlev_bufrInquire_(nreal,nchan,ilat,ilon,maxrec) - implicit none - - integer(kind=i_kind), intent(out):: nreal ! number of real parameters per record - integer(kind=i_kind), intent(out):: nchan ! number of channels or levels per record - integer(kind=i_kind), intent(out):: ilat ! index to latitude in nreal parameters. - integer(kind=i_kind), intent(out):: ilon ! index to longitude in nreal parameters. - - integer(kind=i_kind), intent(out):: maxrec ! extimated input record count - - character(len=*), parameter:: myname_=myname//'::ozlev_bufrInquire_' - -! Configure the record, they are not (dfile,dtype,dplat) dependent in this case. - nreal = 12 - nchan = 1 - ilat=4 - ilon=3 - - maxrec = MAXOBS_ -end subroutine ozlev_bufrInquire_ - -subroutine ozlev_bufrread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, & - jsatid, gstime,twind) - - use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons - use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar - - use constants, only: deg2rad,zero,r60inv - use ozinfo, only: jpch_oz,nusis_oz,iuse_oz - use radinfo, only: dec2bin - - use mpeu_util, only: warn,tell -! use mpeu_util, only: mprefix,stdout - - implicit none - character(len=*), intent(in):: dfile ! obs_input filename - character(len=*), intent(in):: dtype ! observation type (see gsiparm.anl:&OBSINPUT/) - character(len=*), intent(in):: dsis ! sensor/instrument/satellite tag (see gsiparm.anl:&OBSINPUT/ and ozinfo) - - real (kind=r_kind), dimension(:,:), intent(out):: ozout - integer(kind=i_kind), intent(out):: nmrecs ! count of actual records read - integer(kind=i_kind), intent(out):: ndata ! count of processed data - integer(kind=i_kind), intent(out):: nodata ! count of retained data - - character(len=*) , intent(in):: jsatid ! platform ID (verification) - real (kind=r_kind), intent(in):: gstime ! analysis time (minute) from reference date - real (kind=r_kind), intent(in):: twind ! input group time window (hour) - - - character(len=*),parameter:: myname_=myname//'::ozlev_bufrread_' - - integer(kind=i_kind),parameter:: nloz =37 - integer(kind=i_kind),parameter:: lunin=10 - - character(len=*), parameter:: mlstr ='SAID CLAT CLON YEAR MNTH DAYS HOUR MINU SECO SOZA CONV MLST PCCF' - character(len=*), parameter:: mlstrl2='PRLC OZMX OZMP OZME' - - integer(kind=i_kind):: idate,jdate,ksatid,iy,iret,im,ihh,idd - integer(kind=i_kind):: mpos - character(len= 8):: subset - - integer(kind=i_kind):: maxobs - integer(kind=i_kind):: nmind - integer(kind=i_kind):: k - integer(kind=i_kind):: kidsat - integer(kind=i_kind):: idate5(5) - integer(kind=i_kind):: ikx - integer(kind=i_kind):: decimal,binary_mls(18) - - real(kind=r_kind):: tdiff,sstime,dlon,dlat,t4dv - real(kind=r_kind):: slons0,slats0,rsat,dlat_earth,dlon_earth - real(kind=r_kind):: dlat_earth_deg,dlon_earth_deg - real(kind=r_double),dimension(13):: hdrmls - real(kind=r_double),dimension(4,37):: hdrmlsl2 - real(kind=r_double):: hdrmls13 - real(kind=r_kind),allocatable,dimension(:):: mlspres,mlsoz,mlsozpc,usage1 - integer(kind=i_kind),allocatable,dimension(:):: ipos - integer(kind=i_kind):: iprofs,nprofs - logical:: outside - - maxobs=size(ozout,2) - rsat=999._r_kind - - open(lunin,file=dfile,form='unformatted') - - call openbf(lunin,'IN',lunin) - call datelen(10) - call readmg(lunin,subset,idate,iret) - - ! If it has failed at the first read(), ... - if (iret/=0 .or. subset /= 'GM008015') then - call closbf(lunin) - close(lunin) - - call warn(myname_,'Failed at reading BUFR file, dfile =',trim(dfile)) - call warn(myname_,' dtype =',trim(dtype)) - call warn(myname_,' jsatid =',trim(jsatid)) - call warn(myname_,' lunin =',lunin) - call warn(myname_,' iret =',iret) - call warn(myname_,' subset =',trim(subset)) - - nprofs = 0 - nodata = 0 - - return - endif - - call tell(myname_,'MLS o3lev data type, subset =',trim(subset)) - write(6,*)'read_ozone: GMAO MLS o3lev data type, subset=',subset - - ! Q: o3lev data has 37 levels. Then, why is size(ipos) of 44? - j.jin - ! A: Because there are 44 entries in ozinfo.txt at the time - j.guo - - mpos=max(nloz,jpch_oz) - allocate (ipos(mpos)) ! 44? 37? - allocate (usage1(nloz)) - - ipos=999 - nmrecs=0 - nprofs=0 - nodata=0 - -! Set dependent variables and allocate arrays - - allocate (mlspres(nloz)) - allocate (mlsoz(nloz)) - allocate (mlsozpc(nloz)) - - ikx=0 - do k=1,jpch_oz - if(index(nusis_oz(k),'o3lev')/=0) then ! all "o3lev" in ozinfo.txt - ikx=ikx+1 - ipos(ikx)=k - end if - end do - - iy=0 - im=0 - idd=0 - ihh=0 - -!! This is the top of the profile loop - ndata=0 - iprofs=0 - obsloop: do - call readsb(lunin,iret) - if (iret/=0) then !JJJ, end of the subset - call readmg(lunin,subset,jdate,iret) !JJJ open a new mg - if (iret/=0) exit obsloop !JJJ, no more mg, EOF - cycle obsloop - endif - - do k=1,nloz - if (iuse_oz(ipos(k)) < 0) then - usage1(k) = 100._r_kind - else - usage1(k) = zero - endif - end do - -! extract header information - call ufbint(lunin,hdrmls,13,1,iret,mlstr) - rsat = hdrmls(1); ksatid=rsat - - if(jsatid == 'aura')kidsat = 785 - if (ksatid /= kidsat) cycle obsloop - - nmrecs=nmrecs+nloz - -! Convert observation location to radians - slats0= hdrmls(2) - slons0= hdrmls(3) - if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle obsloop - if(slons0< zero) slons0=slons0+r360 - if(slons0==r360) slons0=zero - dlat_earth_deg = slats0 - dlon_earth_deg = slons0 - dlat_earth = slats0 * deg2rad - dlon_earth = slons0 * deg2rad - - if(regional)then - call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) - if(outside) cycle obsloop - else - dlat = dlat_earth - dlon = dlon_earth - call grdcrd1(dlat,rlats,nlat,1) - call grdcrd1(dlon,rlons,nlon,1) - endif - -! convert observation time to relative time - idate5(1) = hdrmls(4) !year - idate5(2) = hdrmls(5) !month - idate5(3) = hdrmls(6) !day - idate5(4) = hdrmls(7) !hour - idate5(5) = hdrmls(8) !minute - call w3fs21(idate5,nmind) - - t4dv=real((nmind-iwinbgn),r_kind)*r60inv - if (l4dvar.or.l4densvar) then - if (t4dvwinlen) then - write(6,*)'read_ozone: mls obs time idate5=',idate5,', t4dv=',& - t4dv,' is outside time window, sstime=',sstime*r60inv - cycle obsloop - endif - else - sstime=real(nmind,r_kind) - tdiff=(sstime-gstime)*r60inv - if (abs(tdiff) > twind) then - write(6,*)'read_ozone: mls obs time idate5=',idate5,', tdiff=',& - tdiff,' is outside time window=',twind - cycle obsloop - endif - end if - -! v2.2 data screening, only accept: -! Pressure range: 215-0.02mb -! Precision: positive OZMP; -! Status flag: only use even number -! Quality(PCCF): use >1.2 for data at 215-100mb & low latitude, -! use >0.4 for data elsewhere -! Convergence: use <1.8 - -! Bit 1 in MLST represents data should not be used -! Note: in BUFR bits are defined from left to right as: 123456789... -! whereas in HDF5 (and the nasa document) bits are defined from right to left as: ...876543210 - - decimal=int(hdrmls(12)) - call dec2bin(decimal,binary_mls,18) - if (binary_mls(1) == 1 ) cycle obsloop - - if(hdrmls(11) >= 1.8_r_kind) cycle obsloop - -! extract pressure, ozone mixing ratio and precision - call ufbrep(lunin,hdrmlsl2,4,nloz,iret,mlstrl2) - - iprofs=iprofs+1 ! counting the profiles - do k=1,nloz - mlspres(k)=log(hdrmlsl2(1,k)*0.001_r_kind) ! mls pressure in Pa, coverted to log(cb) - mlsoz(k)=hdrmlsl2(2,k) ! ozone mixing ratio in ppmv - mlsozpc(k)=hdrmlsl2(3,k) ! ozone mixing ratio precision in ppmv - if (dsis /= 'mls_aura_ozpc') mlsozpc(k)=hdrmlsl2(4,k) ! use obserr if (dsis /= 'mls_aura_ozpc') - end do - - do k=1,nloz - if(hdrmlsl2(1,k)>21600._r_kind .or. hdrmlsl2(1,k)<2._r_kind) usage1(k)=1000._r_kind - if(hdrmlsl2(3,k)<=0._r_kind) usage1(k)=1000._r_kind - end do - - hdrmls13=hdrmls(13)*0.1_r_kind - if (abs(slats0)<30._r_kind) then - do k=1,nloz - if(hdrmlsl2(1,k)>10000._r_kind .and. hdrmlsl2(1,k)<21600._r_kind) then - if(hdrmls13 <= 1.2_r_kind) usage1(k)=1000._r_kind - else - if(hdrmls13 <= 0.4_r_kind) usage1(k)=1000._r_kind - endif - end do - else - if(hdrmls13 <= 0.4_r_kind) then - do k=1,nloz - usage1(k)=1000._r_kind - end do - end if - end if - - do k=1,nloz - - ndata=ndata+1 - - if(ndata <= maxobs) then - ozout( 1,ndata)=rsat - ozout( 2,ndata)=t4dv - ozout( 3,ndata)=dlon ! grid relative longitude - ozout( 4,ndata)=dlat ! grid relative latitude - ozout( 5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) - ozout( 6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) - ozout( 7,ndata)=hdrmls(10) ! solar zenith angle - - ozout( 8,ndata)=usage1(k) ! - ozout( 9,ndata)=mlspres(k) ! mls pressure in log(cb) - ozout(10,ndata)=mlsozpc(k) ! ozone mixing ratio precision in ppmv - ozout(11,ndata)=float(ipos(k)) ! pointer of obs level index in ozinfo.txt - ozout(12,ndata)=nloz ! # of mls vertical levels - ozout(13,ndata)=mlsoz(k) ! ozone mixing ratio in ppmv - endif - end do - - end do obsloop -! End of o3lev bufr loop - - nodata=min(ndata,maxobs) ! count of retained data - if(nodata < ndata) then ! warning if all are not properly retained - call warn(myname_,'insufficient buffer space, expecting =',ndata) - call warn(myname_,' actual retained =',nodata) - call warn(myname_,' size(ozout,2) =',maxobs) - endif - -! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & -! ' obstype,nmrecs,ndata,nodata,no/ndata = ',dtype,nmrecs,ndata,nodata,real(nodata)/ndata - -end subroutine ozlev_bufrread_ - -!.................................................................................. - subroutine check(status,stat) -!.................................................................................. - use netcdf, only: nf90_noerr - use netcdf, only: nf90_strerror - use mpeu_util, only: die,perr,warn - implicit none - integer(i_kind), intent (in) :: status - integer(i_kind), optional, intent(out):: stat - character(len=*),parameter:: myname_=myname//'::check' - - if(present(stat)) stat=status - - if(status /= nf90_noerr) then - if(present(stat)) then - call warn(myname_,'ignored, nf90_strerror =',trim(nf90_strerror(status))) - else - call perr(myname_,'nf90_strerror =',trim(nf90_strerror(status))) - call die(myname_) - endif - endif - end subroutine check -end module m_extOzone diff --git a/src/gsi/m_extozone.F90 b/src/gsi/m_extozone.F90 new file mode 100644 index 0000000000..69d13f9665 --- /dev/null +++ b/src/gsi/m_extozone.F90 @@ -0,0 +1,1452 @@ +module m_extozone +!#define NO_EXTRA_ +!$$$ subprogram documentation block +! . . . . +! subprogram: module m_extozone +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 900.3 +! date: 2013-09-27 +! +! abstract: a module for reading extra ozone observation data +! +! program history log: +! 2012-09-10 wargan - add omi text option for omi with efficiency factors +! 2012-12-06 jin - read mls o3lev data from bufr files. +! keep's kris's ascii reading and rename the type to 'o3levtext'. +! 2012-12-14 jin/Wargan - fixed mflg error under "o3levtext" handling block. +! 2013-01-18 Wargan - read omieffnc from a netcdf file, and removed omieff text. +! 2013-02-05 guo - stop statements were replaced with die() to signal an _abort_. +! - dfile_format() is used to handle formats of "omieff" and "o3lev". +! - code is slightly cleaned up for possible restructuring +! with alternative format handling. +! 2013-09-27 guo - added this document block. +! 2014-02-03 guo - restructured from gmao extended read_ozone() routine +! as a separated module for extended ozone obstypes, in +! particular, if not in bufr format. +! - removed "text" option of "o3lev", for it is not been +! used anymore. +! - moved history log messages here from read_ozone. +! 2015-09-17 Thomas - add l4densvar and thin4d to data selection procedure +! 2015-10-01 guo - consolidate use of ob location (in deg) +! 2016-09-19 guo - moved function dfile_format() here from obsmod.F90. +! 2018-05-21 j.jin - added time-thinning. Moved the checking of thin4d into satthin.F90. +! 2018-05-25 wargan - added omps and hooks for lims, uars mls, mipas +! +! 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 kinds, only: i_kind,r_kind,r_double + use obsmod, only: time_window_max + implicit none + private ! except + + public:: is_extozone + public:: extozone_read + +!! public:: extozone_setupoz +!! public:: extozone_setupozlev +!! public:: extozone_setupoztot + +!! public:: extozone_intoz +!! public:: extozone_intozlev +!! public:: extozone_intoztot + + interface is_extozone; module procedure is_extozone_; end interface + interface extozone_read; module procedure read_; end interface + + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + character(len=*) , parameter:: myname='m_extozone' + + integer(kind=i_kind), parameter:: maxobs_ = 1000000 + + real (kind=r_kind), parameter:: rmiss = -9999.9_r_kind + real (kind=r_kind), parameter:: badoz = 10000.0_r_kind + real (kind=r_kind), parameter:: r6 = 6.0_r_kind + real (kind=r_kind), parameter:: r360 = 360.0_r_kind + real (kind=r_kind) :: ptime,timeinflat,crit0 + integer(kind=i_kind) :: ithin_time,n_tbin,it_mesh + +contains +function is_extozone_(dfile,dtype,dplat,class) + + use mpeu_util, only: die,perr + implicit none + logical:: is_extozone_ ! this is the function return variable + character(len=*),intent(in):: dfile ! observation input filename (see gsiparm.anl:&obsinput/) + character(len=*),intent(in):: dtype ! observation type (see gsiparm.anl:&obsinput/) + character(len=*),intent(in):: dplat ! platform (see gsiparm.anl:&obsinput/) + character(len=*),optional,intent(in):: class ! specify either "level" or "layer" for sub-class + +! Use case 1: (in read_obs()) +! +! elseif(is_extozone(dfile,obstype,dplat)) then +! call read_ozone(obstype,dplat,dsis,dfile,...) +! +! +! Use case 2: (in setuprhsall()) +! +! elseif(is_extozone(dfile,obstype,dplat,class='level')) then +! call setupozlev(obstype,dplat,dsis,...) +! elseif(is_extozone(dfile,obstype,dplat,class='layer')) then +! call setupozlay(obstype,dplat,dsis,...) +! elseif(is_extozone(dfile,obstype,dplat,class='total')) then +! call setupozlay(obstype,dplat,dsis,...) +! +! Use case 3: (where intoz() or stpoz() are called) +! +! elseif(is_extozone(dfile,obstype,dplat)) then +! call intoz(obstype,dplat,dsis,...) + + character(len=*), parameter:: myname_=myname//'::is_extozone_' + + integer(kind=i_kind),parameter:: iany = 0 + integer(kind=i_kind),parameter:: iunknown=-1 + + integer(kind=i_kind),parameter:: ilevel = 1 + integer(kind=i_kind),parameter:: ilayer = 2 + integer(kind=i_kind),parameter:: itotal = 3 + + integer(kind=i_kind),parameter:: itext = 1 + integer(kind=i_kind),parameter:: ibufr = 2 + integer(kind=i_kind),parameter:: inc = 3 + + integer(kind=i_kind):: class_,ifile_ + + ifile_=iunknown + select case(dfile_format(dfile)) + case("text") + ifile_=itext + case("bufr") + ifile_=ibufr + case("nc") + ifile_=inc + end select + + class_=iany + if(present(class)) then + select case(class) + case('level') + class_=ilevel + case('layer') + class_=ilayer + case('total') + class_=itotal + case default + class_=iunknown + call perr(myname_,'unknown ozone class, class =',class) + call perr(myname_,' dfile =',dfile) + call perr(myname_,' dtype =',dtype) + call perr(myname_,' dplat =',dplat) + call die(myname_) + end select + endif + + is_extozone_= .false. +#ifndef NO_EXTRA_ + select case(class_) + case(iany) + is_extozone_= & + ifile_==ibufr .and. dtype == 'o3lev' .or. & + ifile_==inc .and. dtype == 'mls55' .or. & + ifile_==inc .and. dtype == 'ompslpvis' .or. & + ifile_==inc .and. dtype == 'ompslpuv' .or. & + ifile_==inc .and. dtype == 'ompslp' .or. & + ifile_==inc .and. dtype == 'lims' .or. & + ifile_==inc .and. dtype == 'uarsmls' .or. & + ifile_==inc .and. dtype == 'mipas' .or. & + ifile_==inc .and. dtype == 'omieff' .or. & + ifile_==inc .and. dtype == 'tomseff' + + case(ilevel) + is_extozone_= & + ifile_==ibufr .and. dtype == 'o3lev' .or. & + ifile_==inc .and. dtype == 'mls55' .or. & + ifile_==inc .and. dtype == 'ompslpvis' .or. & + ifile_==inc .and. dtype == 'ompslpuv' .or. & + ifile_==inc .and. dtype == 'ompslp' .or. & + ifile_==inc .and. dtype == 'lims' .or. & + ifile_==inc .and. dtype == 'uarsmls' .or. & + ifile_==inc .and. dtype == 'mipas' + + case(ilayer) + is_extozone_= .false. + + case(itotal) + is_extozone_= & + ifile_==inc .and. dtype == 'omieff' .or. & + ifile_==inc .and. dtype == 'tomseff' + + case default + is_extozone_= .false. + end select +#endif + +end function is_extozone_ + +! ---------------------------------------------------------------------- +function dfile_format(dfile) result(dform) +!$$$ subprogram documentation block +! . . . . +! subprogram: function dfile_format +! prgmmr: j guo +! org: NASA/GSFC, Global Modeling and Assimilation Office, 610.1 +! date: 2013-02-04 +! +! abstract: - check filename suffix to guess its format +! +! program history log: +! 2013-02-04 j guo - added this document block +! +! 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 + +! function interface: + + implicit none + + character(len=len('unknown')):: dform ! a 2-4 byte code for a format guess, + ! from a list of filename suffixes, 'bufr', 'text' (also for 'txt', + ! 'tcvitle', or 'vitl'), 'nc', or return a default value 'unknown'. One + ! can extend the list to other suffixes, such as 'hdf', 'hdf4', 'hdf5', + ! etc., if they are needed in the future. + character(len=*),intent(in):: dfile ! a given filename + +!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + character(len=*),parameter :: myname_=myname//'::dfile_format' + integer(i_kind):: i,l + + dform='unknown' + l=len_trim(dfile) + + i=max(0,l-6)+1 ! 6 byte code? + select case(dfile(i:l)) + case ('tcvitl') + dform='text' + end select + if(dform/='unknown') return + + i=max(0,l-4)+1! 4 byte code? + select case(dfile(i:l)) + case ('bufr','avhr') + dform='bufr' + case ('text','vitl') + dform='text' + end select + if(dform/='unknown') return + + i=max(0,l-3)+1 ! 3 byte code? + select case(dfile(i:l)) + case ('txt') ! a short + dform='text' + end select + if(dform/='unknown') return + + i=max(0,l-2)+1 ! 2 byte code? + select case(dfile(i:l)) + case ('nc') + dform='nc' + end select + if(dform/='unknown') return + + dform='bufr' + return +end function dfile_format + +subroutine read_(dfile,dtype,dplat,dsis, & ! intent(in), keys for type managing + nread,npuse,nouse, & ! intent(out), beside other implicit output variables + jsatid,gstime,lunout,twind,ithin,rmesh,nobs) ! intent(in), beside other implicit input variables + + use constants, only: zero + use mpeu_util, only: die,perr,tell + use mpimod, only: npe + use satthin, only: satthin_time_info => radthin_time_info +! use mpeu_util, only: mprefix,stdout +! nobs - array of observations on each subdomain for each processor + + implicit none + character(len=*), intent(in):: dfile ! obs_input filename + character(len=*), intent(in):: dtype ! observation type (see gsiparm.anl:&obsinput/) + character(len=*), intent(in):: dplat ! platform (see gsiparm.anl:&obsinput/) + character(len=*), intent(in):: dsis ! sensor/instrument/satellite tag (see gsiparm.anl:&obsinput/ and ozinfo) + +! Use case: (in read_ozone()) +! +! elseif(is_extozone(dfile,dtype,dplat)) then +! call extozone_read(dtype,dplat,dsis,dfile,...) +! + + character(len=*), parameter:: myname_=myname//'::read_' + + integer(kind=i_kind), intent(out):: nread ! number of obs record reads in this call + integer(kind=i_kind),dimension(npe), intent(inout):: nobs ! number of obs record reads in this call + integer(kind=i_kind), intent(out):: npuse ! nnmber of "preofiles" retained in this call + integer(kind=i_kind), intent(out):: nouse ! nnmber of obs retained in this call + + ! jsatid,gstime,lunout,twind,ithin,rmesh + character(len=*) , intent(in ):: jsatid ! platform id (verification) + real (kind=r_kind), intent(in ):: gstime ! analysis time (minute) from reference date + integer(kind=i_kind), intent(in ):: lunout ! logical unit to send obs output. + real (kind=r_kind), intent(in ):: twind ! input group time window (hour) + + integer(kind=i_kind), intent(in ):: ithin ! flag to thin data + real (kind=r_kind), intent(in ):: rmesh ! thining mesh size (km) + + real(kind=r_kind),pointer,dimension(:,:):: p_out + + integer(kind=i_kind):: nreal,nchan,ilat,ilon + integer(kind=i_kind):: maxobs + integer(kind=i_kind):: i,k + + nread=-1 + npuse=-1 + nouse=-1 + + if(.not.is_extozone(dfile,dtype,dplat)) then + call perr(myname_,'unexpected use, dfile =',dfile) + call perr(myname_,' dtype =',dtype) + call perr(myname_,' dplat =',dplat) + call die (myname_) + endif + call satthin_time_info(dtype, dplat, dsis, ptime, ithin_time) + if( ptime > 0.0_r_kind) then + n_tbin=nint(2.0_r_kind*time_window_max/ptime) + else + n_tbin=1 + endif + + select case(dtype) + case('omieff','tomseff') ! layer-ozone or total-ozone types + select case(dfile_format(dfile)) + case('nc') + call oztot_ncinquire_(nreal,nchan,ilat,ilon, & + ithin,rmesh,maxobs) + + allocate(p_out(nreal+nchan,maxobs)) + p_out(:,:)=rmiss + + call oztot_ncread_(dfile,dtype,dsis, p_out,nread,npuse,nouse, & + gstime, twind, ithin) + + ! Skip all "thinned" data records, and reset values of + ! nread, npuse, and nouse, as they are required by + ! upper level read_obs(). + + k=0 + do i=1,maxobs + if(p_out(1,i)>zero) then + k=k+1 + if(i>k) p_out(:,k)=p_out(:,i) + endif + enddo + nouse=k + npuse=k + end select + + case('o3lev') ! level-ozone types + select case(dfile_format(dfile)) +! case('text') +! call ozlev_textinquire(dfile,dtype,dplat, & +! nreal,nchan,ilat,ilon, maxobs) +! +! allocate(p_out(nreal+nchan,maxobs)) +! p_out(:,:)=rmiss +! +! call ozlev_textread_(dfile,dtype,dplat,dsis, p_out,nread,npuse,nouse, & +! gstime,twind, nreal,nchan,ilat,ilon) +! + case('bufr') + call ozlev_bufrinquire_(nreal,nchan,ilat,ilon,maxobs) + + allocate(p_out(nreal+nchan,maxobs)) + p_out(:,:)=rmiss + + call ozlev_bufrread_(dfile,dtype,dsis, p_out,nread,npuse,nouse, & + jsatid, gstime,twind) + end select + + case('mls55','ompslpvis','ompslpuv','ompslp','lims','uarsmls','mipas') + select case(dfile_format(dfile)) + case('nc') + call ozlev_ncinquire_( nreal,nchan,ilat,ilon,maxobs) + + allocate(p_out(nreal+nchan,maxobs)) + p_out(:,:)=rmiss + + call ozlev_ncread_(dfile,dtype, p_out,nread,npuse,nouse, gstime,twind) + + end select + + end select + + if(nouse<0 .or. .not.associated(p_out)) then + call perr(myname_,'can not process, dtype =',trim(dtype)) + call perr(myname_,' dplat =',trim(dplat)) + call perr(myname_,' dsis =',trim(dsis)) + call perr(myname_,' dfile =',trim(dfile)) + call perr(myname_,' associated(p_out) =',associated(p_out)) + if(associated(p_out)) then + call perr(myname_,' actual size(p_out,1) =',size(p_out,1)) + call perr(myname_,' actual size(p_out,2) =',size(p_out,2)) + call perr(myname_,' nread =',nread) + call perr(myname_,' npuse =',npuse) + call perr(myname_,' nouse =',nouse) + endif + call die(myname_) + endif + + if(nreal+nchan/=size(p_out,1)) then + call perr(myname_,'unexpected size(p_out,1), nreal+nchan =',nreal+nchan) + call perr(myname_,' nreal =',nreal) + call perr(myname_,' nchan =',nchan) + call perr(myname_,' actual size(p_out,1) =',size(p_out,1)) + call perr(myname_,' actual size(p_out,2) =',size(p_out,2)) + call die(myname_) + endif + +! Output candidate observations from _dfile_ to _obsfile_ (pre-opened in lunout). + +! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & +! ' obstype,nread,npuse,nouse,no/npuse = ',dtype,nread,npuse,nouse,real(nouse)/npuse + + ! While nreal+nchan defines the leading dimension of p_out(:,:), it is + ! obviously a bad idea that n_out has been missing from the header for + ! its second dimension. + + call count_obs(npuse,nreal+nchan,ilat,ilon,p_out,nobs) + write(lunout) dtype, dsis, nreal, nchan, ilat, ilon + write(lunout) p_out(:,1:npuse) + + deallocate(p_out) + +end subroutine read_ + +subroutine oztot_ncinquire_( nreal,nchan,ilat,ilon, ithin,rmesh,maxrec) + use satthin, only: satthin_itxmax => itxmax + use satthin, only: satthin_makegrids => makegrids + implicit none + + integer(kind=i_kind), intent(out):: nreal ! number of real parameters per record + integer(kind=i_kind), intent(out):: nchan ! number of channels or levels per record + integer(kind=i_kind), intent(out):: ilat ! index to latitude in nreal parameters. + integer(kind=i_kind), intent(out):: ilon ! index to longitude in nreal parameters. + + integer(kind=i_kind), intent(in ):: ithin ! flag to thin data + real (kind=r_kind), intent(in ):: rmesh ! size (km) of the thinning mesh + + integer(kind=i_kind), intent(out):: maxrec ! extimated input record count + + character(len=*), parameter:: myname_=myname//'::oztot_ncinquire_' + +! Configure the record buffer for this obs class + nreal=36 + nchan=1 + ilat=4 + ilon=3 + +! Make thinning grids, and to define the total record size. + call satthin_makegrids(rmesh,ithin,n_tbin=n_tbin) + maxrec=satthin_itxmax + +end subroutine oztot_ncInquire_ + +!.................................................................................. +subroutine oztot_ncread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, & + gstime,twind, ithin) +!.................................................................................. + + use netcdf, only: nf90_open + use netcdf, only: nf90_nowrite + use netcdf, only: nf90_noerr + use netcdf, only: nf90_inq_dimid + use netcdf, only: nf90_inquire_dimension + use netcdf, only: nf90_inq_varid + use netcdf, only: nf90_get_var + use netcdf, only: nf90_close + + use satthin, only: satthin_makegrids => makegrids + use satthin, only: satthin_tdiff2crit => tdiff2crit + use satthin, only: satthin_map2tgrid => map2tgrid + use satthin, only: satthin_finalcheck => finalcheck + use satthin, only: satthin_destroygrids => destroygrids + + use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons + use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar + use obsmod, only: nloz_omi + + use constants, only: deg2rad,zero,r60inv +! use mpeu_util, only: mprefix,stdout + + implicit none + character(len=*), intent(in):: dfile ! obs_input filename + character(len=*), intent(in):: dtype ! observation type (see gsiparm.anl:&obsinput/) + character(len=*), intent(in):: dsis ! sensor/instrument/satellite tag (see gsiparm.anl:&obsinput/ and ozinfo) + + real (kind=r_kind), dimension(:,:), intent(out):: ozout + integer(kind=i_kind), intent(out):: nmrecs ! read count in ozout + integer(kind=i_kind), intent(out):: ndata ! "good" profile count in ozout + integer(kind=i_kind), intent(out):: nodata ! "good" data count in ozout + + real (kind=r_kind), intent(in):: gstime ! analysis time (minute) from reference date + real (kind=r_kind), intent(in):: twind ! input group time window (hour) + + + integer(kind=i_kind), intent(in ):: ithin ! flag to thin data + + character(len=*), parameter:: myname_=myname//'::oztot_ncread_' + +! parameters for output bookkeeping + integer(kind=i_kind):: i, irec + integer(kind=i_kind):: itx,itt + +! parameters for NetCDF arrays + integer(kind=i_kind), allocatable :: iya(:),ima(:),idda(:),ihha(:),imina(:),fovna(:) + integer(kind=i_kind), allocatable :: toqfa(:),alqfa(:) + real (kind=r_kind), allocatable :: rseca(:),slatsa(:),slonsa(:),totoza(:),szaa(:) + real (kind=r_kind), allocatable :: aprioria(:,:),efficiencya(:,:) + integer(kind=i_kind) nrecdimid,nomilevsdimid,lonvarid,latvarid,yyvarid,mmvarid + integer(kind=i_kind) ddvarid,hhvarid,minvarid,ssvarid,fovnvarid,toqfvarid,alqfvarid + integer(kind=i_kind) szavarid,totozvarid,apriorivarid,efficiencyvarid + integer(kind=i_kind) ier, ncid, nomilevs + integer(kind=i_kind):: iy,im,idd,ihh,imin,nmind + integer(kind=i_kind):: toqf,alqf,fovn + real (kind=r_kind):: rsec,slats,slons + real (kind=r_kind),dimension(nloz_omi):: apriori, efficiency + integer(kind=i_kind):: binary(17) + + real (kind=r_kind):: dlon,dlon_earth,dlon_earth_deg + real (kind=r_kind):: dlat,dlat_earth,dlat_earth_deg + real (kind=r_kind):: tdiff,sstime,t4dv,crit1,dist1,rsat + integer(kind=i_kind):: idate5(5) + + real (kind=r_double):: totoz, sza + + logical:: outside, removescans, iuse + integer(kind=i_kind):: maxobs + + + removescans = .true. ! remove omi scan numbers >= 25? + maxobs=size(ozout,2) + rsat=999._r_kind + +! Using omi/toms with efficience factors in netcdf format + nodata = 0 + ndata = 0 + + call check(nf90_open(trim(dfile),nf90_nowrite,ncid),stat=ier) + + ! ignore if the file is not actually present. + if(ier/=nf90_noerr) then + call satthin_destroygrids() + return + end if + + ! Get dimensions from omi input file + call check(nf90_inq_dimid(ncid, "nrec", nrecdimid),stat=ier) + + ! ignore if the file header is empty + if(ier/=nf90_noerr) then + call check(nf90_close(ncid),stat=ier) + call satthin_destroygrids() + return + endif + + ! Get dimensions from omi/toms input file +! nmrecs=0 + call check(nf90_inquire_dimension(ncid, nrecdimid, len = nmrecs),stat=ier) + + ! ignore if the file header is empty + if(ier/=nf90_noerr .or. nmrecs==0) then + call check(nf90_close(ncid),stat=ier) + call satthin_destroygrids() + return + endif + + ! Continue the input + call check(nf90_inq_dimid(ncid, "nlevs", nomilevsdimid)) + call check(nf90_inquire_dimension(ncid, nomilevsdimid, len = nomilevs)) + + ! We have dimensions so we can allocate arrays + allocate(iya(nmrecs),ima(nmrecs),idda(nmrecs),ihha(nmrecs),imina(nmrecs), & + rseca(nmrecs),fovna(nmrecs),slatsa(nmrecs),slonsa(nmrecs),totoza(nmrecs), & + toqfa(nmrecs),alqfa(nmrecs),szaa(nmrecs)) + allocate(aprioria(nomilevs,nmrecs),efficiencya(nomilevs,nmrecs)) + + ! Read variables and store them in these arrays + call check(nf90_inq_varid(ncid, "lon", lonvarid)) + call check(nf90_get_var(ncid, lonvarid, slonsa)) + + call check(nf90_inq_varid(ncid, "lat", latvarid)) + call check(nf90_get_var(ncid, latvarid, slatsa)) + + call check(nf90_inq_varid(ncid, "yy", yyvarid)) + call check(nf90_get_var(ncid, yyvarid, iya)) + + call check(nf90_inq_varid(ncid, "mm", mmvarid)) + call check(nf90_get_var(ncid, mmvarid, ima)) + + call check(nf90_inq_varid(ncid, "dd", ddvarid)) + call check(nf90_get_var(ncid, ddvarid, idda)) + + call check(nf90_inq_varid(ncid, "hh", hhvarid)) + call check(nf90_get_var(ncid, hhvarid, ihha)) + + call check(nf90_inq_varid(ncid, "min", minvarid)) + call check(nf90_get_var(ncid, minvarid, imina)) + + call check(nf90_inq_varid(ncid, "ss", ssvarid)) + call check(nf90_get_var(ncid, ssvarid, rseca)) + + call check(nf90_inq_varid(ncid, "fov", fovnvarid)) + call check(nf90_get_var(ncid, fovnvarid, fovna)) + + call check(nf90_inq_varid(ncid, "sza", szavarid)) + call check(nf90_get_var(ncid, szavarid, szaa)) + + call check(nf90_inq_varid(ncid, "toqf", toqfvarid)) + call check(nf90_get_var(ncid, toqfvarid, toqfa)) + + call check(nf90_inq_varid(ncid, "alqf", alqfvarid)) + call check(nf90_get_var(ncid, alqfvarid, alqfa)) + + call check(nf90_inq_varid(ncid, "toz", totozvarid)) + call check(nf90_get_var(ncid, totozvarid, totoza)) + + call check(nf90_inq_varid(ncid, "apriori", apriorivarid)) + call check(nf90_get_var(ncid, apriorivarid, aprioria)) + + call check(nf90_inq_varid(ncid, "efficiency", efficiencyvarid)) + call check(nf90_get_var(ncid, efficiencyvarid, efficiencya)) + + ! close the data file + call check(nf90_close(ncid)) + + ! now screen the data and put them into the right places + recloop: do irec = 1, nmrecs + iy = iya(irec) + im = ima(irec) + idd = idda(irec) + ihh = ihha(irec) + imin = imina(irec) + rsec = rseca(irec) + fovn = fovna(irec) + slats = slatsa(irec) + slons = slonsa(irec) + totoz = totoza(irec) + toqf = toqfa(irec) + alqf = alqfa(irec) + sza = szaa(irec) + do i = 1, nomilevs + apriori(i) = aprioria(i, irec) + ! Reported efficiencies from layer 4 up (counting from the surface + ! so from 1 - 8 here) are incorrect (too large) for high sza + ! Setting them all to 1.0 per pk bhartia''s recommendation + if (i <= 8) then + efficiency(i) = 1._r_kind + else + efficiency(i) = efficiencya(i, irec) + endif + end do + +! nmrecs=nmrecs+1 + ! Convert observation location to radians + if(abs(slats)>90._r_kind .or. abs(slons)>r360) cycle recloop + if(slons< zero) slons=slons+r360 + if(slons==r360) slons=zero + dlat_earth_deg = slats + dlon_earth_deg = slons + dlat_earth = slats * deg2rad + dlon_earth = slons * deg2rad + + if(regional)then + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) + if(outside) cycle recloop + else + dlat = dlat_earth + dlon = dlon_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + +! convert observation time to relative time + idate5(1) = iy !year + idate5(2) = im !month + idate5(3) = idd !day + idate5(4) = ihh !hour + idate5(5) = imin !minute + call w3fs21(idate5,nmind) + + t4dv=real((nmind-iwinbgn),r_kind)*r60inv + sstime=real(nmind,r_kind) + tdiff=(sstime-gstime)*r60inv + if (l4dvar.or.l4densvar) then + if (t4dvwinlen) cycle recloop + else + if(abs(tdiff) > twind) cycle recloop + end if + + if (totoz > badoz ) cycle recloop + +! Apply data screening based on quality flags +! Bit 10 (from the left) in TOQF represents row anomaly. All 17 bits in toqf is +! supposed to converted into array elements of binary(:), either for "tomseff" or +! "omieff". + binary(:)=0 + select case(dtype) + case('omieff') + + if (toqf /= 0 .and. toqf /= 1) cycle recloop + +! Remove obs at high solar zenith angles + if (sza > 84.0_r_kind) cycle recloop + +! remove the bad scan position data: fovn beyond 25 + if (removescans) then + if (fovn >=25) cycle recloop + endif + if (fovn <=2 .or. fovn >=58) cycle recloop + +! remove the data in which the c-pair algorithm ((331 and 360 nm) is used. + if (alqf == 3 .or. alqf == 13) cycle recloop + + case('tomseff') + ! The meaning of quality flags for toms version 8 is similar to that + ! for sbuv: + ! 0 - good data, 1 - good data with sza > 84 deg + if (toqf /= 0) cycle recloop + + case default + end select + +! thin omi/toms data + + crit0 = 0.01_r_kind + timeinflat=r6 + call satthin_tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh) + if (ithin /= 0) then + call satthin_map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,dsis,it_mesh=it_mesh) + if(.not. iuse)cycle recloop + call satthin_finalcheck(dist1,crit1,itx,iuse) + if(.not. iuse)cycle recloop + ndata=ndata+1 + nodata=ndata + else + ndata=ndata+1 + nodata=ndata + itx= ndata + end if + +! ASSERT_(size(ozout,2)>=itx) + + if(itx <= maxobs) then + ozout(1,itx)=rsat + ozout(2,itx)=t4dv + ozout(3,itx)=dlon ! grid relative longitude + ozout(4,itx)=dlat ! grid relative latitude + ozout(5,itx)=dlon_earth_deg ! earth relative longitude (degrees) + ozout(6,itx)=dlat_earth_deg ! earth relative latitude (degrees) + ozout(7,itx)=real(toqf) ! - total ozone quality code (not used) + ozout(8,itx)=real(sza) ! solar zenith angle + ozout(9,itx)=binary(10) ! row anomaly flag, is actually fixed to 0 + ozout(10,itx)=0._r_kind ! - cloud amount (not used) + ozout(11,itx)=0._r_kind ! - vzan (not used) + ozout(12,itx)=0._r_kind ! - aerosol index (not used) + ozout(13,itx)=0._r_kind ! - ascending/descending (not used) + ozout(14,itx)=real(fovn) ! scan position + ! "(not used)" flags above imply that they + ! are not used in setupozlay(). + +! Added apriori and efficiency profiles + ozout(15:25,itx)=apriori + ozout(26:36,itx)=efficiency + ozout(37,itx)=totoz + endif + +! ASSERT_(size(ozout,1)==37) + + end do recloop + + deallocate(iya,ima,idda,ihha,imina, & + rseca,fovna,slatsa,slonsa,totoza, & + toqfa,alqfa,szaa,aprioria,efficiencya) +! end +! End of loop over observations +! End of omi block with efficiency factors in netcdf format + + call satthin_destroygrids() + + nodata=min(ndata,maxobs) +! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & +! ' obstype,nmrecs,ndata,nodata,no/ndata = ',dtype,nmrecs,ndata,nodata,real(nodata)/ndata + return + +end subroutine oztot_ncread_ +!.................................................................................. + + +subroutine ozlev_ncinquire_( nreal,nchan,ilat,ilon, maxrec) + implicit none + + integer(kind=i_kind), intent(out):: nreal ! number of real parameters per record + integer(kind=i_kind), intent(out):: nchan ! number of channels or levels per record + integer(kind=i_kind), intent(out):: ilat ! index to latitude in nreal parameters. + integer(kind=i_kind), intent(out):: ilon ! index to longitude in nreal parameters. + + integer(kind=i_kind), intent(out):: maxrec ! extimated input record count + + character(len=*), parameter:: myname_=myname//'::ozlev_ncinquire_' + +! Configure the record, they are not (dfile,dtype,dplat) dependent in this case. + nreal = 12 + nchan = 1 ! There are 'levs' levels but each is treated + ! as a separate observation so that nchanl = 1 + ilat=4 + ilon=3 + + maxrec = maxobs_ +end subroutine ozlev_ncinquire_ + +!.................................................................................. +subroutine ozlev_ncread_(dfile,dtype,ozout,nmrecs,ndata,nodata, gstime,twind) +!.................................................................................. + use netcdf, only: nf90_open + use netcdf, only: nf90_nowrite + use netcdf, only: nf90_noerr + use netcdf, only: nf90_inq_dimid + use netcdf, only: nf90_inquire_dimension + use netcdf, only: nf90_inq_varid + use netcdf, only: nf90_get_var + use netcdf, only: nf90_close + + use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons + use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar + + use constants, only: deg2rad,zero,one_tenth,r60inv + use ozinfo, only: jpch_oz,nusis_oz,iuse_oz + use mpeu_util, only: perr,die +! use mpeu_util, only: mprefix,stdout + + implicit none + character(len=*), intent(in):: dfile ! obs_input filename + character(len=*), intent(in):: dtype ! obs_input dtype + + real (kind=r_kind), dimension(:,:), intent(out):: ozout + integer(kind=i_kind), intent(out):: nmrecs ! count of records read + integer(kind=i_kind), intent(out):: ndata ! count of processed + integer(kind=i_kind), intent(out):: nodata ! count of retained + + real (kind=r_kind), intent(in):: gstime ! analysis time (minute) from reference date + real (kind=r_kind), intent(in):: twind ! input group time window (hour) + + + character(len=*), parameter:: myname_=myname//'::ozlev_ncread_' + + integer(kind=i_kind):: ier, iprof, nprofs, maxobs + integer(kind=i_kind):: i, ilev, ikx, ncid, k0 + integer(kind=i_kind),allocatable,dimension(:):: ipos + + integer(kind=i_kind):: nrecdimid,lonvarid,latvarid,yyvarid,mmvarid + integer(kind=i_kind):: ddvarid,hhvarid,minvarid,ssvarid + integer(kind=i_kind):: pressvarid + integer(kind=i_kind):: convvarid, errvarid, ozonevarid + integer(kind=i_kind):: levsdimid,levs + + integer(kind=i_kind), allocatable :: iya(:),ima(:),idda(:),ihha(:),imina(:),iseca(:) + real (kind=r_kind), allocatable :: slatsa(:),slonsa(:) + real (kind=r_kind), allocatable :: press(:), ozone(:,:), qual(:), press2d(:,:) + real (kind=r_kind), allocatable :: conv(:), err(:,:) + + integer(kind=i_kind):: nmind + + real (kind=r_kind):: slons0,slats0 + real (kind=r_kind):: ppmv, pres, pob, obserr, usage + + real (kind=r_kind):: dlon,dlon_earth,dlon_earth_deg + real (kind=r_kind):: dlat,dlat_earth,dlat_earth_deg + real (kind=r_kind):: tdiff,sstime,t4dv,rsat + integer(kind=i_kind):: idate5(5) + + logical:: outside + logical:: first + + maxobs=size(ozout,2) + rsat=999._r_kind + ndata = 0 + nmrecs=0 + nodata=-1 +!.................................................................................. + ! ---------------mls and other limb ozone data in netcdf format ---------- + ! Open file and read dimensions + call check(nf90_open(trim(dfile),nf90_nowrite,ncid),stat=ier) + + ! ignore if the file is not actually present. + if(ier/=nf90_noerr) return + + ! Get dimensions from omi input file + call check(nf90_inq_dimid(ncid, "nprofiles", nrecdimid),stat=ier) + + ! ignore if the file header is empty + if(ier/=nf90_noerr) then + call check(nf90_close(ncid),stat=ier) + return + endif + + ! Get dimensions from the input file: # of profiles and # of levels + nprofs=0 + call check(nf90_inquire_dimension(ncid, nrecdimid, len = nprofs),stat=ier) + + ! ignore if the file header is empty + if(ier/=nf90_noerr) then + call check(nf90_close(ncid),stat=ier) + return + endif + + if(nprofs==0) then + nodata=0 + call check(nf90_close(ncid),stat=ier) + return + endif + + ! Continue the input + call check(nf90_inq_dimid(ncid, "nlevs", levsdimid)) + call check(nf90_inquire_dimension(ncid, levsdimid, len = levs)) + + ! Note: Make sure that 'ozinfo' has the same number of levels + ! for nrt it is 55 + allocate(ipos(levs)) + ipos=999 + + ! Process limb data in netcdf format + ikx = 0 + first=.false. + do i=1,jpch_oz + if( (.not. first) .and. index(nusis_oz(i), trim(dtype))/=0) then + k0=i + first=.true. + end if + if(first.and.index(nusis_oz(i),trim(dtype))/=0) then + ikx=ikx+1 + ipos(ikx)=k0+ikx-1 + end if + end do + + if (ikx/=levs) call die(myname_//': inconsistent levs for '//dtype) + + nmrecs=0 + ! Allocate space and read data + + allocate(iya(nprofs),ima(nprofs),idda(nprofs),ihha(nprofs),imina(nprofs), & + iseca(nprofs),slatsa(nprofs),slonsa(nprofs), ozone(levs,nprofs), & + err(levs,nprofs),qual(nprofs), conv(nprofs)) + if (index(dtype, 'ompslp') == 0) then + allocate(press(levs)) + else + allocate(press2d(levs,nprofs)) + endif + + ! Read variables and store them in these arrays + call check(nf90_inq_varid(ncid, "lon", lonvarid)) + call check(nf90_get_var(ncid, lonvarid, slonsa)) + + call check(nf90_inq_varid(ncid, "lat", latvarid)) + call check(nf90_get_var(ncid, latvarid, slatsa)) + + call check(nf90_inq_varid(ncid, "yy", yyvarid)) + call check(nf90_get_var(ncid, yyvarid, iya)) + + call check(nf90_inq_varid(ncid, "mm", mmvarid)) + call check(nf90_get_var(ncid, mmvarid, ima)) + + call check(nf90_inq_varid(ncid, "dd", ddvarid)) + call check(nf90_get_var(ncid, ddvarid, idda)) + + call check(nf90_inq_varid(ncid, "hh", hhvarid)) + call check(nf90_get_var(ncid, hhvarid, ihha)) + + call check(nf90_inq_varid(ncid, "min", minvarid)) + call check(nf90_get_var(ncid, minvarid, imina)) + + call check(nf90_inq_varid(ncid, "ss", ssvarid)) + call check(nf90_get_var(ncid, ssvarid, iseca)) + + if (index(dtype, 'ompslp') == 0) then + call check(nf90_inq_varid(ncid, "press", pressvarid)) + call check(nf90_get_var(ncid, pressvarid, press)) + else + call check(nf90_inq_varid(ncid, "press", pressvarid)) + call check(nf90_get_var(ncid, pressvarid, press2d)) + endif + + if (index(dtype, 'ompslp') == 0) then + call check(nf90_inq_varid(ncid, "conv", convvarid)) + call check(nf90_get_var(ncid, convvarid, conv)) + end if + + !call check(nf90_inq_varid(ncid, "qual", qualvarid)) + !call check(nf90_get_var(ncid, qualvarid, qual)) + + call check(nf90_inq_varid(ncid, "oberr", errvarid)) + call check(nf90_get_var(ncid, errvarid, err)) + + call check(nf90_inq_varid(ncid, "ozone", ozonevarid)) + call check(nf90_get_var(ncid, ozonevarid, ozone)) + + ! close the data file + call check(nf90_close(ncid)) + + ! 'Unpack' the data + nmrecs = 0 + nodata = 0 + do iprof = 1, nprofs + do ilev = 1, levs + ! note that most of the quality control is done at the + ! pre-processing stage +! if (press(ilev) > 262.0_r_kind .or. press(ilev) < 0.1_r_kind ) cycle ! undefined + if (ozone(ilev, iprof) < -900.0_r_kind) cycle ! undefined + if (err(ilev, iprof) < -900.0_r_kind) cycle ! undefined + if (iuse_oz(ipos(ilev)) < 0) then + usage = 100._r_kind + else + usage = zero + endif + nmrecs=nmrecs+1 + + ! convert observation location to radians + slons0=slonsa(iprof) + slats0=slatsa(iprof) + if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle + if(slons0< zero) slons0=slons0+r360 + if(slons0==r360) slons0=zero + dlat_earth_deg = slats0 + dlon_earth_deg = slons0 + dlat_earth = slats0 * deg2rad + dlon_earth = slons0 * deg2rad + + if(regional)then + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) + if(outside) cycle + else + dlat = dlat_earth + dlon = dlon_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + + idate5(1) = iya(iprof) !year + idate5(2) = ima(iprof) !month + idate5(3) = idda(iprof) !day + idate5(4) = ihha(iprof) !hour + idate5(5) = imina(iprof) !minute + call w3fs21(idate5,nmind) + t4dv=real((nmind-iwinbgn),r_kind)*r60inv + if (l4dvar.or.l4densvar) then + if (t4dvwinlen) then + write(6,*)'read_ozone: ', dtype,' obs time idate5=',idate5,', t4dv=',& + t4dv,' is outside time window, sstime=',sstime*r60inv + cycle + end if + else + sstime=real(nmind,r_kind) + tdiff=(sstime-gstime)*r60inv + if(abs(tdiff) > twind)then + write(6,*)'read_ozone: ',dtype,' obs time idate5=',idate5,', tdiff=',& + tdiff,' is outside time window=',twind + cycle + end if + end if + + obserr = err(ilev, iprof) + ppmv = ozone(ilev, iprof) + if (index(dtype, 'ompslp') == 0) then + pres = press(ilev) + else + pres = press2d(ilev, iprof) + end if + pob = log(pres * one_tenth) + ndata = ndata+1 + if(ndata<=maxobs) then + nodata = nodata + 1 + ozout(1,ndata)=rsat + ozout(2,ndata)=t4dv + ozout(3,ndata)=dlon ! grid relative longitude + ozout(4,ndata)=dlat ! grid relative latitude + ozout(5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) + ozout(6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) + + ozout(7,ndata)=rmiss ! used to be solar zenith angle + ozout(8,ndata)=usage + ozout(9,ndata)=pob ! pressure + ozout(10,ndata)=obserr ! ozone mixing ratio precision in ppmv + ozout(11,ndata)=real(ipos(ilev)) ! pointer of obs level index in ozinfo.txt + ozout(12,ndata)=levs ! # of vertical levels + ozout(13,ndata)=ppmv ! ozone mixing ratio in ppmv + endif + + end do + end do + + + deallocate(iya,ima,idda,ihha,imina,iseca,slatsa,slonsa, ozone, & + err,qual, conv) + if (index(dtype, 'ompslp') == 0) then + deallocate(press) + else + deallocate(press2d) + end if + deallocate(ipos) + + ! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & + ! ' obstype,nmrecs,ndata,nodata,no/ndata = ',dtype,nmrecs,ndata,nodata,real(nodata)/ndata + + if (ndata > maxobs) then + call perr(myname_,'Number of limb obs reached maxobs = ', maxobs) + call perr(myname_,' ndata = ', ndata) + call perr(myname_,' nodata = ', nodata) + call die(myname_) + endif + + !---------------End mls nrt netcdf--------------------------- +end subroutine ozlev_ncread_ + +subroutine ozlev_bufrinquire_(nreal,nchan,ilat,ilon,maxrec) + implicit none + + integer(kind=i_kind), intent(out):: nreal ! number of real parameters per record + integer(kind=i_kind), intent(out):: nchan ! number of channels or levels per record + integer(kind=i_kind), intent(out):: ilat ! index to latitude in nreal parameters. + integer(kind=i_kind), intent(out):: ilon ! index to longitude in nreal parameters. + + integer(kind=i_kind), intent(out):: maxrec ! extimated input record count + + character(len=*), parameter:: myname_=myname//'::ozlev_bufrinquire_' + +! Configure the record, they are not (dfile,dtype,dplat) dependent in this case. + nreal = 12 + nchan = 1 + ilat=4 + ilon=3 + + maxrec = maxobs_ +end subroutine ozlev_bufrinquire_ + +subroutine ozlev_bufrread_(dfile,dtype,dsis, ozout,nmrecs,ndata,nodata, & + jsatid, gstime,twind) + + use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons + use gsi_4dvar, only: l4dvar,iwinbgn,winlen,l4densvar + + use constants, only: deg2rad,zero,r60inv + use ozinfo, only: jpch_oz,nusis_oz,iuse_oz + use radinfo, only: dec2bin + + use mpeu_util, only: warn,tell +! use mpeu_util, only: mprefix,stdout + + implicit none + character(len=*), intent(in):: dfile ! obs_input filename + character(len=*), intent(in):: dtype ! observation type (see gsiparm.anl:&obsinput/) + character(len=*), intent(in):: dsis ! sensor/instrument/satellite tag (see gsiparm.anl:&obsinput/ and ozinfo) + + real (kind=r_kind), dimension(:,:), intent(out):: ozout + integer(kind=i_kind), intent(out):: nmrecs ! count of actual records read + integer(kind=i_kind), intent(out):: ndata ! count of processed data + integer(kind=i_kind), intent(out):: nodata ! count of retained data + + character(len=*) , intent(in):: jsatid ! platform id (verification) + real (kind=r_kind), intent(in):: gstime ! analysis time (minute) from reference date + real (kind=r_kind), intent(in):: twind ! input group time window (hour) + + + character(len=*),parameter:: myname_=myname//'::ozlev_bufrread_' + + integer(kind=i_kind),parameter:: nloz =37 + integer(kind=i_kind),parameter:: lunin=10 + + character(len=*), parameter:: mlstr ='SAID CLAT CLON YEAR MNTH DAYS HOUR MINU SECO SOZA CONV MLST PCCF' + character(len=*), parameter:: mlstrl2='PRLC OZMX OZMP OZME' + + integer(kind=i_kind):: idate,jdate,ksatid,iy,iret,im,ihh,idd + integer(kind=i_kind):: mpos + character(len= 8):: subset + + integer(kind=i_kind):: maxobs + integer(kind=i_kind):: nmind + integer(kind=i_kind):: k + integer(kind=i_kind):: kidsat + integer(kind=i_kind):: idate5(5) + integer(kind=i_kind):: ikx + integer(kind=i_kind):: decimal,binary_mls(18) + + real(kind=r_kind):: tdiff,sstime,dlon,dlat,t4dv + real(kind=r_kind):: slons0,slats0,rsat,dlat_earth,dlon_earth + real(kind=r_kind):: dlat_earth_deg,dlon_earth_deg + real(kind=r_double),dimension(13):: hdrmls + real(kind=r_double),dimension(4,37):: hdrmlsl2 + real(kind=r_double):: hdrmls13 + real(kind=r_kind),allocatable,dimension(:):: mlspres,mlsoz,mlsozpc,usage1 + integer(kind=i_kind),allocatable,dimension(:):: ipos + integer(kind=i_kind):: iprofs,nprofs + logical:: outside + + maxobs=size(ozout,2) + rsat=999._r_kind + + open(lunin,file=dfile,form='unformatted') + + call openbf(lunin,'IN',lunin) + call datelen(10) + call readmg(lunin,subset,idate,iret) + + ! If it has failed at the first read(), ... + if (iret/=0 .or. subset /= 'GM008015') then + call closbf(lunin) + close(lunin) + + call warn(myname_,'Failed at reading BUFR file, dfile =',trim(dfile)) + call warn(myname_,' dtype =',trim(dtype)) + call warn(myname_,' jsatid =',trim(jsatid)) + call warn(myname_,' lunin =',lunin) + call warn(myname_,' iret =',iret) + call warn(myname_,' subset =',trim(subset)) + + nprofs = 0 + nodata = 0 + + return + endif + + call tell(myname_,'MLS o3lev data type, subset =',trim(subset)) + write(6,*)'read_ozone: GMAO MLS o3lev data type, subset=',subset + + ! Q: o3lev data has 37 levels. Then, why is size(ipos) of 44? + ! A: Because there are 44 entries in ozinfo.txt at the time + + mpos=max(nloz,jpch_oz) + allocate (ipos(mpos)) ! 44? 37? + allocate (usage1(nloz)) + + ipos=999 + nmrecs=0 + nprofs=0 + nodata=0 + +! Set dependent variables and allocate arrays + + allocate (mlspres(nloz)) + allocate (mlsoz(nloz)) + allocate (mlsozpc(nloz)) + + ikx=0 + do k=1,jpch_oz + if(index(nusis_oz(k),'o3lev')/=0) then ! all "o3lev" in ozinfo.txt + ikx=ikx+1 + ipos(ikx)=k + end if + end do + + iy=0 + im=0 + idd=0 + ihh=0 + +! This is the top of the profile loop + ndata=0 + iprofs=0 + obsloop: do + call readsb(lunin,iret) + if (iret/=0) then !JJJ, end of the subset + call readmg(lunin,subset,jdate,iret) !JJJ open a new mg + if (iret/=0) exit obsloop !JJJ, no more mg, EOF + cycle obsloop + endif + + do k=1,nloz + if (iuse_oz(ipos(k)) < 0) then + usage1(k) = 100._r_kind + else + usage1(k) = zero + endif + end do + +! extract header information + call ufbint(lunin,hdrmls,13,1,iret,mlstr) + rsat = hdrmls(1); ksatid=rsat + + if(jsatid == 'aura')kidsat = 785 + if (ksatid /= kidsat) cycle obsloop + + nmrecs=nmrecs+nloz + +! Convert observation location to radians + slats0= hdrmls(2) + slons0= hdrmls(3) + if(abs(slats0)>90._r_kind .or. abs(slons0)>r360) cycle obsloop + if(slons0< zero) slons0=slons0+r360 + if(slons0==r360) slons0=zero + dlat_earth_deg = slats0 + dlon_earth_deg = slons0 + dlat_earth = slats0 * deg2rad + dlon_earth = slons0 * deg2rad + + if(regional)then + call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside) + if(outside) cycle obsloop + else + dlat = dlat_earth + dlon = dlon_earth + call grdcrd1(dlat,rlats,nlat,1) + call grdcrd1(dlon,rlons,nlon,1) + endif + +! convert observation time to relative time + idate5(1) = hdrmls(4) !year + idate5(2) = hdrmls(5) !month + idate5(3) = hdrmls(6) !day + idate5(4) = hdrmls(7) !hour + idate5(5) = hdrmls(8) !minute + call w3fs21(idate5,nmind) + + t4dv=real((nmind-iwinbgn),r_kind)*r60inv + if (l4dvar.or.l4densvar) then + if (t4dvwinlen) then + write(6,*)'read_ozone: mls obs time idate5=',idate5,', t4dv=',& + t4dv,' is outside time window, sstime=',sstime*r60inv + cycle obsloop + endif + else + sstime=real(nmind,r_kind) + tdiff=(sstime-gstime)*r60inv + if (abs(tdiff) > twind) then + write(6,*)'read_ozone: mls obs time idate5=',idate5,', tdiff=',& + tdiff,' is outside time window=',twind + cycle obsloop + endif + end if + +! v2.2 data screening, only accept: +! Pressure range: 215-0.02mb +! Precision: positive ozmp; +! Status flag: only use even number +! Quality(pccf): use >1.2 for data at 215-100mb & low latitude, +! use >0.4 for data elsewhere +! Convergence: use <1.8 + +! Bit 1 in mlst represents data should not be used +! Note: in bufr bits are defined from left to right as: 123456789... +! whereas in hdf5 (and the nasa document) bits are defined from right to left as: ...876543210 + + decimal=int(hdrmls(12)) + call dec2bin(decimal,binary_mls,18) + if (binary_mls(1) == 1 ) cycle obsloop + + if(hdrmls(11) >= 1.8_r_kind) cycle obsloop + +! extract pressure, ozone mixing ratio and precision + call ufbrep(lunin,hdrmlsl2,4,nloz,iret,mlstrl2) + + iprofs=iprofs+1 ! counting the profiles + do k=1,nloz + mlspres(k)=log(hdrmlsl2(1,k)*0.001_r_kind) ! mls pressure in Pa, coverted to log(cb) + mlsoz(k)=hdrmlsl2(2,k) ! ozone mixing ratio in ppmv + mlsozpc(k)=hdrmlsl2(3,k) ! ozone mixing ratio precision in ppmv + if (dsis /= 'mls_aura_ozpc') mlsozpc(k)=hdrmlsl2(4,k) ! use obserr if (dsis /= 'mls_aura_ozpc') + end do + + do k=1,nloz + if(hdrmlsl2(1,k)>21600._r_kind .or. hdrmlsl2(1,k)<2._r_kind) usage1(k)=1000._r_kind + if(hdrmlsl2(3,k)<=0._r_kind) usage1(k)=1000._r_kind + end do + + hdrmls13=hdrmls(13)*0.1_r_kind + if (abs(slats0)<30._r_kind) then + do k=1,nloz + if(hdrmlsl2(1,k)>10000._r_kind .and. hdrmlsl2(1,k)<21600._r_kind) then + if(hdrmls13 <= 1.2_r_kind) usage1(k)=1000._r_kind + else + if(hdrmls13 <= 0.4_r_kind) usage1(k)=1000._r_kind + endif + end do + else + if(hdrmls13 <= 0.4_r_kind) then + do k=1,nloz + usage1(k)=1000._r_kind + end do + end if + end if + + do k=1,nloz + + ndata=ndata+1 + + if(ndata <= maxobs) then + ozout( 1,ndata)=rsat + ozout( 2,ndata)=t4dv + ozout( 3,ndata)=dlon ! grid relative longitude + ozout( 4,ndata)=dlat ! grid relative latitude + ozout( 5,ndata)=dlon_earth_deg ! earth relative longitude (degrees) + ozout( 6,ndata)=dlat_earth_deg ! earth relative latitude (degrees) + ozout( 7,ndata)=hdrmls(10) ! solar zenith angle + + ozout( 8,ndata)=usage1(k) ! + ozout( 9,ndata)=mlspres(k) ! mls pressure in log(cb) + ozout(10,ndata)=mlsozpc(k) ! ozone mixing ratio precision in ppmv + ozout(11,ndata)=real(ipos(k)) ! pointer of obs level index in ozinfo.txt + ozout(12,ndata)=nloz ! # of mls vertical levels + ozout(13,ndata)=mlsoz(k) ! ozone mixing ratio in ppmv + endif + end do + + end do obsloop +! End of o3lev bufr loop + + nodata=min(ndata,maxobs) ! count of retained data + if(nodata < ndata) then ! warning if all are not properly retained + call warn(myname_,'insufficient buffer space, expecting =',ndata) + call warn(myname_,' actual retained =',nodata) + call warn(myname_,' size(ozout,2) =',maxobs) + endif + +! write(stdout,'(3a,3i8,f8.2)') mprefix('read_ozone'), & +! ' obstype,nmrecs,ndata,nodata,no/ndata = ',dtype,nmrecs,ndata,nodata,real(nodata)/ndata + +end subroutine ozlev_bufrread_ + +!.................................................................................. +subroutine check(status,stat) +!.................................................................................. + use netcdf, only: nf90_noerr + use netcdf, only: nf90_strerror + use mpeu_util, only: die,perr,warn + implicit none + integer(i_kind), intent (in) :: status + integer(i_kind), optional, intent(out):: stat + character(len=*),parameter:: myname_=myname//'::check' + + if(present(stat)) stat=status + + if(status /= nf90_noerr) then + if(present(stat)) then + call warn(myname_,'ignored, nf90_strerror =',trim(nf90_strerror(status))) + else + call perr(myname_,'nf90_strerror =',trim(nf90_strerror(status))) + call die(myname_) + endif + endif +end subroutine check +end module m_extozone diff --git a/src/gsi/m_find.f90 b/src/gsi/m_find.f90 index 3a667bf31c..657e5302ff 100644 --- a/src/gsi/m_find.f90 +++ b/src/gsi/m_find.f90 @@ -2,15 +2,15 @@ module m_find !$$$ module documentation block ! . . . . -! module: m_find similar to IDL's "where" function, but needs to know the +! module: m_find similar to idl's "where" function, but needs to know the ! dimension of the outupt e.g., index=find(mask,count) ! prgmmr: eliu ! -! abstract: module similar to IDL's "where" function +! abstract: module similar to idl's "where" function ! ! program history log: -! 1997-08-13 Joiner - initial coding from NASA/GMAO -! 2012-02-15 eliu - reformat to use in GSI +! 1997-08-13 Joiner - initial coding from nasa/gmao +! 2012-02-15 eliu - reformat to use in gsi ! ! subroutines included: ! @@ -25,8 +25,8 @@ module m_find use kinds,only : i_kind, r_kind interface find - module procedure find2 - module procedure find1 + module procedure find2 + module procedure find1 end interface contains @@ -107,10 +107,10 @@ function find1(mask ) result (index) index=temp(1) else !if (present(iprt)) then - ! if (iprt >= 3) then - ! print *,'find: WARNING, found more than 1 value' - ! print *,'taking the first' - ! endif + ! if (iprt >= 3) then + ! print *,'find: WARNING, found more than 1 value' + ! print *,'taking the first' + ! endif !endif if (n > 1) then allocate(temp2(n))