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 ! ! 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: ione,zero,one,cp,rd,pi,two,r1000 use gridmod, only: nsig,nlat,nlon,lon2,lat2,istart,rlats,ak5,bk5 use mpimod, only: mype use guess_grids, only: ges_ps,ntguessig use state_vectors implicit none type(state_vector), intent(in ) :: xst real(r_quad) , intent( out) :: enorm type(state_vector), intent(inout) :: yst ! Declare local variables real(r_kind) :: tfact, pfact, pref, tref, gridfac, zps real(r_kind) :: coslat(lat2), dsig(lat2,lon2,nsig), akk(nsig) integer(i_kind) :: ii,jj,kk,ijk,ilat real(r_kind), parameter :: Pa_per_kPa = r1000 ! ---------------------------------------------------------------------- tref=280.0_r_kind pref=1.0e5_r_kind tfact=cp/tref pfact=rd*tref/(pref*pref) gridfac=one/(nlat*nlon) do kk=1,nsig akk(kk) = Pa_per_kPa * (ak5(kk)-ak5(kk+ione)) enddo dsig=HUGE(dsig) do jj=2,lon2-ione do ii=2,lat2-ione zps = Pa_per_kPa * ges_ps(ii,jj,ntguessig) do kk=1,nsig dsig(ii,jj,kk) = akk(kk) + (bk5(kk)-bk5(kk+ione)) * 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+ione,r1000*ak5(nsig+ione),bk5(nsig+ione),r1000*ak5(nsig+ione)+bk5(nsig+ione)*pref write(6,*)'enorm_state: negative dsig' call stop2(123) endif enddo enddo coslat=HUGE(coslat) do ii=2,lat2-ione ilat=istart(mype+ione)+ii-2_i_kind if (ilatnlat) then write(6,*)'enorm_state: error ilat',ilat call stop2(124) end if coslat(ii)=cos(rlats(ilat)) if (coslat(ii)