module obs_sensitivity
!$$$ module documentation block
!           .      .    .                                       .
! module:   obs_sensitivity
!   prgmmr: tremolet
!
! abstract: Contains variables and routines for computation of
!           forecast sensitivity to observations.
!
! program history log:
!   2007-06-26 tremolet
!   2007-07-19 tremolet - increment sensitivity to observations
!   2009-08-07 lueken   - updated documentation
!   2010-04-30 tangborn - add pointer to carbon monoxide
!   2010-05-27 todling  - remove all user-specific TL-related references
!   2010-07-16 todling  - add reference to aero and aerol
!   2010-08-19 lueken   - add only to module use;no machine code, so use .f90
!   2011-03-29 todling  - add reference to pm2_5
!   2012-04-15 todling  - add reference to gust, vis, pblh
!   2015-07-10 pondeca  - add reference to wspd10m, td2m ,mxtm ,mitm ,pmsl,
!                         howv ,tcamt, lcbas, cldch
!   2016-02-20 pagowski - add pm10
!   2016-05-05 pondeca  - add reference to uwnd10m, vwnd10m
!   2017-05-12 Y. Wang and X. Wang - add reflectivity (dBZ), POC: xuguang.wang@ou.edu
!   2017-01-16 Apodaca  - add reference to lightning
!   2017-05-06 todling  - reload ensemble when FSOI calc doing EnVarFSOI-like
!   2017-05-21 todling  - implement ability to do 2nd EnVarFSOI-like calc
!
! Subroutines Included:
!   init_fc_sens  - Initialize computations
!
! Variable Definitions:
!   fcsens - forecast sensitivity gradient
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
! ------------------------------------------------------------------------------
use kinds, only: r_kind,i_kind,r_quad
use constants, only: zero, zero_quad, two
use gsi_4dvar, only: nobs_bins, l4dvar, lsqrtb, nsubwin
use gsi_4dvar, only: tau_fcst,evfsoi_order,evfsoi_afcst,evfsoi_ana
use jfunc, only: jiter, miter, niter, iter, jiterend

use gsi_obOperTypeManager, only: nobs_type => obOper_count
use gsi_obOperTypeManager, only: obOper_typeinfo
use mpimod, only: mype
use control_vectors, only: control_vector,allocate_cv,read_cv,deallocate_cv, &
    dot_product,assignment(=)
use state_vectors, only: allocate_state,deallocate_state
use gsi_bundlemod, only: self_add
use gsi_bundlemod, only: assignment(=)
use gsi_bundlemod, only: gsi_bundle
use bias_predictors, only: predictors,allocate_preds,deallocate_preds, &
    assignment(=)
use mpl_allreducemod, only: mpl_allreduce
use gsi_4dcouplermod, only: gsi_4dcoupler_getpert
use hybrid_ensemble_parameters,only : l_hyb_ens,ntlevs_ens,destroy_hybens_localization_parameters
use hybrid_ensemble_isotropic, only: create_ensemble,load_ensemble,destroy_ensemble
use hybrid_ensemble_isotropic, only: hybens_localization_setup
use mpeu_util, only: perr,die
! ------------------------------------------------------------------------------
implicit none
save
private
public lobsensfc,lobsensjb,lobsensincr,lobsensadj,&
       lobsensmin,iobsconv,llancdone,lsensrecompute, &
       fcsens, sensincr, &
       init_obsens, init_fc_sens, save_fc_sens, dot_prod_obs
public evfsoi_o2_update

public:: obsensCounts_realloc
public:: obsensCounts_set
public:: obsensCounts_dealloc

logical lobsensfc,lobsensjb,lobsensincr, &
        lobsensadj,lobsensmin,llancdone,lsensrecompute
integer(i_kind) :: iobsconv

! ------------------------------------------------------------------------------
type(control_vector) :: fcsens
real(r_kind), allocatable :: sensincr(:,:,:)
integer(i_kind):: orig_tau_fcst=-1

integer(i_kind),save,allocatable:: obscounts(:,:)

character(len=*),parameter:: myname="obs_sensitivity"
! ------------------------------------------------------------------------------
contains
!>> object obsensCounts_:
!>> this object was public obsmod::obscounts(:,:), but now private module
!>> variable in this module.  It is accessed through following module procedures
!>> []_alloc(), []_set() and []_dealloc().

subroutine obsensCounts_realloc(ntype,nbin)
!>> was implemented in setuprhsall()
  implicit none
  integer(i_kind),intent(in):: ntype
  integer(i_kind),intent(in):: nbin
  character(len=*),parameter:: myname_=myname//"::obsensCounts_realloc"
  if(allocated(obscounts)) deallocate(obscounts)
  allocate(obscounts(ntype,nbin))
end subroutine obsensCounts_realloc

subroutine obsensCounts_set(iobsglb)
!>> was implemented in evaljo()
  implicit none
  integer(i_kind),dimension(:,:),intent(in):: iobsglb
  character(len=*),parameter:: myname_=myname//"::obsensCounts_set"
  if(.not.allocated(obscounts)) then
    call perr(myname_,'not allocated, obscounts')
    call perr(myname_,'was evaljo() exception 125')
    call  die(myname_)
  endif
  if(any(shape(obscounts)/=shape(iobsglb))) then
    call perr(myname_,'mismatched, storage size(obscounts,1) =',size(obscounts,1))
    call perr(myname_,'             argument size(iobsglb,1) =',size(iobsglb,1))
    call perr(myname_,'            storage size(obscounts,2) =',size(obscounts,2))
    call perr(myname_,'             argument size(iobsglb,2) =',size(iobsglb,2))
    call  die(myname_)
  endif
  obscounts(:,:)=iobsglb(:,:)
end subroutine obsensCounts_set

subroutine obsensCounts_dealloc()
!>> was a part of obsmod::destroyobs_().
  implicit none
  character(len=*),parameter:: myname_=myname//"::obsensCounts_dealloc"
  if(allocated(obscounts)) deallocate(obscounts)
end subroutine obsensCounts_dealloc

! ------------------------------------------------------------------------------
subroutine init_obsens
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_obsens
!   prgmmr:
!
! abstract:
!
! program history log:
!   2009-08-07  lueken - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
implicit none

lobsensfc=.false.
lobsensjb=.false.
lobsensincr=.false.
lobsensadj=.false.
lobsensmin=.false.
lsensrecompute=.false.
llancdone=.false.
iobsconv=0

end subroutine init_obsens
! ------------------------------------------------------------------------------
subroutine init_fc_sens
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_fc_sens
!   prgmmr:      tremolet
!
! abstract: Read forecast sensitivity gradient
!
! program history log:
!   2007-06-26  tremolet - initial code
!   2009-08-07  lueken - added subprogram doc block
!   2010-05-27  todling - gsi_4dcoupler; remove dependence on GMAO specifics
!   2012-05-22  todling - update interface to getpert
!   2015-12-01  todling - add several obs-types that Pondeca forget to add here
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

implicit none

character(len=12) :: clfile
type(gsi_bundle) :: fcgrad(nsubwin)
type(gsi_bundle) :: eval(ntlevs_ens)
type(predictors) :: zbias
type(control_vector) :: xwork
real(r_kind) :: zjx
integer(i_kind) :: ii
character(len=80),allocatable,dimension(:)::fname
integer :: mydate(5)

if (mype==0) then
   write(6,*)'init_fc_sens: lobsensincr,lobsensfc,lobsensjb=', &
                            lobsensincr,lobsensfc,lobsensjb
   write(6,*)'init_fc_sens: lobsensadj,lobsensmin,iobsconv=', &
                            lobsensadj,lobsensmin,iobsconv
   write(6,*)'init_fc_sens: lsensrecompute=',lsensrecompute
endif

call allocate_cv(fcsens)
fcsens=zero

if (lobsensadj.and.lobsensmin) then
   write(6,*)'init_fc_sens: unknown method',lobsensadj,lobsensmin
   call stop2(155)
end if

if (iobsconv>=2) then
   allocate(sensincr(nobs_bins,nobs_type,niter(jiter)))
else
   allocate(sensincr(nobs_bins,nobs_type,1))
endif
sensincr=zero

! Initialize fcsens
if (lobsensfc) then
   if (lobsensincr) then
      clfile='xhatsave.ZZZ'
      write(clfile(10:12),'(I3.3)') jiter
      call read_cv(fcsens,clfile)
      if (jiter>1) then
         clfile='xhatsave.ZZZ'
         write(clfile(10:12),'(I3.3)') jiter-1
         call allocate_cv(xwork)
         call read_cv(xwork,clfile)
         do ii=1,fcsens%lencv
            fcsens%values(ii) = fcsens%values(ii) - xwork%values(ii)
         end do
         call deallocate_cv(xwork)
      endif
   else
      if (jiter==miter) then
         if (lobsensjb) then
            clfile='xhatsave.ZZZ'
            write(clfile(10:12),'(I3.3)') miter
            call read_cv(fcsens,clfile)
         else
!           read and convert output of GCM adjoint
            allocate(fname(nsubwin))
            fname='NULL'
            if (tau_fcst>0) then
               fname='B'
            endif
            do ii=1,nsubwin
               call allocate_state(fcgrad(ii))
            end do
            call allocate_preds(zbias)
            zbias=zero
            call gsi_4dcoupler_getpert(fcgrad,nsubwin,'adm',fname)
            if (tau_fcst>0) then ! .and. miter/=jiterend) then ! note: when jiterend=miter forecast ens already loaded
               mydate(1)=1776;mydate(2)=7;mydate(3)=4;mydate(4)=12;mydate(5)=0
               call advect_cv(mype,mydate,tau_fcst,fcgrad)
            endif
            if (lsqrtb) then
               call control2model_ad(fcgrad,zbias,fcsens)
            else
               if (l_hyb_ens) then
                  if (tau_fcst>0) then ! .and. miter/=jiterend) then ! note: when jiterend=miter forecast ens already loaded
                     call destroy_hybens_localization_parameters
                     call destroy_ensemble
                     call create_ensemble
                     ! now restore actual ens of backgrounds
                     call load_ensemble(tau_fcst)
                     call hybens_localization_setup
                  endif
                  do ii=1,ntlevs_ens
                     call allocate_state(eval(ii))
                  end do
                  eval(1)=fcgrad(1)
                  fcgrad(1)=zero
                  call ensctl2state_ad(eval,fcgrad(1),fcsens)
                  call control2state_ad(fcgrad,zbias,fcsens)
                  do ii=1,ntlevs_ens
                     call deallocate_state(eval(ii))
                  end do
                  if (tau_fcst>0) then
!                    mydate(1)=1776;mydate(2)=7;mydate(3)=4;mydate(4)=12;mydate(5)=0
!                    call advect_cv(mype,mydate,tau_fcst,fcsens)
                     call destroy_hybens_localization_parameters
                     call destroy_ensemble
                     call create_ensemble
                     ! now restore actual ens of backgrounds
                     evfsoi_afcst=.false.
                     orig_tau_fcst=tau_fcst
                     tau_fcst=-1
                     call load_ensemble(tau_fcst)
                     call hybens_localization_setup
                  endif
               else
                  call control2state_ad(fcgrad,zbias,fcsens)
               end if
            endif
            do ii=1,nsubwin
               call deallocate_state(fcgrad(ii))
            end do
            call deallocate_preds(zbias)
            deallocate(fname)
         endif
      else
!        read gradient from outer loop jiter+1
         clfile='fgsens.ZZZ'
         WRITE(clfile(8:10),'(I3.3)') jiter+1
         call read_cv(fcsens,clfile)
      endif
   endif
   zjx=dot_product(fcsens,fcsens)
   if (mype==0) write(6,888)'init_fc_sens: Norm fcsens=',sqrt(zjx)
endif
888 format(A,3(1X,ES25.18))

return
end subroutine init_fc_sens
! ------------------------------------------------------------------------------
subroutine evfsoi_o2_update(sval)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    order2_fc_sens
!   prgmmr:      todling
!
! abstract: Read forecast sensitivity gradient
!
! program history log:
!   2007-06-26  tremolet - initial code
!   2009-08-07  lueken - added subprogram doc block
!   2010-05-27  todling - gsi_4dcoupler; remove dependence on GMAO specifics
!   2012-05-22  todling - update interface to getpert
!   2015-12-01  todling - add several obs-types that Pondeca forget to add here
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

implicit none
type(gsi_bundle),intent(inout)  :: sval(nobs_bins)

character(len=80),allocatable,dimension(:)::fname
type(gsi_bundle) :: fcgrad(nsubwin)
type(gsi_bundle) :: eval(ntlevs_ens)
type(gsi_bundle) :: mval(nsubwin)
type(predictors) :: zbias
real(r_kind) :: zjx
integer(i_kind) :: ii

tau_fcst=orig_tau_fcst
if (tau_fcst<=0) return
if (evfsoi_order/=2) return
if (.not.l_hyb_ens) return

if (mype==0) then
   print *, 'in order2_fc_sens ...'
endif

if (lobsensadj.and.lobsensmin) then
   write(6,*)'order2_fc_sens: unknown method',lobsensadj,lobsensmin
   call stop2(155)
end if

! Zero out whatever is in fcsens
 fcsens=zero

! read and convert output of GCM adjoint
 allocate(fname(nsubwin))
 fname='A'
 do ii=1,nsubwin
    call allocate_state(fcgrad(ii))
 end do
 call allocate_preds(zbias)
 zbias=zero
 do ii=1,ntlevs_ens
    call allocate_state(eval(ii))
 end do
 do ii=1,nsubwin
    fcgrad(ii)=zero
 end do
 call gsi_4dcoupler_getpert(fcgrad,nsubwin,'adm',fname)
! Wipe out loaded ensemble members to allow reloading ensemble of analysis
 call destroy_ensemble
! Read in ens forecasts issued from "analysis"
 evfsoi_afcst=.true.
 call create_ensemble
 call load_ensemble(tau_fcst)
 do ii=1,ntlevs_ens
    call allocate_state(eval(ii))
 end do
 eval(1)=fcgrad(1)
 fcgrad(1)=zero
 call ensctl2state_ad(eval,fcgrad(1),fcsens)
 call control2state_ad(fcgrad,zbias,fcsens)
 do ii=1,ntlevs_ens
    call deallocate_state(eval(ii))
 end do
 do ii=1,nsubwin
    call deallocate_state(fcgrad(ii))
 end do
 call deallocate_preds(zbias)
 deallocate(fname)
! Wipeout ensemble from memory
 call destroy_ensemble

! RT: I think somewhere around here is where we should advect the resulting
! sensitivities ... TBD (similar to advecting localization scales, is it?)
! could also be a few lines above ... not sure yet

! Report magnitude of input vector
 zjx=dot_product(fcsens,fcsens)
 if (mype==0) write(6,888)'order2_fc_sens: Norm fcsens=',sqrt(zjx)
888 format(A,3(1X,ES25.18))

! Read in ensemble of analysis (these are EnKF, not GSI, analyses obviously)
 evfsoi_ana=.true.
 call create_ensemble
 call load_ensemble(-1)

! Allocate local variables
 do ii=1,nsubwin
    call allocate_state(mval(ii))
 end do
 do ii=1,ntlevs_ens
    call allocate_state(eval(ii))
 end do
 call allocate_preds(zbias)

! Convert from control variable to state space
 call control2state(fcsens,mval,zbias)

! Convert ensemble control variable to state space and update from previous value
 call ensctl2state(fcsens,mval(1),eval)
 do ii=1,ntlevs_ens
    call self_add(sval(ii),eval(ii))
 end do

! Release memory
call deallocate_preds(zbias)
do ii=1,ntlevs_ens
   call deallocate_state(eval(ii))
end do
do ii=1,nsubwin
   call deallocate_state(mval(ii))
end do

return
end subroutine evfsoi_o2_update
subroutine save_fc_sens
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    save_fc_sens
!   prgmmr:      tremolet
!
! abstract: Compute and save forecast sensitivity to observations
!
! program history log:
!   2007-06-26  tremolet - initial code
!   2009-08-07  lueken - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
implicit none

real(r_kind) :: zz
integer(i_kind) :: ii,jj,kk

! Save statistics
if (mype==0) then

!  Full stats
   do jj=1,nobs_type
      write(6,'(A,2X,I3,2X,A)')'Obs types:',jj,obOper_typeinfo(jj)
   enddo
   write(6,'(A,2X,I4)')'Obs bins:',nobs_bins
   write(6,*)'Obs Count Begin'
   if (.not.allocated(obscounts)) then
      write(6,*)'save_fc_sens: obscounts not allocated'
      call stop2(156)
   end if
   do jj=1,nobs_type
      write(6,'((1X,I12))')(obscounts(jj,ii),ii=1,nobs_bins)
   enddo
   write(6,*)'Obs Count End'

   write(6,*)'Obs Impact Begin'
   do kk=1,SIZE(sensincr,3)
      if (SIZE(sensincr,3)==1) then
         write(6,'(A,I4)')'Obs Impact iteration= ',niter(jiter)
      else
         write(6,'(A,I4)')'Obs Impact iteration= ',kk
      endif
      do jj=1,nobs_type
         write(6,'((1X,ES12.5))')(sensincr(ii,jj,kk),ii=1,nobs_bins)
      enddo
   enddo
   write(6,*)'Obs Impact End'

   kk=SIZE(sensincr,3)
!  Summary by obs type
   do jj=1,nobs_type
      zz=zero
      do ii=1,nobs_bins
         zz=zz+sensincr(ii,jj,kk)
      enddo
      if (zz/=zero) write(6,'(A,2X,A3,2X,ES12.5)')'Obs Impact type',obOper_typeinfo(jj),zz
   enddo

!  Summary by obs bins
   do ii=1,nobs_bins
      zz=zero
      do jj=1,nobs_type
         zz=zz+sensincr(ii,jj,kk)
      enddo
      if (zz/=zero) write(6,'(A,2X,I3,2X,ES12.5)')'Obs Impact bin',ii,zz
   enddo

endif

deallocate(sensincr)
call deallocate_cv(fcsens)

return
end subroutine save_fc_sens
! ------------------------------------------------------------------------------
real(r_kind) function dot_prod_obs()
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_fc_sens
!   prgmmr:      tremolet
!
! abstract: Computes scalar product in observation space
!           (based on evaljo)
!
! program history log:
!   2007-06-27  tremolet
!   2009-01-18  todling - carry summations in quad precision
!   2009-08-07  lueken - added subprogram doc block
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block
use m_obsdiagNode, only: obs_diag
use m_obsdiags, only: obsdiags
implicit none

integer(i_kind) :: ii,jj,ij,it
real(r_quad)    :: zzz
real(r_quad)    :: zprods(nobs_type*nobs_bins)
type(obs_diag),pointer:: obsptr
! ----------------------------------------------------------

zprods(:)=zero_quad

ij=0
do ii=1,size(obsdiags,2)
   do jj=1,size(obsdiags,1)
      ij=ij+1

      obsptr => obsdiags(jj,ii)%head
      do while (associated(obsptr))
         if (obsptr%luse.and.obsptr%muse(jiter)) then
            zprods(ij) = zprods(ij) + obsptr%nldepart(jiter) * obsptr%obssen(jiter)
         endif
         obsptr => obsptr%next
      enddo

   enddo
enddo

! Gather contributions
call mpl_allreduce(nobs_type*nobs_bins,qpvals=zprods)

! Save intermediate values
it=-1
if (iobsconv>=2) then
   if (iter>=1.and.iter<=niter(jiter)) it=iter
else
   it=1
endif

if (it>0) then
   ij=0
   do ii=1,nobs_bins
      do jj=1,nobs_type
         ij=ij+1
         sensincr(ii,jj,it)=zprods(ij)
      enddo
   enddo
endif

! Sum
zzz=zero_quad

ij=0
do ii=1,nobs_bins
   do jj=1,nobs_type
      ij=ij+1
      zzz=zzz+zprods(ij)
   enddo
enddo

dot_prod_obs=zzz

return
end function dot_prod_obs
! ------------------------------------------------------------------------------
end module obs_sensitivity


text:
      subroutine srchne(nrn1a,nrn2a,ns,nsta_asv,unam,uspeca)
c
c     This subroutine matches the species name unam with the
c     corresponding entry in the nrn1a-th through nrn2a-th range of
c     the species name array uspeca. Only the first 24 characters are
c     compared (uspeca has 48 characters, unam only 24). This
c     subroutine returns the species index ns. If there is no match,
c     the species name is converted to all lower case and the code again
c     searches for a match. If none is found, the species name is
c     coverted to all upper case and the code again searches for a
c     match. If no match is found, ns is returned with a value of 0.
c
c     This subroutine is much like EQLIB/srchn.f, except that if no
c     match is found, the above-noted case conversions are tried in
c     an effort to find a match.
c
c     This subroutine is called by:
c
c       None
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       nrn1a  = start of the range of species to search
c       nrn2a  = end of the range of species to search
c       unam   = name of the species whose index is to be found
c       uspeca = array of species names
c
c     Principal output:
c
c       ns     = index of the species whose name is unam
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nsta_asv
c
      integer nrn1a,nrn2a,ns
c
      character*48 uspeca(nsta_asv)
      character*24 unam
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer k,nfa,nfacap,nfdif,nfz,nfzcap,nlettr
c
      character*24 unamlc,unamuc
      character*1 ulettr
c
c-----------------------------------------------------------------------
c
c     Determine the ranges of upper and lower case letters.
c
      nfa = ichar('a')
      nfz = ichar('z')
      nfacap = ichar('A')
      nfzcap = ichar('Z')
      nfdif = ichar('A') - ichar('a')
c
c     Try the species name as is.
c
      call srchn(nrn1a,nrn2a,ns,nsta_asv,unam,uspeca)
c
      if (ns .le. 0) then
c
c       Convert to lower case.
c
        unamlc = unam
        do k = 1,24
          ulettr = unam(k:k)
          nlettr = ichar(ulettr)
          if (nlettr.ge.nfacap .and. nlettr.le.nfzcap)
     $    unamlc(k:k) = char(nlettr - nfdif)
        enddo
c
c       Test to see if this is different. If so, try it.
c
        if (unamlc .ne. unam) then
c
c         Calling sequence substitutions:
c           unamlc for unam
c
          call srchn(nrn1a,nrn2a,ns,nsta_asv,unamlc,uspeca)
        endif
      endif
c
      if (ns .le. 0) then
c
c       Convert to upper case.
c
        unamuc = unam
        do k = 1,24
          ulettr = unam(k:k)
          nlettr = ichar(ulettr)
          if (nlettr.ge.nfa .and. nlettr.le.nfz)
     $    unamuc(k:k) = char(nlettr + nfdif)
        enddo
c
c       Test to see if this is different. If so, try it.
c
        if (unamuc. ne. unam) then
c
c         Calling sequence substitutions:
c           unamuc for unam
c
          call srchn(nrn1a,nrn2a,ns,nsta_asv,unamuc,uspeca)
        endif
      endif
c
      end

ext: f


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
!  'partisn' problem
!  use of a module with only clause for name that is both
!  a module variable and a member in the module
	module part
	integer :: ib           ! ib is a module variable
	type btype
	   integer ib           ! ib is also a member
	endtype
	type(btype)::bb
	end module

	subroutine sub1
	use part, only: ib       ! ib appears in the only clause
	implicit none
	ib = 99                 ! warning & we create a local ib
	end

	subroutine sub2
	use part, only: ib=>bb, jb=>ib       ! ib appears in the only clause
	ib%ib = jb
	end

	program p
	use part
	integer result(2),expect(2)
	data expect/99,99/
	bb%ib = 0
	ib = 0
	call sub1
	call sub2
!	print *,ib
!	print *,bb%ib
	result(1) = ib
	result(2) = bb%ib
	call check(result,expect,2)
	end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
subroutine lmn_turnof(a)
   ! DESCRIPTION
   !    This routine closes the run for the incompressible NS equations
   !-----------------------------------------------------------------------
   use typre
   use Mod_iofile
   use Mod_LowMach
   implicit none
   class(LowMachProblem) :: a

   integer(ip) :: ndime,npoin,nelem,ncomp,ncsgs,pgaus,iaux,ielem,pnode,nboun,nbody,vesgs_coun,tesgs_coun,prsgs_coun,resid_coun,ibody

   !Unknowns
   call a%Mesh%GetNdime(ndime)
   call a%Mesh%GetNpoin(npoin)
   call a%Mesh%GetNelem(nelem)
   call a%Mesh%GetNboun(nboun)
   call a%Mesh%GetNbody(nbody)
   ncomp = a%ncomp

   !Write tail for formatted files
   if (a%MPIrank == a%MPIroot) then
      write(a%lun_outph,'(//,a,/)') '     * * * END OF LOW MACH RUN * * *'
      write(a%lun_solve,'(//,a,/)') '     * * * END OF LOW MACH RUN * * *'
   endif
   call a%Memor%dealloc(ndime,npoin,a%ncomp,a%veloc,'veloc','lmn_turnof')
   call a%Memor%dealloc(npoin,a%ncomp      ,a%press,'press','lmn_turnof')
   call a%Memor%dealloc(      a%ncomp      ,a%pther,'pther','lmn_turnof')
   call a%Memor%dealloc(npoin,a%ncomp      ,a%tempe,'tempe','lmn_turnof')
   call a%Memor%dealloc(npoin              ,a%itemp,'itemp','lmn_turnof')
   if(a%npp_stepi(4)>0) then
      call a%Memor%dealloc(npoin              ,a%densf,'densf','lmn_turnof')
   end if

   if (a%kfl_sourc == 2) then
      call a%Memor%dealloc(npoin,a%PointwiseSource,'PointwiseSource','lmn_reabcs')
   endif

   !Transient subgrid scales
   if(a%kfl_trasg/=0) then
      ncsgs=1
      vesgs_coun = 0
      tesgs_coun = 0
      prsgs_coun = 0
      if(a%kfl_tacsg>0) ncsgs=a%ncomp-1
      do ielem=1,nelem
         call a%Mesh%GetElemArraySize(ielem,pnode,pgaus)
         deallocate(a%vesgs(ielem)%a)
         deallocate(a%tesgs(ielem)%a)
         deallocate(a%prsgs(ielem)%a)
         vesgs_coun = vesgs_coun + ndime*ncsgs*pgaus
         tesgs_coun = tesgs_coun + ncsgs*pgaus
         prsgs_coun = prsgs_coun + pgaus
      end do
      call a%Memor%deallocObj(0,'vesgs%a','lmn_turnof',vesgs_coun*rp)
      call a%Memor%deallocObj(0,'tesgs%a','lmn_turnof',tesgs_coun*rp)
      call a%Memor%deallocObj(0,'prsgs%a','lmn_turnof',prsgs_coun*rp)
      call a%Memor%dealloc(nelem,a%vesgs,'vesgs','lmn_turnof')
      call a%Memor%dealloc(nelem,a%tesgs,'tesgs','lmn_turnof')
      call a%Memor%dealloc(nelem,a%prsgs,'prsgs','lmn_turnof')
   end if

   !Residual Projection
   if (a%kfl_repro >=1) then
      iaux = max(ndime+2,a%kfl_repro*(ndime+2))
      call a%Memor%dealloc(iaux,npoin,a%repro,'repro','lmn_turnof')
   endif

   !Gradient Orthogonal projection
   if(a%kfl_adapsgs == 1) then
      call a%Memor%dealloc(ndime+1,ndime,npoin,a%grprj,'grprj','lmn_turnof')
   end if

   !Postprocess the residual at the gauss points
   if (a%npp_stepi(11) /= 0) then
      resid_coun = 0
      do ielem=1,nelem
         call a%Mesh%GetElemArraySize(ielem,pnode,pgaus)
         deallocate(a%residualU(ielem)%a)
         deallocate(a%residualP(ielem)%a)
         deallocate(a%residualT(ielem)%a)
         resid_coun = resid_coun + (ndime+1)*pgaus
      end do
      call a%Memor%deallocObj(0,'residual%a','lmn_turnof',resid_coun*rp)
      call a%Memor%dealloc(nelem,a%residualU,'residual','lmn_turnof')
      call a%Memor%dealloc(nelem,a%residualP,'residual','lmn_turnof')
      call a%Memor%dealloc(nelem,a%residualT,'residual','lmn_turnof')
   endif

   if(a%kfl_outfm==1)then
      call a%Memor%dealloc(ndime,nbody,a%force,'force','lmn_turnof')
      call a%Memor%dealloc(3,nbody,a%momen,'momen','lmn_turnof')
      call a%Memor%dealloc(nbody,a%heatf,'heatf','lmn_turnof')
      do ibody=1,nbody
         close(a%lun_force(ibody))
      end do
      call a%Memor%dealloc(nbody,a%lun_force,'lun_force','lmn_turnof')
   end if

   if (a%MPIrank == a%MPIroot) then
      call iofile(two,a%lun_trap,a%fil_trap,'LMN TRACKING OF POINTS')
   endif

end subroutine lmn_turnof

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
        SUBROUTINE ZRCCRE(CCRE,GMLAT,NLAT,GMLONG,NLONG,BRATES,NBRATE,
     $			  SCALE,PEDGRE,ISTAT)
*
*  Module number:
*
*  Module name: ZRCCRE
*
*  Keyphrase:
*  ----------
*       Index table CCRE containing predicted background count rates
*
*  Description:
*  ------------
*       This routine reads table CCRE and extracts the geomagnetic
*	latitudes and longitudes and predicted background rates.
*       Sorted lists of latitude, longitude, and predicted background rates
*	are returned.  In addition, scale factors which are stored as
*       table keywords are returned.  The scale factors are needed to
*       adjust appropriately the model for the specific GHRS detector.
*
*  FORTRAN name: ZRCCRE.for
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*       CCRE                    I       table containing CCRE coeffecients
*  Subroutines Called:
*  -------------------
*  CDBS
*       zmsput
*  SDAS:
*       uttopn, utpgti, utcfnd, utrgt*, uttclo
*
*  History:
*  --------
*  Version      Date        Author          Description
*       1       Oct 96  M. De La Pena   Copied and modified from YRCCS8.FOR
*                                       written by S. Hulbert for CALFOS
*       1.1     May 98  M. De La Pena   Modified column name for CCRE table,
*                                       DETECTOR -> FOS_DETECTOR.
*-------------------------------------------------------------------------------
*
* Inputs:
*
*       ccre - table name (character*64)
*
* Outputs:
*
*	gmlat  - pointer to vector of geomagnetic latitudes (integer)
*	nlat   - number of latitudes (integer)
*	gmlong - pointer to vector of geomagnetic longitudes (integer)
*	nlong  - number of longitudes (integer)
*	brates - pointer to vector of predicted background rates (integer)
*	nbrate - number of rates (integer)
*       scale  - scale factors for model to accommodate GHRS detectors (real)
*       istat  - error status (integer)
*
*-----------------------------------------------------------------------------
C
C PASSED PARAMETERS
C
      INTEGER      NLAT,NLONG,NBRATE,ISTAT
      INTEGER      GMLAT,GMLONG,BRATES
      REAL         SCALE(2)
      CHARACTER*64 CCRE
C------------------------------------------------------------------------------
C
C IRAF MEM COMMON
C
      LOGICAL          MEMB(1)
      INTEGER*2        MEMS(1)
      INTEGER*4        MEMI(1)
      INTEGER*4        MEML(1)
      REAL             MEMR(1)
      DOUBLE PRECISION MEMD(1)
      COMPLEX          MEMX(1)
      EQUIVALENCE (MEMB, MEMS, MEMI, MEML, MEMR, MEMD, MEMX)
      COMMON /MEM/ MEMD
C------------------------------------------------------------------------------
C
C ZMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C------------------------------------------------------------------------------
C
C LOCAL VARIABLES
C
      INTEGER TYINT
      PARAMETER (TYINT = 4)
      INTEGER TYDOUB
      PARAMETER (TYDOUB = 7)
C
C FILE I/O ACCESS MODES
C
      INTEGER   RDONLY
      PARAMETER (RDONLY = 1)
C
      INTEGER   TBNROW
      PARAMETER (TBNROW = 21)
C
      REAL         DEFSCL(2),XLAT,XLONG
      INTEGER      I,J,IDIN,COLIDS(4),ISTATS(4),NROWS,ROW,IROW
      INTEGER      INDEX,NB
      CHARACTER*5  DETECTOR
      CHARACTER*15 COLNAM(4)
      CHARACTER*68 PEDGRE,DESCRP
      CHARACTER*80 CONTXT
      LOGICAL      NULL,FOLD
      DATA         DEFSCL/0.93,1.67/
      DATA         COLNAM/'FOS_DETECTOR','GM_LONG','GM_LAT','BACK_RATE'/
C------------------------------------------------------------------------------
C
C Open table
C
        CALL UTTOPN(CCRE,RDONLY,IDIN,ISTAT)
        IF(ISTAT.NE.0)THEN
            CONTXT='ERROR opening CCRE table '//CCRE
            GO TO 998
        ENDIF
C
C Get PEDIGREE, DESCRIP, and SCALE header keywords
C NOTE: PEDIGREE and DESCRIP do not exist in older calibration files
C
        CALL UTHGTT(IDIN,'PEDIGREE',PEDGRE,ISTATS(1))
        CALL UTHGTT(IDIN,'DESCRIP',DESCRP,ISTATS(2))
        CALL UTHGTR(IDIN,'D1SCALE',SCALE(1),ISTATS(3))
        CALL UTHGTR(IDIN,'D2SCALE',SCALE(2),ISTATS(4))
        DO 10 I = 3, 4
            IF(ISTATS(I).NE.0)THEN
                J = I - 2
                WRITE(CONTXT,11) DEFSCL(J)
11              FORMAT('ERROR reading scaling factor keyword; ',
     *                 'keyword set to default value: ',f4.2)
                CALL ZMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
                SCALE(J) = DEFSCL(J)
            ENDIF
10      CONTINUE
C
C If the table contains dummy data, then don't bother reading it
C
	IF (PEDGRE(1:5).EQ.'DUMMY') GO TO 997
C
C Get number of rows
C
        CALL UTPGTI(IDIN,TBNROW,NROWS,ISTAT)
        IF(ISTAT.NE.0)THEN
            CONTXT='ERROR reading CCRE table '//CCRE
            GO TO 999
        ENDIF
	IF(NROWS.LT.2)THEN
            CONTXT='ERROR -- not enough rows in CCRE table '//CCRE
            GO TO 999
	ENDIF
C
C Get column ids
C
        CALL UTCFND(IDIN,COLNAM,4,COLIDS,ISTAT)
        IF(ISTAT.NE.0)THEN
            CONTXT='ERROR locating needed columns in CCRE table '//CCRE
            GO TO 999
	ENDIF
C
C Allocate memory for arrays
C
        CALL UDMGET (NROWS, TYDOUB, GMLAT, ISTATS(1))
        CALL UDMGET (NROWS, TYDOUB, GMLONG, ISTATS(2))
        CALL UDMGET (NROWS, TYDOUB, BRATES, ISTATS(3))
        CALL UDMGET (NROWS, TYINT, INDEX, ISTATS(4))
        DO 20 I = 1, 4
            IF(ISTATS(I).NE.0)THEN
                CONTXT='ERROR allocating dynamic memory'
                GO TO 999
            ENDIF
20      CONTINUE
C
C Sort the table on the first three columns
C
	DO 30 I=1,NROWS
	    MEMI(INDEX+I-1)=I
30	CONTINUE
	FOLD = .TRUE.
	CALL TBTSRT (IDIN, 3, COLIDS, FOLD, NROWS, MEMI(INDEX))
C
C Load the arrays
C
	NB=0
        DO 40 ROW=1,NROWS
	    IROW=MEMI(INDEX+ROW-1)
            CALL UTRGTT(IDIN,COLIDS(1),1,IROW,DETECTOR,NULL,ISTAT)
            IF(ISTAT.NE.0)THEN
                CONTXT='ERROR reading DETECTOR from CCRE table '//CCRE
                GO TO 999
            ENDIF
	    NB=NB+1
            CALL UTRGTD(IDIN,COLIDS(3),1,IROW,MEMD(GMLAT+NB-1),
     $			NULL,ISTAT)
            IF(ISTAT.NE.0)THEN
                CONTXT='ERROR reading GM_LAT from CCRE table '//CCRE
                GO TO 999
            ENDIF
            CALL UTRGTD(IDIN,COLIDS(2),1,IROW,MEMD(GMLONG+NB-1),
     $			NULL,ISTAT)
            IF(ISTAT.NE.0)THEN
                CONTXT='ERROR reading GM_LONG from CCRE table '//
     $				CCRE
                GO TO 999
            ENDIF
            CALL UTRGTD(IDIN,COLIDS(4),1,IROW,MEMD(BRATES+NB-1),
     $		        NULL,ISTAT)
            IF(ISTAT.NE.0)THEN
                CONTXT='ERROR reading BACK_RATE from CCRE table '//
     $				CCRE
                GO TO 999
            ENDIF
40      CONTINUE
	NBRATE=NB
C
C Sort for unique latitudes and longitudes and return just these values
C
	NLAT=1
	XLAT=MEMD(GMLAT+NLAT-1)
	DO 50 I=2,NBRATE
	    IF(MEMD(GMLAT+I-1).GT.XLAT)THEN
		NLAT=NLAT+1
		MEMD(GMLAT+NLAT-1)=MEMD(GMLAT+I-1)
		XLAT=MEMD(GMLAT+NLAT-1)
	    ENDIF
50	CONTINUE
	NLONG=1
	XLONG=MEMD(GMLONG+NLONG-1)
	DO 60 I=2,NBRATE
	    IF(MEMD(GMLONG+I-1).GT.XLONG)THEN
		NLONG=NLONG+1
		MEMD(GMLONG+NLONG-1)=MEMD(GMLONG+I-1)
		XLONG=MEMD(GMLONG+NLONG-1)
	    ENDIF
60	CONTINUE
C
        CALL UDMFRE (INDEX, TYINT, ISTAT)
        IF(ISTAT.NE.0)THEN
            CONTXT='ERROR allocating dynamic memory'
            GO TO 999
        ENDIF
C
997     CALL UTTCLO(IDIN,ISTAT)
        ISTAT=0
        GO TO 1000
999     CALL UTTCLO(IDIN,ISTAT)
998     CALL ZMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
        ISTAT=1
1000    RETURN
        END

max_stars_repo_name: iraf-community/stsdas


text:
      SUBROUTINE PLT_CNTR_RANGE
     &          (amp,mxlat,nrlat,mxlon,nrlon,
     &           ampinc,ampmin,ampmax)

c Determines the range of contour levels found in input array

c Change History

c*******************!***************************************************

      dimension     amp(mxlon,mxlat)

c     Get data range
      ampmin=amp(1,1)
      ampmax=amp(1,1)
      do l=2,nrlat
         do k=2,nrlon
            ampmin=AMIN1(ampmin,amp(k,l))
            ampmax=AMAX1(ampmax,amp(k,l))
         end do
      end do

      ampmin=INT(ampmin/ampinc-.9)*ampinc
      if (ampmin .gt. 0.) ampmin=ampmin+ampinc

      ampmax=INT(ampmax/ampinc+.9)*ampinc
      if (ampmax .lt. 0.) ampmax=ampmax-ampinc

      RETURN
      END      ! PLT_CNTR_RANGE

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


text:
C-------GF 28.08.1980-------------------------------- M E M S
      INTEGER*2 FMEM,MEMHIB
      INTEGER*2 EOP
      INTEGER*2 MEM(XXMEMH),MEMNUC(XXMEMH)
      COMMON /MEMS/
     =  FMEM      ! 1ST FREE MEMBER
     = ,MEMHIB
     = ,EOP      ! SYMBOL AT THE END OF A PRODUCTION
     = ,MEM      ! SYMBOL-NUMBER, =1: END OF PRODUCTION
      COMMON /MEMS/
     =  MEMNUC      ! ANCHOR FOR 'NUC'-LISTS

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


text:
I am charming,and i am very entertaing. I am new to Davis Wiki and I like it!!

welcome! Daubert

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! Copyright (c) 2015-2021, the ELSI team.
! All rights reserved.
!
! This file is part of ELSI and is distributed under the BSD 3-clause license,
! which may be found in the LICENSE file in the ELSI root directory.

!>
!! Interface to SIPS.
!!
module ELSI_SIPS

   use ELSI_CONSTANT, only: UNSET,GET_DM
   use ELSI_DATATYPE, only: elsi_param_t,elsi_basic_t
   use ELSI_MALLOC, only: elsi_allocate,elsi_deallocate
   use ELSI_MPI
   use ELSI_OUTPUT, only: elsi_say,elsi_get_time
   use ELSI_PRECISION, only: r8,i4
   use M_SIPS, only: sips_initialize,sips_load_ham_ovlp,sips_load_ham,&
       sips_update_ham,sips_set_eps,sips_update_eps,sips_set_slices,&
       sips_get_slices,sips_get_inertias,sips_get_slices_from_inertias,&
       sips_solve_eps,sips_get_eigenvalues,sips_get_eigenvectors,&
       sips_get_dm_edm,sips_finalize

   implicit none

   private

   public :: elsi_init_sips
   public :: elsi_cleanup_sips
   public :: elsi_solve_sips
   public :: elsi_build_dm_edm_sips

   interface elsi_solve_sips
      module procedure elsi_solve_sips_real
   end interface

   interface elsi_build_dm_edm_sips
      module procedure elsi_build_dm_edm_sips_real
   end interface

contains

!>
!! Initialize SLEPc-SIPs.
!!
subroutine elsi_init_sips(ph,bh)

   implicit none

   type(elsi_param_t), intent(inout) :: ph
   type(elsi_basic_t), intent(inout) :: bh

   character(len=*), parameter :: caller = "elsi_init_sips"

   if(.not. ph%sips_started) then
      call sips_initialize()

      if(ph%sips_n_slices == UNSET) then
         ! TODO: Number of slices
         ph%sips_n_slices = bh%n_procs
      end if

      ! 1D block distribution
      bh%n_lcol_sp1 = ph%n_basis/bh%n_procs

      ! The last process holds all remaining columns
      if(bh%myid == bh%n_procs-1) then
         bh%n_lcol_sp1 = ph%n_basis-(bh%n_procs-1)*bh%n_lcol_sp1
      end if

      if(bh%n_lcol_sp == UNSET) then
         bh%n_lcol_sp = bh%n_lcol_sp1
      end if

      ph%sips_started = .true.
   end if

end subroutine

!>
!! Interface to SLEPc-SIPs.
!!
subroutine elsi_solve_sips_real(ph,bh,row_ind,col_ptr,ham,ovlp,eval,evec)

   implicit none

   type(elsi_param_t), intent(inout) :: ph
   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i4), intent(inout) :: row_ind(bh%nnz_l_sp1)
   integer(kind=i4), intent(inout) :: col_ptr(bh%n_lcol_sp1+1)
   real(kind=r8), intent(in) :: ham(bh%nnz_l_sp1)
   real(kind=r8), intent(in) :: ovlp(bh%nnz_l_sp1)
   real(kind=r8), intent(inout) :: eval(ph%n_states)
   real(kind=r8), intent(out) :: evec(bh%n_lcol_sp1,ph%n_states)

   real(kind=r8) :: t0
   real(kind=r8) :: t1
   real(kind=r8) :: diff
   integer(kind=i4) :: i_val
   integer(kind=i4) :: n_solved
   logical :: inertia_ok
   character(len=200) :: msg

   real(kind=r8), allocatable :: eval_save(:)
   real(kind=r8), allocatable :: slices(:)
   integer(kind=i4), allocatable :: inertias(:)

   character(len=*), parameter :: caller = "elsi_solve_sips_real"

   write(msg,"(A)") "Starting SLEPc-SIPs eigensolver"
   call elsi_say(bh,msg)

   call elsi_get_time(t0)

   if(ph%sips_first) then
      if(.not. ph%unit_ovlp) then
         ! Load H and S
         call sips_load_ham_ovlp(ph%n_basis,bh%n_lcol_sp1,bh%nnz_l_sp1,row_ind,&
              col_ptr,ham,ovlp)
         call sips_set_eps(0)
      else
         ! Load H
         call sips_load_ham(ph%n_basis,bh%n_lcol_sp1,bh%nnz_l_sp1,row_ind,&
              col_ptr,ham)
         call sips_set_eps(1)
      end if
   else
      ! Update H matrix
      call sips_update_ham(bh%n_lcol_sp1,bh%nnz_l_sp1,row_ind,col_ptr,ham)
      call sips_update_eps(ph%sips_n_slices)
   end if

   call elsi_get_time(t1)

   write(msg,"(A)") "Finished loading matrices"
   call elsi_say(bh,msg)
   write(msg,"(A,F10.3,A)") "| Time :",t1-t0," s"
   call elsi_say(bh,msg)

   call elsi_allocate(bh,slices,ph%sips_n_slices+1,"slices",caller)
   call elsi_allocate(bh,eval_save,ph%n_states,"eval_save",caller)

   eval_save(:) = eval

   if(ph%sips_first .or. ph%sips_do_inertia) then
      call elsi_get_time(t0)

      call elsi_allocate(bh,inertias,ph%sips_n_slices+1,"inertias",caller)

      inertia_ok = .false.
      eval(1) = eval(1)-ph%sips_buffer
      eval(ph%n_states) = eval(ph%n_states)+ph%sips_buffer

      do while(.not. inertia_ok)
         inertia_ok = .true.

         call sips_get_slices(0,ph%n_states,ph%sips_n_slices,0.0_r8,1.0e-5_r8,&
              eval,slices)
         call sips_get_inertias(ph%sips_n_slices,slices,inertias)

         if(inertias(1) > 0) then
            eval(1) = eval(1)-ph%sips_buffer
            inertia_ok = .false.
         else
            do i_val = 1,ph%sips_n_slices
               if(inertias(i_val+1) > 0) then
                  eval(1) = slices(i_val)

                  if(i_val > 1) then
                     inertia_ok = .false.
                  end if

                  exit
               end if
            end do
         end if

         if(inertias(ph%sips_n_slices+1) < ph%n_states) then
            eval(ph%n_states) = eval(ph%n_states)+ph%sips_buffer
            inertia_ok = .false.
         else
            do i_val = ph%sips_n_slices+1,2,-1
               if(inertias(i_val-1) < ph%n_states) then
                  eval(ph%n_states) = slices(i_val)

                  if(i_val < ph%sips_n_slices+1) then
                     inertia_ok = .false.
                  end if

                  exit
               end if
            end do
         end if
      end do

      call sips_get_slices_from_inertias(ph%n_states,ph%sips_n_slices,inertias,&
           slices)

      call elsi_deallocate(bh,inertias,"inertias")

      call elsi_get_time(t1)

      write(msg,"(A)") "Finished inertia counting"
      call elsi_say(bh,msg)
      write(msg,"(A,F10.3,A)") "| Time :",t1-t0," s"
      call elsi_say(bh,msg)
   else
      call sips_get_slices(2,ph%n_states,ph%sips_n_slices,&
           ph%sips_inertia_tol*2,1.0e-5_r8,eval,slices)
   end if

   call sips_set_slices(ph%sips_n_slices,slices)

   call elsi_get_time(t0)

   ! Solve
   call sips_solve_eps(n_solved)

   if(n_solved < ph%n_states) then
      write(msg,"(A)") "SLEPc-SIPs solver failed"
      call elsi_stop(bh,msg,caller)
   end if

   call elsi_get_time(t1)

   write(msg,"(A)") "Finished solving generalized eigenproblem"
   call elsi_say(bh,msg)
   write(msg,"(A,F10.3,A)") "| Time :",t1-t0," s"
   call elsi_say(bh,msg)

   call elsi_get_time(t0)

   ! Get solutions
   call sips_get_eigenvalues(ph%n_states,eval(1:ph%n_states))
   call sips_get_eigenvectors(ph%n_states,bh%n_lcol_sp1,evec)

   diff = maxval(abs(eval_save-eval))

   if(diff > ph%sips_inertia_tol) then
      ph%sips_do_inertia = .true.
   else
      ph%sips_do_inertia = .false.
   end if

   call elsi_deallocate(bh,eval_save,"eval_save")
   call elsi_deallocate(bh,slices,"slices")

   call elsi_get_time(t1)

   write(msg,"(A)") "Finished retrieving eigensolutions"
   call elsi_say(bh,msg)
   write(msg,"(A,F10.3,A)") "| Time :",t1-t0," s"
   call elsi_say(bh,msg)

   ph%sips_first = .false.

end subroutine

!>
!! Construct density matrix or energy-weighted density matrix from eigenvectors.
!!
subroutine elsi_build_dm_edm_sips_real(ph,bh,row_ind,col_ptr,occ,dm,which)

   implicit none

   type(elsi_param_t), intent(in) :: ph
   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i4), intent(inout) :: row_ind(bh%nnz_l_sp1)
   integer(kind=i4), intent(inout) :: col_ptr(bh%n_lcol_sp1+1)
   real(kind=r8), intent(in) :: occ(ph%n_states)
   real(kind=r8), intent(out) :: dm(bh%nnz_l_sp1)
   integer(kind=i4), intent(in) :: which

   real(kind=r8) :: t0
   real(kind=r8) :: t1
   character(len=200) :: msg

   character(len=*), parameter :: caller = "elsi_build_dm_edm_sips_real"

   call elsi_get_time(t0)

   call sips_get_dm_edm(bh%n_lcol_sp1,bh%nnz_l_sp1,row_ind,col_ptr,ph%n_states,&
        occ,dm,which)

   call elsi_get_time(t1)

   if(which == GET_DM) then
      write(msg,"(A)") "Finished density matrix calculation"
   else
      write(msg,"(A)") "Finished energy density matrix calculation"
   end if

   call elsi_say(bh,msg)
   write(msg,"(A,F10.3,A)") "| Time :",t1-t0," s"
   call elsi_say(bh,msg)

end subroutine

!>
!! Clean up SLEPc-SIPs.
!!
subroutine elsi_cleanup_sips(ph)

   implicit none

   type(elsi_param_t), intent(inout) :: ph

   character(len=*), parameter :: caller = "elsi_cleanup_sips"

   if(ph%sips_started) then
      call sips_finalize()
   end if

   ph%sips_first = .true.
   ph%sips_started = .false.

end subroutine

end module ELSI_SIPS


text:
! Program to check resolution of symbols with the same name
program contained2
   implicit none
   integer var1

   var1 = 42
   if (f1() .ne. 1) call abort
   call f2()
   if (var1 .ne. 42) call abort
contains

function f1 ()
   implicit none
   integer f1
   integer var1
   integer f2

   var1 = 1
   f2 = var1
   f1 = f2
end function

subroutine f2()
   implicit none
   if (f1() .ne. 1) call abort
end subroutine

end program

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
!
! ODZeeman_Predictor
!
! Module containing routines to compute the predictors for calculations
! gaseous absorption affected by Zeeman spilitting in
! the Optical Depth Pressure Space (ODPS).
!
!
! CREATION HISTORY:
!       Written by:     Yong Han, JCSDA, NOAA/NESDIS 20-Sep-2008
!

MODULE ODZeeman_Predictor

  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds,            ONLY: fp
  USE Message_Handler,       ONLY: SUCCESS, FAILURE, Display_Message
  USE ODPS_Predictor_Define, ONLY: ODPS_Predictor_type
  USE ODPS_Define,           ONLY: ODPS_type

  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  ! Everything private by default
  PRIVATE
  ! Public routines
  PUBLIC :: Compute_Predictors_zssmis
  PUBLIC :: Compute_Predictors_zssmis_TL
  PUBLIC :: Compute_Predictors_zssmis_AD
  PUBLIC :: Compute_Predictors_zamsua
  PUBLIC :: Compute_Predictors_zamsua_TL
  PUBLIC :: Compute_Predictors_zamsua_AD

  ! Public parameters
  PUBLIC :: N_ZCOMPONENTS
  PUBLIC :: N_ZABSORBERS
  PUBLIC :: MAX_N_PREDICTORS_ZSSMIS
  PUBLIC :: MAX_N_PREDICTORS_ZAMSUA
  PUBLIC :: ZSSMIS_ChannelMap
  PUBLIC :: ZAMSUA_ChannelMap
  PUBLIC :: ODPS_gINDEX_ZSSMIS
  PUBLIC :: ODPS_gINDEX_ZAMSUA

  ! ----------
  ! Parameters
  ! ----------
  REAL(fp), PARAMETER :: ZERO = 0.0_fp, ONE = 1.0_fp, TWO = 2.0_fp

  INTEGER, PARAMETER :: N_ZCOMPONENTS = 1
  ! set this variable to 1 although it should be 1 to consistant with the ODPS structure (for
  ! Zeeman applications only O2 need to be considered, which will be include as a fixed gas, decribed
  ! by pressure and temperature.
  INTEGER, PARAMETER :: N_ZABSORBERS  = 1

  !----------------------------------------------------------------
  !  ZSSMIS parameters (SSMIS channels 19 - 22)
  !----------------------------------------------------------------
  INTEGER, PARAMETER :: ODPS_gINDEX_ZSSMIS = 4  ! ODPS group index
  INTEGER, PARAMETER :: MAX_N_PREDICTORS_ZSSMIS = 18
  ! Global to ZSSMIS channel index mapping
  INTEGER, PARAMETER :: N_CHANNELS_SSMIS = 24
  INTEGER, PARAMETER :: ZSSMIS_ChannelMap(N_CHANNELS_SSMIS) = (/&
                                      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
                                      0, 0, 0, 0, 0, 0, 0, 0, 1, 2, &
                                      3, 4, 0, 0 /)

  !----------------------------------------------------------------
  !  ZAMSUA parameters (AMSUA channel 14)
  !----------------------------------------------------------------
  INTEGER, PARAMETER :: MAX_N_PREDICTORS_ZAMSUA = 7
  INTEGER, PARAMETER :: ODPS_gINDEX_ZAMSUA = 5  ! ODPS group index
  ! Global to ZAMSUA channel index mapping
  INTEGER, PARAMETER :: N_CHANNELS_AMSUA = 15
  INTEGER, PARAMETER :: ZAMSUA_ChannelMap(N_CHANNELS_AMSUA) = (/&
                                      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
                                      0, 0, 0, 1, 0/)

CONTAINS

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_Predictors_zssmis
!
! PURPOSE:
!       Subroutine to compute predictors for the ZSSMIS sensor (a virtual sensor
!       created for SSMIS Zeeman channels 19, 20, 21 and 22
!
! CALLING SEQUENCE:
!       CALL  Compute_Predictors_zssmis(Temperature,   &
!                                       Be,            &
!                                       CosBK,         &
!                                       Doppler_Shift, &
!                                       Secang,        &
!                                       Predictor )
!
! INPUT ARGUMENTS:
!
!       Temperature:     Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!         Be    :        Earth magnetic filed strength
!                        UNITS:      Guass (range 0.2 - 0.7)
!                        TYPE:       REAL(fp)
!                        DIMENSION:  scalar
!                        ATTRIBUTES: INTENT(IN)
!
!         CosBK :        Cosine of the angle between the Earth magnetic
!                        field and the electromagnetic wave propagation
!                        direction
!                        UNITS:      N/A
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Doppler_Shift :  Doppler shift due Earth ratation, positive when
!                        the emitter moves towards the sensor
!                        UNITS:      KHz
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!             Secang  :  Secant of the zenith angle profile
!                        UNITS:      N/A
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       Predictor:      Predictor structure
!                       UNITS:      N/A
!                       TYPE:       TYPE(ODPS_Predictor_type)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN OUT)
!
!------------------------------------------------------------------------------

  SUBROUTINE Compute_Predictors_zssmis(Temperature,   &
                                       Be,            &
                                       CosBK,         &
                                       Doppler_Shift, &
                                       Secang,        &
                                       Predictor )
    REAL(fp), INTENT( IN ) :: Temperature(:)
    REAL(fp), INTENT( IN ) :: Be
    REAL(fp), INTENT( IN ) :: CosBK
    REAL(fp), INTENT( IN ) :: Doppler_Shift
    REAL(fp), INTENT( IN ) :: Secang(:)
    TYPE(ODPS_Predictor_type), INTENT( INOUT ) :: Predictor

    ! LOCAL
    REAL(fp) :: Be2, Be3, OneOverBe, OneOverBe2, CosBK2, OneOverT
    INTEGER  :: k

    Be2         = Be  * Be
    Be3         = Be2 * Be
    OneOverBe   = ONE / Be
    OneOverBe2  = OneOverBe * OneOverBe
    CosBK2     = CosBK * CosBK

    DO k = 1, Predictor%n_Layers
      ! Predictors associated Temperature
      OneOverT = 300.0_fp / Temperature(k)
      Predictor%X(k, 1,  1)  = OneOverT                ! 1/T
      Predictor%X(k, 2,  1)  = OneOverT * CosBK2       ! COS(BK_ang)**2 / T
      Predictor%X(k, 3,  1)  = OneOverT / Be           ! 1 / (T*Be)
      Predictor%X(k, 4,  1)  = OneOverT * CosBK / Be   ! COS(BK_ang) / (T*Be)
      Predictor%X(k, 5,  1)  = OneOverT * CosBK / Be2  ! COS(BK_ang) / (T*Be**2)
    END DO

    ! Constant predictor
    Predictor%X(:, 6,  1)  = ONE                       ! 1  - constant term
    Predictor%X(:, 7,  1)  = CosBK                     ! COS(BK_ang)
    Predictor%X(:, 8,  1)  = CosBK2                    ! COS(BK_ang)**2
    Predictor%X(:, 9,  1)  = OneOverBe                 ! 1 / Be
    Predictor%X(:, 10, 1)  = OneOverBe2                ! 1 / Be**2
    Predictor%X(:, 11, 1)  = OneOverBe2 * OneOverBe    ! 1 / Be^3
    Predictor%X(:, 12, 1)  = CosBK * Be2               ! COS(BK_ang) * Be**2
    Predictor%X(:, 13, 1)  = CosBK2 * Be2              ! (COS(BK_ang) * Be)**2
    Predictor%X(:, 14, 1)  = CosBK * Be3               ! COS(BK_ang) * Be**3
    Predictor%X(:, 15, 1)  = Be                        ! Be
    Predictor%X(:, 16, 1)  = Be3                       ! Be**3
    Predictor%X(:, 17, 1)  = CosBK  * Be               ! COS(BK_ang)*Be
    Predictor%X(:, 18, 1)  = CosBK2 * Be               ! COS(BK_ang)**2*Be

    Predictor%n_CP(1) = MAX_N_PREDICTORS_ZSSMIS

    ! pass Doppler shift and Secant angle to Predicotr structre
    Predictor%u = Doppler_Shift
    Predictor%Secant_Zenith = Secang

  END SUBROUTINE Compute_Predictors_zssmis

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_Predictors_zssmis_TL
!
! PURPOSE:
!       Subroutine to compute TL predictors for the ZSSMIS sensor (a virtual sensor
!       created for SSMIS Zeeman channels 19, 20, 21 and 22)
!
! CALLING SEQUENCE:
!
!       CALL Compute_Predictors_zssmis_TL(Temperature,   &
!                                         Be,            &
!                                         CosBK,         &
!                                         Temperature_TL,&
!                                         Predictor_TL )
! INPUT ARGUMENTS:
!
!       Temperature:     Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!         Be    :        Earth magnetic filed strength
!                        UNITS:      Guass (range 0.2 - 0.7)
!                        TYPE:       REAL(fp)
!                        DIMENSION:  scalar
!                        ATTRIBUTES: INTENT(IN)
!
!         CosBK :        Cosine of the angle between the Earth magnetic
!                        field and the electromagnetic wave propagation
!                        direction
!                        UNITS:      N/A
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Temperature_TL:  TL Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       Predictor_TL:   TL Predictor structure
!                       UNITS:      N/A
!                       TYPE:       TYPE(ODPS_Predictor_type)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN OUT)
!
!------------------------------------------------------------------------------

  SUBROUTINE Compute_Predictors_zssmis_TL(Temperature,   &
                                          Be,            &
                                          CosBK,         &
                                          Temperature_TL,&
                                          Predictor_TL )
    REAL(fp), INTENT( IN ) :: Temperature(:)
    REAL(fp), INTENT( IN ) :: Be
    REAL(fp), INTENT( IN ) :: CosBK
    REAL(fp), INTENT( IN ) :: Temperature_TL(:)
    TYPE(ODPS_Predictor_type), INTENT( INOUT ) :: Predictor_TL

    ! LOCAL
    REAL(fp) :: Be2, CosBK2, OneOverT, OneOverT_TL
    INTEGER  :: k

    Be2         = Be  * Be
    CosBK2     = CosBK * CosBK

    DO k = 1, Predictor_TL%n_Layers
      ! Predictors associated Temperature
      OneOverT = 300.0_fp / Temperature(k)
      OneOverT_TL = -(OneOverT / Temperature(k))*Temperature_TL(k)
      Predictor_TL%X(k, 1,  1)  = OneOverT_TL
      Predictor_TL%X(k, 2,  1)  = OneOverT_TL * CosBK2
      Predictor_TL%X(k, 3,  1)  = OneOverT_TL / Be
      Predictor_TL%X(k, 4,  1)  = OneOverT_TL * CosBK / Be
      Predictor_TL%X(k, 5,  1)  = OneOverT_TL * CosBK / Be2
    END DO

    ! Constant predictor
    Predictor_TL%X(:, 6:18,  1)  = ZERO

    Predictor_TL%n_CP(1) = MAX_N_PREDICTORS_ZSSMIS

  END SUBROUTINE Compute_Predictors_zssmis_TL

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_Predictors_zssmis_AD
!
! PURPOSE:
!       Subroutine to compute AD predictors for the ZSSMIS sensor (a virtual sensor
!       created for SSMIS Zeeman channels 19, 20, 21 and 22)
!
! CALLING SEQUENCE:
!
!        CALL Compute_Predictors_zssmis_AD(Temperature,   &
!                                          Be,            &
!                                          CosBK,         &
!                                          Predictor_AD,  &
!                                          Temperature_AD )
! INPUT ARGUMENTS:
!
!       Temperature:     Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!         Be    :        Earth magnetic filed strength
!                        UNITS:      Guass (range 0.2 - 0.7)
!                        TYPE:       REAL(fp)
!                        DIMENSION:  scalar
!                        ATTRIBUTES: INTENT(IN)
!
!         CosBK :        Cosine of the angle between the Earth magnetic
!                        field and the electromagnetic wave propagation
!                        direction
!                        UNITS:      N/A
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!       Predictor_AD:    AD Predictor structure
!                        UNITS:      N/A
!                        TYPE:       TYPE(ODPS_Predictor_type)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(INOUT)
!
! OUTPUT ARGUMENTS:
!       Temperature_AD:  AD Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(INOUT)
!
!------------------------------------------------------------------------------
  SUBROUTINE Compute_Predictors_zssmis_AD(Temperature,   &
                                          Be,            &
                                          CosBK,         &
                                          Predictor_AD,  &
                                          Temperature_AD )
    REAL(fp), INTENT( IN ) :: Temperature(:)
    REAL(fp), INTENT( IN ) :: Be
    REAL(fp), INTENT( IN ) :: CosBK
    TYPE(ODPS_Predictor_type), INTENT( INOUT ) :: Predictor_AD
    REAL(fp), INTENT( INOUT ) :: Temperature_AD(:)

    ! LOCAL
    REAL(fp) :: Be2, CosBK2, OneOverT, OneOverT_AD
    INTEGER  :: k

    Be2         = Be  * Be
    CosBK2     = CosBK * CosBK

    DO k = 1, Predictor_AD%n_Layers
      ! Predictors associated Temperature
      OneOverT = 300.0_fp / Temperature(k)

      OneOverT_AD = Predictor_AD%X(k, 1,  1) + &
                    Predictor_AD%X(k, 2,  1)*CosBK2 + &
                    Predictor_AD%X(k, 3,  1) / Be + &
                    Predictor_AD%X(k, 4,  1) * CosBK / Be + &
                    Predictor_AD%X(k, 5,  1) * CosBK / Be2
      Temperature_AD(k) = Temperature_AD(k) - &
                     (OneOverT / Temperature(k)) * OneOverT_AD
    END DO

    ! Constant predictor
    Predictor_AD%X(:, :,  1)  = ZERO

    Predictor_AD%n_CP(1) = MAX_N_PREDICTORS_ZSSMIS

  END SUBROUTINE Compute_Predictors_zssmis_AD

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_Predictors_zamsua
!
! PURPOSE:
!       Subroutine to compute predictors for the ZAMSUA sensor (a virtual sensor
!       with AMSUA channel 14)
!
! CALLING SEQUENCE:
!       CALL  Compute_Predictors_zamsua(Temperature,     &
!                                       Ref_Temperature, &
!                                       Be,            &
!                                       CosBK,         &
!                                       Secang,        &
!                                       Predictor )
!
! INPUT ARGUMENTS:
!
!       Temperature:     Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!   Ref_Temperature:     Reference temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!         Be    :        Earth magnetic filed strength
!                        UNITS:      Guass (range 0.2 - 0.7)
!                        TYPE:       REAL(fp)
!                        DIMENSION:  scalar
!                        ATTRIBUTES: INTENT(IN)
!
!         CosBK :        Cosine of the angle between the Earth magnetic
!                        field and the electromagnetic wave propagation
!                        direction
!                        UNITS:      N/A
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(IN)
!
!             Secang  :  Secant of the zenith angle profile
!                        UNITS:      N/A
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       Predictor:      Predictor structure
!                       UNITS:      N/A
!                       TYPE:       TYPE(ODPS_Predictor_type)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN OUT)
!
!------------------------------------------------------------------------------

  SUBROUTINE Compute_Predictors_zamsua(Temperature,     &
                                       Ref_Temperature, &
                                       Be,              &
                                       CosBK,           &
                                       Secang,          &
                                       Predictor )
    REAL(fp), INTENT( IN ) :: Temperature(:)
    REAL(fp), INTENT( IN ) :: Ref_Temperature(:)
    REAL(fp), INTENT( IN ) :: Be
    REAL(fp), INTENT( IN ) :: CosBK
    REAL(fp), INTENT( IN ) :: Secang(:)
    TYPE(ODPS_Predictor_type), INTENT( INOUT ) :: Predictor

    ! LOCAL
    REAL(fp) :: Be3,  CosBK2, CosBK2Be2, tk
    INTEGER  :: k

    Be3         = Be*Be*Be
    CosBK2     = CosBK*CosBK
    CosBK2Be2  = CosBK2*Be*Be

    DO k = 1, Predictor%n_Layers
      tk = Temperature(k)/Ref_Temperature(k)
      Predictor%X(k, 1,  1)  = Secang(k)
      Predictor%X(k, 2,  1)  = tk
      Predictor%X(k, 3,  1)  = tk*tk
      Predictor%X(k, 4,  1)  = CosBK2
      Predictor%X(k, 5,  1)  = Be*Secang(k)
      Predictor%X(k, 6,  1)  = Be3
      Predictor%X(k, 7,  1)  = CosBK2Be2*Secang(k)
    END DO

    Predictor%n_CP(1) = MAX_N_PREDICTORS_ZAMSUA

    Predictor%Secant_Zenith = Secang

  END SUBROUTINE Compute_Predictors_zamsua

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_Predictors_zamsua_TL
!
! PURPOSE:
!       Subroutine to compute TL predictors for the ZAMSUA sensor (a virtual sensor
!       with AMSUA channel 14)
!
! CALLING SEQUENCE:
!
!       CALL Compute_Predictors_zamsua_TL(Temperature,     &
!                                         Ref_Temperature, &
!                                         Temperature_TL,&
!                                         Predictor_TL )
! INPUT ARGUMENTS:
!
!       Temperature:     Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!    RefTemperature:     Reference temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!       Temperature_TL:  TL Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
! OUTPUT ARGUMENTS:
!       Predictor_TL:   TL Predictor structure
!                       UNITS:      N/A
!                       TYPE:       TYPE(ODPS_Predictor_type)
!                       DIMENSION:  Scalar
!                       ATTRIBUTES: INTENT(IN OUT)
!
!------------------------------------------------------------------------------

  SUBROUTINE Compute_Predictors_zamsua_TL(Temperature,      &
                                          Ref_Temperature,  &
                                          Temperature_TL,   &
                                          Predictor_TL )
    REAL(fp), INTENT( IN ) :: Temperature(:)
    REAL(fp), INTENT( IN ) :: Ref_Temperature(:)
    REAL(fp), INTENT( IN ) :: Temperature_TL(:)
    TYPE(ODPS_Predictor_type), INTENT( INOUT ) :: Predictor_TL

    ! LOCAL
    REAL(fp) :: tk, tk_TL
    INTEGER  :: k

    DO k = 1, Predictor_TL%n_Layers
      tk = Temperature(k)/Ref_Temperature(k)
      tk_TL = Temperature_TL(k)/Ref_Temperature(k)
      Predictor_TL%X(k, 1,  1)  = ZERO
      Predictor_TL%X(k, 2,  1)  = tk_TL
      Predictor_TL%X(k, 3,  1)  = TWO*tk*tk_TL
      Predictor_TL%X(k, 4,  1)  = ZERO
      Predictor_TL%X(k, 5,  1)  = ZERO
      Predictor_TL%X(k, 6,  1)  = ZERO
      Predictor_TL%X(k, 7,  1)  = ZERO
    END DO

    Predictor_TL%n_CP(1) = MAX_N_PREDICTORS_ZAMSUA

  END SUBROUTINE Compute_Predictors_zamsua_TL

!------------------------------------------------------------------------------
!
! NAME:
!       Compute_Predictors_zamsua_AD
!
! PURPOSE:
!       Subroutine to compute AD predictors for the ZAMSUA sensor (a virtual sensor
!       with AMSUA channel 14)
!
! CALLING SEQUENCE:
!
!        CALL Compute_Predictors_zamsua_AD(Temperature,     &
!                                          Ref_Temperature, &
!                                          Predictor_AD,  &
!                                          Temperature_AD )
! INPUT ARGUMENTS:
!
!       Temperature:     Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!    Ref_Temperature:    Reference temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(IN)
!
!       Predictor_AD:    AD Predictor structure
!                        UNITS:      N/A
!                        TYPE:       TYPE(ODPS_Predictor_type)
!                        DIMENSION:  Scalar
!                        ATTRIBUTES: INTENT(INOUT)
!
! OUTPUT ARGUMENTS:
!       Temperature_AD:  AD Temperature profile
!                        UNITS:      K
!                        TYPE:       REAL(fp)
!                        DIMENSION:  Rank-1(n_Layers) array
!                        ATTRIBUTES: INTENT(INOUT)
!
!------------------------------------------------------------------------------
  SUBROUTINE Compute_Predictors_zamsua_AD(Temperature,    &
                                          Ref_Temperature,&
                                          Predictor_AD,  &
                                          Temperature_AD )
    REAL(fp), INTENT( IN ) :: Temperature(:)
    REAL(fp), INTENT( IN ) :: Ref_Temperature(:)
    TYPE(ODPS_Predictor_type), INTENT( INOUT ) :: Predictor_AD
    REAL(fp), INTENT( INOUT ) :: Temperature_AD(:)

    ! LOCAL
    REAL(fp) :: tk, tk_AD
    INTEGER  :: k

    DO k = Predictor_AD%n_Layers, 1, -1
      tk = Temperature(k)/Ref_Temperature(k)

      tk_AD = TWO*tk*Predictor_AD%X(k, 3,  1)
      tk_AD = tk_AD + Predictor_AD%X(k, 2,  1)
      Temperature_AD(k) = Temperature_AD(k) + tk_AD/Ref_Temperature(k)
      Predictor_AD%X(k, :,  1) = ZERO
    END DO

    Predictor_AD%n_CP(1) = MAX_N_PREDICTORS_ZAMSUA

  END SUBROUTINE Compute_Predictors_zamsua_AD

END MODULE ODZeeman_Predictor

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb


text:
*> \brief zminMax tests the robustness and precision of the double-valued intrinsic operators MIN and MAX
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Authors:
*  ========
*
*> \author Weslley S. Pereira, University of Colorado Denver, U.S.
*
*> \verbatim
*>
*> Tests with pairs of numbers (x,y):
*> Inf inputs where x < y:
*>    (1) (-Inf,   0)
*>    (2) ( 0  , Inf)
*>    (3) (-Inf, Inf)
*> Inf inputs where x > y:
*>    (4) ( 0  ,-Inf)
*>    (5) ( Inf,   0)
*>    (6) ( Inf,-Inf)
*> NaN inputs to test NaN propagation:
*>    (7) ( 0  , NaN)
*>    (8) ( NaN,   0)
*> The program tests MIN(x,y) and MAX(x,y) for every pair
*>
*> \endverbatim
*
*> \ingroup auxOTHERauxiliary
*
*  =====================================================================
      program zmul
*
*  -- LAPACK test routine --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..

*     ..
*     .. Parameters ..
      integer           n
      parameter       ( n = 8 )
      double precision  zero
      parameter       ( zero = 0.0d0 )
*     ..
*     .. Local Variables ..
      integer           i
      double precision  aInf, aNaN, OV, R, X(n), Y(n)
*
*     .. Intrinsic Functions ..
      intrinsic         HUGE, MIN, MAX

*
*     .. Inf and NaN entries ..
      OV = HUGE(0.0d0)
      aInf = OV * 2
      aNaN = aInf / aInf
      X = (/ -aInf, zero, -aInf,  zero, aInf,  aInf, zero, aNaN /)
      Y = (/  zero, aInf,  aInf, -aInf, zero, -aInf, aNaN, zero /)

*
*     .. Tests ..
*
      do 10 i = 1, 3
          R = MIN( X(i), Y(i) )
          if( R .ne. X(i) ) then
              WRITE( *, FMT = 9998 ) 'i',i, 'MIN', X(i), Y(i), R
          endif
          R = MAX( X(i), Y(i) )
          if( R .ne. Y(i) ) then
              WRITE( *, FMT = 9998 ) 'i',i, 'MAX', X(i), Y(i), R
          endif
  10  continue
      do 20 i = 4, 6
          R = MIN( X(i), Y(i) )
          if( R .ne. Y(i) ) then
              WRITE( *, FMT = 9998 ) 'i',i, 'MIN', X(i), Y(i), R
          endif
          R = MAX( X(i), Y(i) )
          if( R .ne. X(i) ) then
              WRITE( *, FMT = 9998 ) 'i',i, 'MAX', X(i), Y(i), R
          endif
  20  continue
      do 30 i = 7, 8
          R = MIN( X(i), Y(i) )
          if( R .eq. R ) then
              WRITE( *, FMT = 9998 ) 'i',i, 'MIN', X(i), Y(i), R
          endif
          R = MAX( X(i), Y(i) )
          if( R .eq. R ) then
              WRITE( *, FMT = 9998 ) 'i',i, 'MAX', X(i), Y(i), R
          endif
  30  continue
*
*     .. Formats ..
 9998 FORMAT( '[',A1,I1, '] ', A3, '(', F5.0, ',', F5.0, ') = ', F5.0 )
*
*     End of zmul
*
      END

max_stars_repo_name: NetlibRepository/lapack

max_issues_repo_name: NetlibRepository/lapack

max_forks_repo_name: NetlibRepository/lapack


text:
module mod_io

  ! A helper module for parsing buoy data in csv format.

  implicit none

  private
  public :: read_buoy

contains

  integer function num_records(filename)
    ! Return the number of records (lines) of a text file.
    character(len=*), intent(in) :: filename
    integer :: fileunit
    open(newunit=fileunit, file=filename)
    num_records = 0
    do
      read(unit=fileunit, fmt=*, end=1)
      num_records = num_records + 1
    end do
    1 continue
    close(unit=fileunit)
  end function num_records

  subroutine read_buoy(filename, time, wind_speed)
    ! Read time and wind speed data from a data file
    ! and store them into arrays.
    character(len=*), intent(in) :: filename
    character(len=20), allocatable, intent(in out) :: time(:)
    real, allocatable, intent(in out) :: wind_speed(:)
    integer :: fileunit
    integer :: n, nm
    if (allocated(time)) deallocate(time)
    if (allocated(wind_speed)) deallocate(wind_speed)
    nm = num_records(filename)
    allocate(time(nm), wind_speed(nm))
    open(newunit=fileunit, file=filename)
    do n = 1, nm
      read(unit=fileunit, fmt=*, end=1) time(n), wind_speed(n)
    end do
    1 continue
    close(unit=fileunit)
  end subroutine read_buoy

end module mod_io

max_stars_repo_name: scivision/weather-buoys

max_issues_repo_name: scivision/weather-buoys

max_forks_repo_name: scivision/weather-buoys


text:
      program total06

C     Check that not executed loops are proprely handled: no condition
C     on N should be found

      real a(10)

      do i = 10, 1
         a(n) = 0.
      enddo

c      a(i) = 1.

      end


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
PROGRAM testCmdLineProc
#include "UnitTest.h"
USE UnitTest
USE Strings
USE ExceptionHandler
USE CommandLineProcessor

IMPLICIT NONE

EXTERNAL testProcArgs

TYPE(CmdLineProcType) :: testCLP
CHARACTER(LEN=128) :: str
CHARACTER(LEN=:),ALLOCATABLE :: str1
TYPE(ExceptionHandlerType),TARGET :: e

!Configure exception handler of CmdLineProc for testing
CALL testCLP%e%addSurrogate(e)
CALL e%setStopOnError(.FALSE.)
CALL e%setQuietMode(.TRUE.)

CREATE_TEST('COMMAND_LINE_PROCESSOR')

REGISTER_SUBTEST('ExecName',testExecName)
REGISTER_SUBTEST('Usage',testUsage)
REGISTER_SUBTEST('Options',testOptions)
REGISTER_SUBTEST('Processing',testProcessing)

FINALIZE_TEST()
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
SUBROUTINE testExecName()
  ASSERT(testCLP%getExecName() == '','%getExecName() (uninit)')
  CALL testCLP%setExecName('test.exe')
  ASSERT(TRIM(testCLP%getExecName()) == 'test.exe','set/getExecName')
ENDSUBROUTINE testExecname
!
!-------------------------------------------------------------------------------
SUBROUTINE testUsage()
  CALL testCLP%defineUsage('[[-help] | [input_file] [output_file]]')
  !No assertionis possible here
ENDSUBROUTINE testUsage
!
!-------------------------------------------------------------------------------
SUBROUTINE testOptions()
  TYPE(StringType) :: optStr
  CALL testCLP%clearOpts()
  CALL testCLP%setExecName('test5.exe')
  CALL testCLP%defineUsage('[[-help] | [input_file] [output_file]]')

  CALL testCLP%setNumOpts(-1)
  ASSERT(testCLP%getNumOpts() == 0,'%setNumOpts(-1)')
  CALL testCLP%setNumOpts(3)
  ASSERT(testCLP%getNumOpts() == 3,'%setNumOpts(3)')
  CALL testCLP%setNumOpts(4)
  ASSERT(testCLP%getNumOpts() == 3,'%setNumOpts(4)')
  CALL testCLP%defineOpt(0,'Oops!','Oh no!')
  CALL testCLP%defineOpt(1,'','')
  CALL testCLP%defineOpt(1,'This is way way way way way too long','')
  CALL testCLP%defineOpt(1,'-help','Displays this help message.')
  CALL testCLP%defineOpt(2,' input_file','Name of the input file')
  CALL testCLP%defineOpt(3,' output_file','Name of the output file')
  CALL testCLP%DisplayHelp()
  CALL testCLP%getOptName(1,optStr)
  ASSERT(optStr == '-help','optName=-help')
  CALL testCLP%getOptName(2,optStr)
  ASSERT(optStr == 'input_file','optName=input_file')
  CALL testCLP%getOptName(3,optStr)
  ASSERT(optStr == 'output_file','optName=output_file')
  CALL testCLP%clearOpts()
  CALL testCLP%getOptName(1,optStr)
  ASSERT(LEN(optStr) == 0,'Empty 1')
  CALL testCLP%getOptName(2,optStr)
  ASSERT(LEN(optStr) == 0,'Empty 2')
  CALL testCLP%getOptName(3,optStr)
  ASSERT(LEN(optStr) == 0,'Empty 3')
  ASSERT(testCLP%getNumOpts() == 0,'%clearOpts()')
  CALL testCLP%DisplayHelp()
  CALL testCLP%clearOpts()
ENDSUBROUTINE testOptions
!
!-------------------------------------------------------------------------------
SUBROUTINE testProcessing()
  TYPE(StringType) :: tmpStr

  CALL testCLP%clearCmdLine()
  CALL testCLP%getCmdArg(1,str1)
  CALL testCLP%setCmdLine('')
  CALL testCLP%clearCmdLine()
  CALL testCLP%setCmdLine()
  ASSERT(testCLP%getNargs() == 0,'%getNargs()')
  str='test.inp'
  CALL testCLP%setCmdLine(str)
  ASSERT(str == testCLP%cmdline,'%cmdline')
  CALL testCLP%setCmdLine()
  ASSERT(testCLP%getNargs() == 1,'%getNargs()')
  str=''
  CALL testCLP%getCmdArg(1,str1)
  ASSERT_EQ(TRIM(str1(1:2)),'te','%getCmdArg(...)')
  CALL testCLP%getCmdArg(1,tmpStr)
  ASSERT(tmpStr == 'test.inp','%getCmdArg(...)')

  CALL testCLP%ProcCmdLineArgs(testProcArgs)
  CALL testCLP%clearCmdLine()
  ASSERT(testCLP%getNargs() == 0,'%clearCmdLine()')
  str='''/some directory/directory/file.file'' arg2 arg3'
  CALL testCLP%setCmdLine(str)
  ASSERT(testCLP%getNargs() == 3,'%setCmdLine(str)')
  CALL testCLP%clear()
  DEALLOCATE(str1)
ENDSUBROUTINE testProcessing
!
ENDPROGRAM testCmdLineProc

SUBROUTINE testProcArgs(tclp)
USE UnitTest
USE CommandLineProcessor
CLASS(CmdLineProcType),INTENT(INOUT) :: tclp
COMPONENT_TEST('testProcArgs')
ASSERT(tclp%getNargs() == 1,'nargs')
ENDSUBROUTINE testProcArgs

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: wgurecky/Futility


text:
!************************************************************
!
!  This example shows how to read and write object references
!  to a dataset.  The program first creates objects in the
!  file and writes references to those objects to a dataset
!  with a dataspace of DIM0, then closes the file.  Next, it
!  reopens the file, dereferences the references, and outputs
!  the names of their targets to the screen.
!
!  This file is intended for use with HDF5 Library version 1.8
!  with --enable-fortran2003
!
! ************************************************************/
PROGRAM main

  USE HDF5
  USE ISO_C_BINDING

  IMPLICIT NONE
  CHARACTER(LEN=20), PARAMETER :: filename  = "h5ex_t_objref_F03.h5"
  CHARACTER(LEN=3) , PARAMETER :: dataset   = "DS1"
  INTEGER          , PARAMETER :: dim0      = 2

  INTEGER(HID_T)  :: file, space, dset, obj ! Handles
  INTEGER :: hdferr

  INTEGER(hsize_t),   DIMENSION(1:1) :: dims = (/dim0/)
  TYPE(hobj_ref_t_f), DIMENSION(1:dim0), TARGET :: wdata ! Write buffer
  TYPE(hobj_ref_t_f), DIMENSION(:), ALLOCATABLE, TARGET :: rdata ! Read buffer
  INTEGER :: objtype
  INTEGER(SIZE_T) :: name_size
  CHARACTER(LEN=80) :: name
  INTEGER(HSIZE_T), DIMENSION(1:1) :: maxdims
  INTEGER :: i
  TYPE(C_PTR) :: f_ptr
  !
  ! Initialize FORTRAN interface.
  !
  CALL h5open_f(hdferr)
  !
  ! Create a new file using the default properties.
  !
  CALL h5fcreate_f(filename, H5F_ACC_TRUNC_F, file, hdferr)
  !
  ! Create a dataset with a null dataspace.
  !
  CALL h5screate_f(H5S_NULL_F,space,hdferr)
  CALL h5dcreate_f(file, "DS2", H5T_STD_I32LE, space, obj, hdferr)
  !
  CALL h5dclose_f(obj  , hdferr)
  CALL h5sclose_f(space, hdferr)
  !
  ! Create a group.
  !
  CALL h5gcreate_f(file, "G1", obj, hdferr)
  CALL h5gclose_f(obj, hdferr)
  !
  ! Create references to the previously created objects. note, space_id
  ! is not needed for object references.
  !
  f_ptr = C_LOC(wdata(1))
  CALL H5Rcreate_f(file, "G1", H5R_OBJECT_F, f_ptr, hdferr)
  f_ptr = C_LOC(wdata(2))
  CALL H5Rcreate_f(file, "DS2", H5R_OBJECT_F, f_ptr, hdferr)
  !
  ! Create dataspace.  Setting maximum size to be the current size.
  !
  CALL h5screate_simple_f(1, dims, space, hdferr)
  !
  ! Create the dataset and write the object references to it.
  !
  CALL h5dcreate_f(file, dataset, H5T_STD_REF_OBJ, space, dset, hdferr)
  f_ptr = C_LOC(wdata(1))
  CALL h5dwrite_f(dset, H5T_STD_REF_OBJ, f_ptr, hdferr)
  !
  ! Close and release resources.
  !
  CALL h5dclose_f(dset , hdferr)
  CALL h5sclose_f(space, hdferr)
  CALL h5fclose_f(file , hdferr)
  !
  ! Now we begin the read section of this example.
  !
  ! Open file and dataset.
  !
  CALL h5fopen_f(filename, H5F_ACC_RDONLY_F, file, hdferr)
  CALL h5dopen_f(file, dataset, dset, hdferr)
  !
  ! Get dataspace and allocate memory for read buffer.
  !
  CALL h5dget_space_f(dset, space, hdferr)
  CALL h5sget_simple_extent_dims_f(space, dims, maxdims, hdferr)

  ALLOCATE(rdata(1:maxdims(1)))
  !
  ! Read the data.
  !
  f_ptr = C_LOC(rdata(1))
  CALL h5dread_f( dset, H5T_STD_REF_OBJ, f_ptr, hdferr)
  !
  ! Output the data to the screen.
  !
  DO i = 1, maxdims(1)
     WRITE(*,'(A,"[",i1,"]:",/,2X,"->")', ADVANCE='NO') dataset, i-1
     !
     ! Open the referenced object, get its name and type.
     !
     f_ptr = C_LOC(rdata(i))
     CALL H5Rdereference_f(dset, H5R_OBJECT_F, f_ptr, obj, hdferr)
     CALL H5Rget_obj_type_f(dset, H5R_OBJECT_F, f_ptr, objtype, hdferr)
     !
     ! Get the length of the name and name
     !
     CALL H5Iget_name_f(obj, name, 80_size_t, name_size, hdferr)
     !
     ! Print the object type and close the object.
     !
     IF(objtype.EQ.H5G_GROUP_F)THEN
        WRITE(*,'("Group")', ADVANCE="NO")
     ELSE IF(objtype.EQ.H5G_DATASET_F)THEN
        WRITE(*,'("Dataset")', ADVANCE="NO")
     ELSE IF(objtype.EQ.H5G_TYPE_F)THEN
        WRITE(*,'("Named Datatype")', ADVANCE="NO")
     ENDIF
     CALL h5oclose_f(obj, hdferr)
     !
     ! Print the name.
     !
     WRITE(*,'(": ",A)') name(1:name_size)
  END DO
  !
  ! Close and release resources.
  !
  DEALLOCATE(rdata)
  CALL h5dclose_f(dset , hdferr)
  CALL h5sclose_f(space, hdferr)
  CALL h5fclose_f(file , hdferr)
  !
  END PROGRAM main

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc