subroutine compute_derived(mype,init_pass)
!$$$ subprogram documentation block
! . . . .
! subprogram: compute_derived compute derived quantites from current solution
! prgmmr: derber org: np2 date: 2005-11-29
!
! abstract: This routine performs various functions, all related in one
! way or the other to the model guess or current solution.
!
! Functions performed in this routine include the following
! a) compute guess-derived fields required by pcp forward model
! b) compute 3d pressure grids
! c) compute saturation specific humidity. on first outer iteration
! save qs for normalization in background error. qs for limq
! is updated each outer iteration.
! d) compute 3d geopotential height
! e) compute 2d tropopause pressure map
!
! program history log:
! 2005-11-21 derber - new routine from read_*.f90 routines
! 2005-11-22 wu - bug fix qoption=2 for regional nmm on dqdp and
! set dqdt=0 above sigma=0.15; for regional mass
! set dqdt=dqdp=0 above sigma=0.15
! 2005-12-09 guo - remove GMAO derivative computation code. Use
! unified NCEP compact_diff procedures.
! 2006-01-09 derber - include calculation of sigsum add capability of set_nrh_var
! 2006-01-30 kleist - correct tropprs unit error in qoption=2 q/t decoupling
! 2006-02-02 treadon - consolidate/unify use of guess pressure arrays
! 2006-02-03 derber - clean up RH statistics printout
! 2006-03-07 treadon - remove ges_prslk and related code
! 2006-03-27 treadon - remove guess bias correction arrays since not used
! 2006-04-17 treadon - replace sigi with bk5; replace sigl with
! ges_prslavg/ges_psfcavg
! 2006-04-21 kleist - modify call to calctends
! 2006-07-31 kleist - changes to use ps instead of ln(ps)
! 2006-09-29 treadon - add option to compute 10m wind factor fields
! 2007-03-13 derber - add changes to make qoption=2 variances work as others
! 2007-05-08 kleist - remove jcdivt from use list
! 2007-06-21 rancic - add pbl code
! 2007-07-26 cucurull - call gesprs, add ges_3dp and remove ps
! in calctends argument list
! 2007-08-08 derber - pass ges_teta to calctends rather than calculate seperately
! 2008-06-05 safford - rm unused uses
! 2008-10-10 derber - add calculation of fact_tv
! 2008-11-03 sato - add anisotropic mode procedures
! 2008-12-08 todling - move 3dprs/geop-hght calculation from here into setuprhsall
! 2009-08-19 guo - add verifications of drv_initialized and tnd_initialized
! before the use of related module variables.
! 2009-10-15 parrish - add rescale of ensemble rh perturbations
! (currently for internal generated ensemble only)
! 2010-03-11 derber/zhu - add qvar3d to prewgt and prewgt_reg, remove rescale_ensemble_rh_perturbations
! 2010-04-15 hou - add importing gfs co2 historical data through calling
! read_gfsco2.
!
! input argument list:
! mype - mpi task id
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use jfunc, only: qsatg,qgues,jiter,jiterstart,&
dqdt,dqdp,qoption,switch_on_derivatives,&
tendsflag,varq
use control_vectors, only: nrf3_q,nrf3_loc
use mpimod, only: levs_id
use guess_grids, only: ges_z,ges_ps,ges_u,ges_v,&
ges_tv,ges_q,ges_oz,ges_cwmr,ges_tsen,sfct,&
tropprs,ges_prsl,ntguessig,&
nfldsig,&
ges_teta,fact_tv, &
ges_u_lon,ges_v_lon,ges_tvlon,ges_ps_lon,ges_qlon,ges_ozlon,ges_cwmr_lon, &
ges_u_lat,ges_v_lat,ges_tvlat,ges_ps_lat,ges_qlat,ges_ozlat,ges_cwmr_lat
use guess_grids, only: ges_u_ten,ges_v_ten,ges_tv_ten,ges_prs_ten,ges_q_ten,&
ges_oz_ten,ges_cwmr_ten
use guess_grids, only: ges_prslavg,ges_psfcavg
use guess_grids, only: tnd_initialized
use guess_grids, only: drv_initialized
use guess_grids, only: ges_prsi,ges_co2,igfsco2
use gridmod, only: lat2,lon2,nsig,nnnn1o,aeta2_ll,nsig1o
use gridmod, only: regional
use gridmod, only: twodvar_regional
use gridmod, only: wrf_nmm_regional,wrf_mass_regional,nems_nmmb_regional
use gridmod, only: rlats,istart,nlat
use berror, only: hswgt
use balmod, only: rllat1,llmax
use mod_strong, only: jcstrong,baldiag_full
use obsmod, only: write_diag,lobserver
use hybrid_ensemble_parameters, only: l_hyb_ens,generate_ens
use hybrid_ensemble_isotropic_regional, only: rescale_ensemble_rh_perturbations
use obsmod, only: iadate
use ncepgfs_ghg, only: read_gfsco2
use constants, only: zero,one,one_tenth,half,fv
! for anisotropic mode
use sub2fslab_mod, only: setup_sub2fslab, sub2fslab, sub2fslab_glb, destroy_sub2fslab
use anberror, only: anisotropic, idvar, kvar_start, an_amp0, ngauss, indices, indices_p, &
& filter_all, filter_p2, filter_p3, &
& pf2aP1, pf2aP2, pf2aP3, rtma_subdomain_option
use anisofilter, only: rh0f, corz, ensamp, mlat, rllatf, fact_qopt2
use anisofilter_glb, only: rh2f, rh3f, ensamp0f, ensamp2f, ensamp3f, &
p0ilatf, p2ilatf, p3ilatf, p2ilatfm, p3ilatfm, get_stat_factk
use constants, only: zero,one,one_tenth,half,fv
use mpeu_util, only: die, tell
implicit none
! Declare passed variables
integer(i_kind),intent(in ) :: mype
logical ,intent(in ) :: init_pass
! Declare local variables
logical ice,fullfield
integer(i_kind) i,j,k,it,k150,kpres,n,np,l,l2,iderivative
integer(i_kind) :: iyear, month
real(r_kind) d,dl1,dl2,psfc015,dn1,dn2
real(r_kind),dimension(lat2,lon2,nsig+1):: ges_3dp
real(r_kind),dimension(lat2,lon2,nfldsig):: sfct_lat,sfct_lon
real(r_kind),dimension(lat2,lon2,nsig):: rhgues
real(r_kind),dimension(lat2):: xlats
real(r_kind),dimension(lat2,lon2,nsig+1):: prsi
! for anisotropic mode
integer(i_kind):: k1,ivar,kvar,igauss
real(r_kind):: factor,factk,hswgtsum
if(init_pass .and. (ntguessig<1 .or. ntguessig>nfldsig)) &
call die('compute_derived','invalid init_pass, ntguessig =',ntguessig)
! Limit q to be >= 1.e-10_r_kind
do it=1,nfldsig
do k=1,nsig
do j=1,lon2
do i=1,lat2
ges_q(i,j,k,it)=max(ges_q(i,j,k,it),1.e-10_r_kind)
end do
end do
end do
end do
!-----------------------------------------------------------------------------------
! Compute derivatives for .not. twodvar_regional case
if (.not. twodvar_regional)then
if (switch_on_derivatives) then
if(.not.drv_initialized) &
call die('compute_derived','unexpected drv_initialized =',drv_initialized)
! Instead, update gradients of all guess fields. these will
! be used for forward models that need gradient of background field,
! and for getting time derivatives of prognostic variables for
! time extrapolation and non-linear balance constraints.
it=ntguessig
call get_derivatives(ges_u(1,1,1,it),ges_v(1,1,1,it), &
ges_tv(1,1,1,it),ges_ps,ges_q(1,1,1,it),&
ges_oz(1,1,1,it),sfct(1,1,it),ges_cwmr(1,1,1,it), &
ges_u_lon,ges_v_lon,ges_tvlon,ges_ps_lon,ges_qlon,&
ges_ozlon,sfct_lon,ges_cwmr_lon, &
ges_u_lat,ges_v_lat,ges_tvlat,ges_ps_lat,ges_qlat,&
ges_ozlat,sfct_lat,ges_cwmr_lat, &
nnnn1o,nfldsig)
if(.not. wrf_mass_regional .and. tendsflag)then
if(.not.tnd_initialized) &
call die('compute_derived','unexpected tnd_initialized =',tnd_initialized)
! now that we have derivs, get time tendencies if necessary
if(init_pass) then
it=ntguessig
call getprs(ges_ps(1,1,it),ges_3dp)
call calctends(ges_u(1,1,1,it),ges_v(1,1,1,it),ges_tv(1,1,1,it), &
ges_q(1,1,1,it),ges_oz(1,1,1,it),ges_cwmr(1,1,1,it),&
ges_teta(1,1,1,it),ges_z(1,1,it), &
ges_u_lon,ges_u_lat,ges_v_lon,&
ges_v_lat,ges_tvlon,ges_tvlat,ges_ps_lon(1,1,it), &
ges_ps_lat(1,1,it),ges_qlon,ges_qlat,ges_ozlon,&
ges_ozlat,ges_cwmr_lon,ges_cwmr_lat,&
mype,ges_u_ten,ges_v_ten,ges_tv_ten,ges_prs_ten,ges_q_ten,&
ges_oz_ten,ges_cwmr_ten,ges_3dp)
if(jcstrong .and. write_diag(jiter) .and. baldiag_full) then
fullfield=.true.
call strong_bal_correction(ges_u_ten,ges_v_ten,ges_tv_ten,ges_prs_ten,mype, &
ges_u,ges_v,ges_tv,ges_ps,.true.,fullfield,.false.)
end if
end if ! (init_pass)
end if
end if
if(init_pass) then
! Compute tropopause level (in pressure, hPa). The 'pvoz'
! string means compute tropopause using potential vorticity
! and ozone. The 'temp' string means compute tropopause
! using WMO temperature lapse rate method.
! NOTE: tropopause pressure is not needed for 2dvar option
if(regional)then
call tpause(mype,'temp')
else ! (regional)
call tpause(mype,'pvoz')
end if ! (regional)
endif ! (init_pass)
endif ! (!twodvar_regional)
if(.not. init_pass) return
! Load guess q for use in limq. Initialize saturation array to guess.
do k=1,nsig
do j=1,lon2
do i=1,lat2
qgues(i,j,k)=ges_q(i,j,k,ntguessig) ! q guess
fact_tv(i,j,k)=one/(one+fv*qgues(i,j,k)) ! factor for tv to tsen conversion
end do
end do
end do
! Compute saturation specific humidity.
iderivative = 0
if(qoption == 1)then
if(jiter == jiterstart)iderivative = 1
else
iderivative = 2
end if
ice=.true.
call genqsat(qsatg,ges_tsen(1,1,1,ntguessig),ges_prsl(1,1,1,ntguessig),lat2,lon2, &
nsig,ice,iderivative)
!?????????????????????????? need any of this????
!! qoption 1: use psuedo-RH
! if(qoption==1)then
!
!
!! qoption 2: use normalized RH
! else
if(qoption == 2) then
! Load arrays based on option for moisture background error
! variance update for anisotropic mode
if( anisotropic .and. .not.rtma_subdomain_option ) then
hswgtsum=sum(hswgt(1:ngauss))
call setup_sub2fslab
do k=1,nsig
do j=1,lon2
do i=1,lat2
rhgues(i,j,k)=qgues(i,j,k)/qsatg(i,j,k)
end do
end do
end do
if( regional ) then
allocate(rh0f(pf2aP1%nlatf,pf2aP1%nlonf,nsig1o))
call sub2fslab(rhgues,rh0f)
do k=indices%kps,indices%kpe
ivar=idvar(k)
if(ivar==nrf3_loc(nrf3_q)) then
kvar=k-kvar_start(ivar)+1
do k1=1,nsig1o
if(levs_id(k1)==kvar) exit
end do
do j=indices%jps,indices%jpe
do i=indices%ips,indices%ipe
l =max(min(int(rllatf(i,j)),mlat),1)
l2=min((l+1),mlat)
dl2=rllatf(i,j)-float(l)
dl1=one-dl2
factk=dl1*corz(l,kvar,nrf3_q)+dl2*corz(l2,kvar,nrf3_q)
call fact_qopt2(factk,rh0f(i,j,k1),kvar)
do igauss=1,ngauss
factor=hswgt(igauss)*factk*an_amp0(ivar)/sqrt(hswgtsum)
filter_all(1)%amp(igauss,i,j,k)=factor*filter_all(2)%amp(igauss,i,j,k)
if (allocated(ensamp)) then
filter_all(1)%amp(igauss,i,j,k)=filter_all(1)%amp(igauss,i,j,k)*ensamp(i,j,k1)
end if
end do
end do
end do
end if
end do
deallocate(rh0f)
else
allocate(rh0f(pf2aP1%nlatf,pf2aP1%nlonf,nsig1o))
allocate(rh2f(pf2aP2%nlatf,pf2aP2%nlonf,nsig1o))
allocate(rh3f(pf2aP3%nlatf,pf2aP3%nlonf,nsig1o))
call sub2fslab_glb (rhgues,rh0f,rh2f,rh3f)
do k=indices%kps,indices%kpe
ivar=idvar(k)
if(ivar==nrf3_loc(nrf3_q)) then
kvar=k-kvar_start(ivar)+1
do k1=1,nsig1o
if(levs_id(k1)==kvar) exit
end do
! zonal patch
do j=indices%jps,indices%jpe
do i=indices%ips,indices%ipe
call get_stat_factk(p0ilatf(i),ivar,kvar,factk, &
rh0f(i,j,k1),one)
do igauss=1,ngauss
factor=hswgt(igauss)*factk*an_amp0(ivar)/sqrt(hswgtsum)
filter_all(1)%amp(igauss,i,j,k)=factor*filter_all(2)%amp(igauss,i,j,k)
if (allocated(ensamp0f)) then
filter_all(1)%amp(igauss,i,j,k)=filter_all(1)%amp(igauss,i,j,k)*ensamp0f(i,j,k1)
end if
end do
end do
end do
! polar patches
do j=indices_p%jps,indices_p%jpe
do i=indices_p%ips,indices_p%ipe
! north polar
if(p2ilatf(i,j)/=zero) then
call get_stat_factk(p2ilatf(i,j),ivar,kvar,factk, &
rh2f(i,j,k1),one)
else
call get_stat_factk(p2ilatfm ,ivar,kvar,factk, &
rh2f(i,j,k1),one)
end if
do igauss=1,ngauss
factor=hswgt(igauss)*factk*an_amp0(ivar)/sqrt(hswgtsum)
filter_p2(1)%amp(igauss,i,j,k)=factor*filter_p2(2)%amp(igauss,i,j,k)
if(allocated(ensamp2f)) then
filter_p2(1)%amp(igauss,i,j,k)=filter_p2(1)%amp(igauss,i,j,k)*sqrt(ensamp2f(i,j,k))
end if
end do
! south polar
if(p3ilatf(i,j)/=zero) then
call get_stat_factk(p3ilatf(i,j),ivar,kvar,factk, &
rh3f(i,j,k1),one)
else
call get_stat_factk(p3ilatfm ,ivar,kvar,factk, &
rh3f(i,j,k1),one)
end if
do igauss=1,ngauss
factor=factk*an_amp0(ivar)/sqrt(real(ngauss,r_kind))
filter_p3(1)%amp(igauss,i,j,k)=factor*filter_p3(2)%amp(igauss,i,j,k)
if(allocated(ensamp3f)) then
filter_p3(1)%amp(igauss,i,j,k)=filter_p3(1)%amp(igauss,i,j,k)*sqrt(ensamp3f(i,j,k))
end if
end do
end do
end do
end if
end do
deallocate(rh0f,rh2f,rh3f)
end if
call destroy_sub2fslab
end if
! End of qoption block
endif
call q_diag(mype)
! get subdomain latitude array
j = mype + 1
do i = 1, lat2
n = min(max(1, istart(j)+i-2), nlat)
xlats(i) = rlats(n)
enddo
! get corresponding pressure interface
do k = 1, nsig+1
do j = 1, lon2
do i = 1, lat2
prsi(i,j,k) = ges_prsi(i,j,k,ntguessig)
enddo
enddo
enddo
! import gfs co2 data
iyear = iadate(1)
month = iadate(2)
call read_gfsco2 &
! --- inputs:
( iyear,month,igfsco2,xlats,prsi,lat2,lon2,nsig,mype, &
! --- outputs:
ges_co2 )
! End of routine
return
end subroutine compute_derived