subroutine enorm_state(xst,enorm,yst)
!$$$ subprogram documentation block
! . . . .
! subprogram: enorm_state
! prgmmr: tremolet
!
! abstract: compute energy norm on GSI grid
!
! program history log:
! 2007-10-19 tremolet - initial code
! 2009-01-18 todling - carry summation in quad precision
! 2009-08-14 lueken - update documentation
! 2010-05-13 todling - update to use gsi_bundle
! 2010-09-06 todling - revisit pointers; add Q-term to norm
! 2011-08-01 lueken - replaced F90 with f90 (no machine logic)
! 2013-10-19 todling - metguess now holds background
!
! input argument list:
! xst
! yst
!
! output argument list:
! enorm
! yst
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind,r_quad
use constants, only: zero,one,cp,rd,pi,two,r1000,hvap
use jcmod, only: eps_eer
use gridmod, only: nsig,nlat,nlon,lon2,lat2,istart,rlats,ak5,bk5
use mpimod, only: mype
use guess_grids, only: ntguessig
use state_vectors, only: dot_product
use gsi_bundlemod, only: assignment(=)
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_metguess_mod, only: gsi_metguess_bundle
implicit none
type(gsi_bundle), intent(inout) :: xst
real(r_quad) , intent( out) :: enorm
type(gsi_bundle), intent(inout) :: yst
! Declare local variables
real(r_kind) :: qfact, tfact, pfact, pref, tref, gridfac, zps
real(r_kind) :: coslat(lat2), dsig(lat2,lon2,nsig), akk(nsig)
integer(i_kind) :: ii,jj,kk,ilat,ier
real(r_kind), parameter :: Pa_per_kPa = r1000
real(r_kind),pointer,dimension(:,:,:) :: xst_u,xst_v,xst_t,xst_q
real(r_kind),pointer,dimension(:,:,:) :: yst_u,yst_v,yst_t,yst_q
real(r_kind),pointer,dimension(:,:) :: xst_p
real(r_kind),pointer,dimension(:,:) :: yst_p
real(r_kind),pointer,dimension(:,:) :: ges_ps
! ----------------------------------------------------------------------
tref=280.0_r_kind
pref=1.0e5_r_kind
tfact=cp/tref
pfact=rd*tref/(pref*pref)
qfact=abs(eps_eer)*hvap*hvap/(cp*tref)
gridfac=one/(nlat*nlon)
do kk=1,nsig
akk(kk) = Pa_per_kPa * (ak5(kk)-ak5(kk+1))
enddo
call gsi_bundlegetpointer(gsi_metguess_bundle(ntguessig),'ps', ges_ps,ier)
dsig=HUGE(dsig)
do jj=2,lon2-1
do ii=2,lat2-1
zps = Pa_per_kPa * ges_ps(ii,jj)
do kk=1,nsig
dsig(ii,jj,kk) = akk(kk) + (bk5(kk)-bk5(kk+1)) * zps
enddo
if (ANY(dsig(ii,jj,:)<=zero)) then
do kk=1,nsig
write(6,'(A,I3,4(2X,F18.8))')'enorm ak,bk,pk,dsig=',&
kk,r1000*ak5(kk),bk5(kk),r1000*ak5(kk)+bk5(kk)*pref,dsig(2,2,kk)
enddo
write(6,'(A,I3,4(2X,F18.8))')'enorm ak,bk,pk =',&
nsig+1,r1000*ak5(nsig+1),bk5(nsig+1),r1000*ak5(nsig+1)+bk5(nsig+1)*pref
write(6,*)'enorm_state: negative dsig'
call stop2(123)
endif
enddo
enddo
coslat=HUGE(coslat)
do ii=2,lat2-1
ilat=istart(mype+1)+ii-2
if (ilat<1.or.ilat>nlat) then
write(6,*)'enorm_state: error ilat',ilat
call stop2(124)
end if
coslat(ii)=cos(rlats(ilat))
if (coslat(ii)<zero) then
write(6,*)'enorm_state error coslat:',ii,ilat,rlats(ilat),pi/two,coslat(ii)
coslat(ii)=-coslat(ii)
endif
enddo
! ----------------------------------------------------------------------
yst=zero
! ----------------------------------------------------------------------
call gsi_bundlegetpointer(xst,'u', xst_u,ier)
call gsi_bundlegetpointer(xst,'v', xst_v,ier)
call gsi_bundlegetpointer(xst,'tv',xst_t,ier)
call gsi_bundlegetpointer(xst,'q', xst_q,ier)
call gsi_bundlegetpointer(xst,'ps',xst_p,ier)
call gsi_bundlegetpointer(yst,'u', yst_u,ier)
call gsi_bundlegetpointer(yst,'v', yst_v,ier)
call gsi_bundlegetpointer(yst,'tv',yst_t,ier)
call gsi_bundlegetpointer(yst,'q', yst_q,ier)
call gsi_bundlegetpointer(yst,'ps',yst_p,ier)
! U
do kk=1,nsig
do jj=2,lon2-1
do ii=2,lat2-1
yst_u(ii,jj,kk)=gridfac*coslat(ii)*dsig(ii,jj,kk)*xst_u(ii,jj,kk)
enddo
enddo
enddo
! V
do kk=1,nsig
do jj=2,lon2-1
do ii=2,lat2-1
yst_v(ii,jj,kk)=gridfac*coslat(ii)*dsig(ii,jj,kk)*xst_v(ii,jj,kk)
enddo
enddo
enddo
! T
do kk=1,nsig
do jj=2,lon2-1
do ii=2,lat2-1
yst_t(ii,jj,kk)=gridfac*coslat(ii)*dsig(ii,jj,kk)*tfact*xst_t(ii,jj,kk)
enddo
enddo
enddo
! Q
if(eps_eer>zero) then
do kk=1,nsig
do jj=2,lon2-1
do ii=2,lat2-1
yst_q(ii,jj,kk)=gridfac*coslat(ii)*dsig(ii,jj,kk)*qfact*xst_q(ii,jj,kk)
enddo
enddo
enddo
endif
! P
do jj=2,lon2-1
do ii=2,lat2-1
yst_p(ii,jj)=gridfac*coslat(ii)*pfact*xst_p(ii,jj)
enddo
enddo
! ----------------------------------------------------------------------
enorm=DOT_PRODUCT(yst,xst)
! ----------------------------------------------------------------------
return
end subroutine enorm_state
subroutine enorm_state_red(xst,enorm,yst)
!$$$ subprogram documentation block
! . . . .
! subprogram: enorm_state
! prgmmr: tremolet
!
! abstract: compute energy norm on GSI grid
!
! program history log:
! 2007-10-19 tremolet - initial code
! 2009-01-18 todling - carry summation in quad precision
! 2009-08-14 lueken - update documentation
! 2010-05-13 todling - update to use gsi_bundle
! 2010-09-06 todling - revisit pointers; add Q-term to norm
! 2011-08-01 lueken - replaced F90 with f90 (no machine logic)
! 2013-10-19 todling - metguess now holds background
!
! input argument list:
! xst
! yst
!
! output argument list:
! enorm
! yst
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind,r_quad
use constants, only: zero,one,cp,rd,pi,two,r1000,hvap
use jcmod, only: eps_eer
use gridmod, only: nsig,nlat,nlon,lon2,lat2,istart,rlats,ak5,bk5
use mpimod, only: mype
use guess_grids, only: ntguessig
use state_vectors, only: dot_product_red
use gsi_bundlemod, only: assignment(=)
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_metguess_mod, only: gsi_metguess_bundle
implicit none
type(gsi_bundle), intent(inout) :: xst
real(r_quad) , intent( out) :: enorm
type(gsi_bundle), intent(inout) :: yst
! Declare local variables
real(r_kind) :: qfact, tfact, pfact, pref, tref, gridfac, zps
real(r_kind) :: coslat(lat2), dsig(lat2,lon2,nsig), akk(nsig)
integer(i_kind) :: ii,jj,kk,ilat,ier
real(r_kind), parameter :: Pa_per_kPa = r1000
real(r_kind),pointer,dimension(:,:,:) :: xst_u,xst_v,xst_t,xst_q
real(r_kind),pointer,dimension(:,:,:) :: yst_u,yst_v,yst_t,yst_q
real(r_kind),pointer,dimension(:,:) :: xst_p
real(r_kind),pointer,dimension(:,:) :: yst_p
real(r_kind),pointer,dimension(:,:) :: ges_ps
! ----------------------------------------------------------------------
tref=280.0_r_kind
pref=1.0e5_r_kind
tfact=cp/tref
pfact=rd*tref/(pref*pref)
qfact=abs(eps_eer)*hvap*hvap/(cp*tref)
gridfac=one/(nlat*nlon)
do kk=1,nsig
akk(kk) = Pa_per_kPa * (ak5(kk)-ak5(kk+1))
enddo
call gsi_bundlegetpointer(gsi_metguess_bundle(ntguessig),'ps', ges_ps,ier)
dsig=HUGE(dsig)
do jj=2,lon2-1
do ii=2,lat2-1
zps = Pa_per_kPa * ges_ps(ii,jj)
do kk=1,nsig
dsig(ii,jj,kk) = akk(kk) + (bk5(kk)-bk5(kk+1)) * zps
enddo
if (ANY(dsig(ii,jj,:)<=zero)) then
do kk=1,nsig
write(6,'(A,I3,4(2X,F18.8))')'enorm ak,bk,pk,dsig=',&
kk,r1000*ak5(kk),bk5(kk),r1000*ak5(kk)+bk5(kk)*pref,dsig(2,2,kk)
enddo
write(6,'(A,I3,4(2X,F18.8))')'enorm ak,bk,pk =',&
nsig+1,r1000*ak5(nsig+1),bk5(nsig+1),r1000*ak5(nsig+1)+bk5(nsig+1)*pref
write(6,*)'enorm_state: negative dsig'
call stop2(123)
endif
enddo
enddo
coslat=HUGE(coslat)
do ii=2,lat2-1
ilat=istart(mype+1)+ii-2
if (ilat<1.or.ilat>nlat) then
write(6,*)'enorm_state: error ilat',ilat
call stop2(124)
end if
coslat(ii)=cos(rlats(ilat))
if (coslat(ii)<zero) then
write(6,*)'enorm_state error coslat:',ii,ilat,rlats(ilat),pi/two,coslat(ii)
coslat(ii)=-coslat(ii)
endif
enddo
! ----------------------------------------------------------------------
yst=zero
! ----------------------------------------------------------------------
call gsi_bundlegetpointer(xst,'u', xst_u,ier)
call gsi_bundlegetpointer(xst,'v', xst_v,ier)
call gsi_bundlegetpointer(xst,'tv',xst_t,ier)
call gsi_bundlegetpointer(xst,'q', xst_q,ier)
call gsi_bundlegetpointer(xst,'ps',xst_p,ier)
call gsi_bundlegetpointer(yst,'u', yst_u,ier)
call gsi_bundlegetpointer(yst,'v', yst_v,ier)
call gsi_bundlegetpointer(yst,'tv',yst_t,ier)
call gsi_bundlegetpointer(yst,'q', yst_q,ier)
call gsi_bundlegetpointer(yst,'ps',yst_p,ier)
! U
do kk=1,nsig
do jj=2,lon2-1
do ii=2,lat2-1
yst_u(ii,jj,kk)=gridfac*coslat(ii)*dsig(ii,jj,kk)*xst_u(ii,jj,kk)
enddo
enddo
enddo
! V
do kk=1,nsig
do jj=2,lon2-1
do ii=2,lat2-1
yst_v(ii,jj,kk)=gridfac*coslat(ii)*dsig(ii,jj,kk)*xst_v(ii,jj,kk)
enddo
enddo
enddo
! T
do kk=1,nsig
do jj=2,lon2-1
do ii=2,lat2-1
yst_t(ii,jj,kk)=gridfac*coslat(ii)*dsig(ii,jj,kk)*tfact*xst_t(ii,jj,kk)
enddo
enddo
enddo
! Q
if(eps_eer>zero) then
do kk=1,nsig
do jj=2,lon2-1
do ii=2,lat2-1
yst_q(ii,jj,kk)=gridfac*coslat(ii)*dsig(ii,jj,kk)*qfact*xst_q(ii,jj,kk)
enddo
enddo
enddo
endif
! P
do jj=2,lon2-1
do ii=2,lat2-1
yst_p(ii,jj)=gridfac*coslat(ii)*pfact*xst_p(ii,jj)
enddo
enddo
! ----------------------------------------------------------------------
enorm=DOT_PRODUCT_red(yst,xst,0)
! ----------------------------------------------------------------------
return
end subroutine enorm_state_red