!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS !
!BOP -------------------------------------------------------------------
!
! !MODULE: m_berror_stats_reg - a module of berror_stats input
!
! !DESCRIPTION:
!
! !INTERFACE:
module m_berror_stats_reg
use kinds,only : i_kind,r_kind
use constants, only: zero,one,max_varname_length,half
use gridmod, only: nsig
use chemmod, only : berror_chem,upper2lower,lower2upper
use m_berror_stats, only: usenewgfsberror,berror_stats
implicit none
private ! except
! reconfigurable parameters, via NAMELIST/setup/
public :: berror_stats ! reconfigurable filename
! interfaces to file berror_stats.
public :: berror_set_reg ! set internal parameters
public :: berror_get_dims_reg ! get dimensions, jfunc::createj_func()
public :: berror_read_bal_reg ! get cross-cov.stats., balmod::prebal()
public :: berror_read_wgt_reg ! get auto-cov.stats., prewgt()
! !REVISION HISTORY:
! 25Feb10 - Zhu - adopt code format from m_berror_stats
! - extract from rdgstat_reg
! - change sructure of error file
! - make changes for generalized control variables
! 01Oct15 Zhu - use centralized cloud_names_fwd and n_clouds_fwd to add
! flexibility for clouds (either cw or individual hydrometeros)
! variances/correlations for all-sky radiance assimilation
! 24Jun16 - Guo - replaced the local berror_stats, with a global user
! configurable m_berror_stats::berror_stats, for the
! filename. This ensure the consistency, as well as the
! reconfigurability of this variable through the GSI.
!EOP ___________________________________________________________________
character(len=*),parameter :: myname='m_berror_stats_reg'
! The same default filename is used as in m_berror_stats, to take the
! advantage of the same berror_stats= entry in the /setup/ namelist, to
! override the default value. Otherwise, a special entry would have to
! be defined for this module in the /setup/ namelist.
integer(i_kind),parameter :: default_unit_ = 22
integer(i_kind),parameter :: ERRCODE=2
logical,save:: cwcoveqqcov_
integer(i_kind),allocatable,dimension(:):: lsig
real(r_kind),allocatable,dimension(:):: coef1,coef2
contains
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: berror_get_dims_reg - get dimensions
!
! !DESCRIPTION:
!
! !INTERFACE:
subroutine berror_get_dims_reg(msig,mlat,unit)
implicit none
integer(i_kind) ,intent( out) :: msig ! dimension of levels
integer(i_kind) ,intent( out) :: mlat ! dimension of latitudes
integer(i_kind),optional,intent(in ) :: unit ! logical unit [22]
! !REVISION HISTORY:
!EOP ___________________________________________________________________
character(len=*),parameter :: myname_=myname//'::berror_get_dims_reg'
integer(i_kind) :: inerr
! Read dimension of stats file
inerr=default_unit_
if(present(unit)) inerr = unit
open(inerr,file=berror_stats,form='unformatted',status='old')
rewind inerr
read(inerr) msig,mlat
close(inerr)
end subroutine berror_get_dims_reg
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: berror_set_reg - set (logical) parameter options internal to B
!
! !DESCRIPTION:
!
! !INTERFACE:
subroutine berror_set_reg(opt,value)
implicit none
character(len=*),intent(in) :: opt
logical(i_kind), intent(in) :: value
! !REVISION HISTORY:
! 2014-10-15 - Zhu - adopted from m_berror_stat to make code structure similar
!EOP ___________________________________________________________________
character(len=*),parameter :: myname_=myname//'::berror_set_reg'
logical found
found=.false.
if(trim(opt)=='cwcoveqqcov') then
cwcoveqqcov_=value
found=.true.
endif
if(.not.found) then
write(6,*) myname_,'(PREBAL_reg): ***ERROR*** cannot find:', trim(opt)
call stop2(999)
endif
end subroutine berror_set_reg
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: berror_read_bal_reg - get cross-corr. coefficients
!
! !DESCRIPTION:
!
! !INTERFACE:
subroutine berror_read_bal_reg(msig,mlat,agvi,bvi,wgvi,mype,unit)
use kinds,only : r_single
use guess_grids, only: ges_psfcavg,ges_prslavg
implicit none
integer(i_kind) ,intent(in ) :: msig,mlat
real(r_kind),dimension(0:mlat+1,nsig,nsig),intent( out) :: agvi
real(r_kind),dimension(0:mlat+1,nsig) ,intent( out) :: wgvi
real(r_kind),dimension(0:mlat+1,nsig) ,intent( out) :: bvi
integer(i_kind) ,intent(in ) :: mype ! "my" processor ID
integer(i_kind),optional ,intent(in ) :: unit ! an alternative unit
! !REVISION HISTORY:
! 25Feb10 - Zhu - extracted from rdgstat_reg
! - change the structure of error file
! - make changes for generalized control variables
!EOP ___________________________________________________________________
character(len=*),parameter :: myname_=myname//'::berror_read_bal_reg'
! workspaces/variables for data not returned
integer(i_kind) k,i,m,j,m1,l1,l
integer(i_kind):: nsigstat,nlatstat
integer(i_kind):: inerr
real(r_kind),dimension(nsig) :: rlsig
real(r_single),dimension(:),allocatable:: clat_avn,sigma_avn
real(r_single),dimension(:,:),allocatable:: bv_avn,wgv_avn
real(r_single),dimension(:,:,:),allocatable:: agv_avn
real(r_kind),dimension(:),allocatable:: rlsigo
! Open background error statistics file
inerr=default_unit_
if(present(unit)) inerr=unit
open(inerr,file=berror_stats,form='unformatted',status='old')
! Read header.
rewind inerr
read(inerr) nsigstat,nlatstat
if(mype==0) then
write(6,*) myname_,'(PREBAL_REG): get balance variables', &
'"',trim(berror_stats),'". ', &
'mype,nsigstat,nlatstat =', &
mype,nsigstat,nlatstat
end if
allocate ( clat_avn(mlat) )
allocate ( sigma_avn(1:msig) )
allocate ( rlsigo(1:msig) )
if(usenewgfsberror)then
allocate ( agv_avn(mlat,1:msig,1:msig) )
allocate ( bv_avn(mlat,1:msig),wgv_avn(mlat,1:msig) )
else
allocate ( agv_avn(0:mlat+1,1:msig,1:msig) )
allocate ( bv_avn(0:mlat+1,1:msig),wgv_avn(0:mlat+1,1:msig) )
end if
! Read background error file to get balance variables
read(inerr)clat_avn,(sigma_avn(k),k=1,msig)
read(inerr)agv_avn,bv_avn,wgv_avn
close(inerr)
! compute vertical(pressure) interpolation index and weight
do k=1,nsig
rlsig(k)=log(ges_prslavg(k)/ges_psfcavg)
enddo
do k=1,msig
rlsigo(k)=log(sigma_avn(k))
enddo
allocate(lsig(nsig),coef1(nsig),coef2(nsig))
do k=1,nsig
if(rlsig(k)>=rlsigo(1))then
m=1
m1=2
lsig(k)=1
coef1(k)=one
coef2(k)=zero
else if(rlsig(k)<=rlsigo(msig))then
m=msig-1
m1=msig
lsig(k)=msig-1
coef1(k)=zero
coef2(k)=one
else
m_loop: do m=1,msig-1
m1=m+1
if((rlsig(k)<=rlsigo(m)) .and. &
(rlsig(k)>rlsigo(m1)) )then
lsig(k)=m
exit m_loop
end if
enddo m_loop
coef1(k)=(rlsigo(m1)-rlsig(k))/(rlsigo(m1)-rlsigo(m))
coef2(k)=one-coef1(k)
if(lsig(k)==msig)then
lsig(k)=msig-1
coef2(k)=one
coef1(k)=zero
endif
endif
end do
! Load agv wgv bv
do k=1,nsig
m=lsig(k)
m1=m+1
do i=1,mlat
wgvi(i,k)=wgv_avn(i,m)*coef1(k)+wgv_avn(i,m1)*coef2(k)
bvi (i,k)=bv_avn (i,m)*coef1(k)+bv_avn (i,m1)*coef2(k)
enddo
do j=1,nsig
l=lsig(j)
l1=l+1
do i=1,mlat
agvi(i,j,k)=(agv_avn(i,l,m)*coef1(j)+agv_avn(i,l1,m)*coef2(j))*coef1(k) &
+(agv_avn(i,l,m1)*coef1(j)+agv_avn(i,l1,m1)*coef2(j))*coef2(k)
enddo
enddo
enddo
deallocate (agv_avn,bv_avn,wgv_avn,clat_avn,sigma_avn,rlsigo)
agvi(0,:,:)=agvi(1,:,:)
wgvi(0,:)=wgvi(1,:)
bvi(0,:)=bvi(1,:)
agvi(mlat+1,:,:)=agvi(mlat,:,:)
wgvi(mlat+1,:)=wgvi(mlat,:)
bvi(mlat+1,:)=bvi(mlat,:)
return
end subroutine berror_read_bal_reg
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: berror_read_wgt_reg - read auto-corr. coeffs.
!
! !DESCRIPTION:
!
! !INTERFACE:
subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qoption,varcw,cwoption,mype,unit)
use kinds,only : r_single,r_kind
use gridmod,only : nsig, twodvar_regional
use gsi_io, only : verbose
use control_vectors,only: nrf,nc2d,nc3d,mvars,nvars
use control_vectors,only: cvars => nrf_var
use control_vectors,only: cvars2d,cvars3d,cvarsmd
use guess_grids, only: ges_psfcavg,ges_prslavg
use constants, only: zero,one,ten,three
use mpeu_util,only: getindex
use radiance_mod, only: icloud_cv,n_clouds_fwd,cloud_names_fwd
implicit none
integer(i_kind) ,intent(in ) :: qoption
integer(i_kind) ,intent(in ) :: cwoption
integer(i_kind) ,intent(in ) :: msig,mlat
integer(i_kind) ,intent(in ) :: mype ! "my" processor ID
integer(i_kind),optional ,intent(in ) :: unit ! an alternative unit
real(r_kind),dimension(:,:,:),intent(inout):: corz
real(r_kind),dimension(:,:) ,intent(inout) :: corp
real(r_kind),dimension(0:mlat+1,1:nsig,1:nc3d), intent(inout):: hwll
real(r_kind),dimension(0:mlat+1,nvars-nc3d) , intent(inout):: hwllp
real(r_kind),dimension(nsig,0:mlat+1,1:nc3d),intent(inout):: vz
real(r_kind),dimension(mlat,nsig),intent(inout)::varq
real(r_kind),dimension(mlat,nsig),intent(inout)::varcw
real(r_kind),dimension(nsig),intent(out):: rlsig
! !REVISION HISTORY:
! 25Feb10 - Zhu - extract from rdgstat_reg
! - change the structure of background error file
! - make changes for generalizing control variables
! - move varq,factoz here from prewgt_reg
! 28May10 Todling Obtain variable id's on the fly (add getindex)
! 01Jun10 Todling These are now alloctable: corz,corp,hwll,hwllp,vz
! 22Jun10 Treadon - move nrf3_loc and nrf2_loc allocate outside read loop
! 23Jun10 Treadon - explicitly specify dimensions for hwll,hwllp,vz
! 20Nov10 Pagowski - make var name longer for chemical berror and
! related change in read
! 16Feb11 Zhu - add gust,vis,pblh
! 15Dec12 Zhu - add varcw and cwoption for all-sky radiance assimiation
! 03Feb14 Todling - varq & qoption in arg list (remove dep on jfunc)
! 19Mar14 pondeca - add wspd10m
! 10Apr14 pondeca - add td2m,mxtm,mitm,pmsl
! 07May14 pondeca - add howv
! 10Jun14 Zhu - tune error variance and correlation lengths of cw for
! all-sky radiance assimilation
! 19Jun14 carley/zhu - add tcamt and lcbas
! 10Jul15 pondeca - add cldch
! 01Oct15 Zhu - use centralized cloud_names_fwd and n_clouds_fwd to add
! flexibility for clouds (either cw or individual hydrometeros)
! variances/correlations for all-sky radiance assimilation
! 05May16 pondeca - add uwnd10m, vuwn10m
! 29Aug16 stelios - Update the constants for howv #ww3
! 2016-09-xx CAPS(G. Zhao)
! - add same error stats as q for hydrometeor variables,
! - which can be tuned in sub prewgt_reg.f90.
! - extra caution required for using logarithmic transformation
! 2018-10-22 CAPS(C.Liu) - add w
!
!EOP ___________________________________________________________________
character(len=*),parameter :: myname_=myname//'::berror_read_wgt_reg'
real(r_kind),parameter:: r25 = one/25.0_r_kind
real(r_kind),parameter:: zero_3 = 0.3_r_kind
real(r_kind),parameter:: vz_oz = 0.53333333_r_kind
! workspace variables not returned
real(r_single),dimension(:,:),allocatable:: corqq_avn
real(r_single),dimension(:,:),allocatable:: corz_avn,hwll_avn,vztdq_avn
real(r_single),dimension(1:mlat,msig,nrf):: corz_tmp
real(r_single),dimension(0:mlat+1,msig,nrf):: hwll_tmp
real(r_single),dimension(msig,0:mlat+1,nrf):: vz_tmp
character*5 :: varshort
character(len=max_varname_length) :: var
logical,dimension(nrf):: nrf_err
integer(i_kind) :: nrf3_oz,nrf3_q,nrf3_cw,nrf3_sf,nrf3_vp,nrf2_sst
integer(i_kind) :: nrf3_t,nrf2_gust,nrf2_vis,nrf2_pblh,nrf2_ps,nrf2_wspd10m
integer(i_kind) :: nrf2_td2m,nrf2_mxtm,nrf2_mitm,nrf2_pmsl,nrf2_howv,nrf2_tcamt,nrf2_lcbas,nrf2_cldch
integer(i_kind) :: nrf2_uwnd10m,nrf2_vwnd10m
integer(i_kind) :: nrf3_sfwter,nrf3_vpwter
integer(i_kind) :: nrf3_dbz
integer(i_kind) :: nrf3_ql,nrf3_qi,nrf3_qr,nrf3_qs,nrf3_qg,nrf3_qnr,nrf3_w
integer(i_kind) :: inerr,istat
integer(i_kind) :: nsigstat,nlatstat,isig
integer(i_kind) :: loc,m1,m,i,n,j,k,ivar,ic
integer(i_kind),allocatable,dimension(:) :: nrf2_loc,nrf3_loc,nmotl_loc
real(r_kind) :: factoz
real(r_kind) :: raux
real(r_kind),dimension(nsig):: dlsig
! corz = sqrt(corz)
real(r_kind), parameter :: corz_default=one,hwll_default=100000_r_kind,vz_default=one
!real(r_kind), parameter :: corz_default=one,hwll_default=27000.00000000,vz_default=one*10
logical :: print_verbose
real(r_kind) ,dimension(mlat,1,2) :: cov_dum
! character(256) :: filename
! filename = 'howv_var_berr.bin'
print_verbose=.false.
if(verbose)print_verbose=.true.
do k=1,nsig
rlsig(k)=log(ges_prslavg(k)/ges_psfcavg)
enddo
! Open background error statistics file
inerr=default_unit_
if(present(unit)) inerr=unit
open(inerr,file=berror_stats,form='unformatted',status='old')
! Read header.
rewind inerr
read(inerr) nsigstat,nlatstat
if(mype==0) then
write(6,*) myname_,'(PREWGT_REG): read error amplitudes ', &
'"',trim(berror_stats),'". ', &
'mype,nsigstat,nlatstat =', &
mype,nsigstat,nlatstat
end if
! Read background error file to get past alance variables
read(inerr)
read(inerr)
allocate(nrf3_loc(nc3d),nrf2_loc(nc2d),nmotl_loc(mvars))
do n=1,nc3d
nrf3_loc(n)=getindex(cvars,cvars3d(n))
enddo
do n=1,nc2d
nrf2_loc(n)=getindex(cvars,cvars2d(n))
enddo
do n=1,mvars
nmotl_loc(n)=getindex(cvars,cvarsmd(n))
enddo
! Read amplitudes
nrf_err=.false.
read: do
if (berror_chem) then
read(inerr,iostat=istat) varshort,isig
var=upper2lower(varshort)
if (trim(var) == 'pm25') var = 'pm2_5'
else
read(inerr,iostat=istat) varshort, isig
var=varshort
endif
if (istat /= 0) exit read
do n=1,nrf
if (trim(var)==cvars(n)) then
nrf_err(n)=.true.
loc=n
exit
else
loc=-999
end if
end do
if(loc == -999)then
if(mype == 0)write(6,*) 'variable in input file, but not used in analysis ',var,isig
read(inerr)
read(inerr)
if(isig > 1)read(inerr)
cycle read
end if
allocate ( corz_avn(1:mlat,1:isig) )
if (trim(var)/='q' .or. (trim(var)=='cw' .and. cwoption==2)) then
read(inerr) corz_avn
else
allocate ( corqq_avn(1:mlat,1:isig) )
read(inerr) corz_avn,corqq_avn
end if
if(usenewgfsberror)then
allocate ( hwll_avn(mlat,1:isig) )
else
allocate ( hwll_avn(0:mlat+1,1:isig) )
end if
read(inerr) hwll_avn
if (isig==msig) then
if(usenewgfsberror)then
allocate ( vztdq_avn(mlat,1:isig) )
else
allocate ( vztdq_avn(1:isig,0:mlat+1) )
end if
read(inerr) vztdq_avn
do n=1,nc3d
if (nrf3_loc(n)==loc) then
if ((trim(var)=='q' .and. qoption==2) .or. (trim(var)=='cw' .and. cwoption==2)) then
! choose which q stat to use
do k=1,msig
do i=1,mlat
corz_tmp(i,k,n)=corqq_avn(i,k)
end do
end do
else
do k=1,msig
do i=1,mlat
corz_tmp(i,k,n)=corz_avn(i,k)
end do
end do
end if
do k=1,msig
if(usenewgfsberror)then
do i=1,mlat
hwll_tmp(i,k,n)=hwll_avn(i,k)
vz_tmp(k,i,n)=vztdq_avn(i,k)
end do
hwll_tmp(0,k,n)=hwll_avn(1,k)
hwll_tmp(mlat+1,k,n)=hwll_avn(mlat,k)
vz_tmp(k,0,n)=vztdq_avn(1,k)
vz_tmp(k,mlat+1,n)=vztdq_avn(mlat,k)
else
do i=0,mlat+1
hwll_tmp(i,k,n)=hwll_avn(i,k)
vz_tmp(k,i,n)=vztdq_avn(k,i)
end do
end if
end do
exit
end if
end do
deallocate ( vztdq_avn )
else if(isig == 1)then
do n=1,nc2d
if (nrf2_loc(n)==loc) then
do i=1,mlat
corp(i,n)=corz_avn(i,1)
hwllp(i,n)=hwll_avn(i,1)
end do
if(usenewgfsberror)then
hwllp(0,n)=hwll_avn(1,1)
hwllp(mlat+1,n)=hwll_avn(mlat,1)
else
hwllp(0,n)=hwll_avn(0,1)
hwllp(mlat+1,n)=hwll_avn(mlat+1,1)
end if
exit
end if
end do
end if
deallocate ( corz_avn )
deallocate ( hwll_avn )
if(allocated(corqq_avn)) deallocate ( corqq_avn )
enddo read
close(inerr)
! 3d variable
do n=1,nc3d
loc=nrf3_loc(n)
if (loc <= 0)cycle
if (nrf_err(loc)) then
do k=1,nsig
m=lsig(k)
m1=m+1
do i=1,mlat
corz(i,k,n)=corz_tmp(i,m,n)*coef1(k)+corz_tmp(i,m1,n)*coef2(k)
enddo
do i=0,mlat+1
hwll(i,k,n)=hwll_tmp(i,m,n)*coef1(k)+hwll_tmp(i,m1,n)*coef2(k)
vz(k,i,n)=vz_tmp(m,i,n)*coef1(k)+vz_tmp(m1,i,n)*coef2(k)
enddo
enddo
else
do k=1,nsig
do i=1,mlat
corz(i,k,n)=corz_default
enddo
do i=0,mlat+1
hwll(i,k,n)=hwll_default
vz(k,i,n)=vz_default
enddo
enddo
if(mype==0) then
write(6,*)'Assigned default statistics to variable ',cvars(loc)
endif
end if
enddo
! Get control variable indexes
nrf3_t =getindex(cvars3d,'t')
nrf3_q =getindex(cvars3d,'q')
nrf3_oz =getindex(cvars3d,'oz')
nrf3_cw =getindex(cvars3d,'cw')
nrf3_sf =getindex(cvars3d,'sf')
nrf3_vp =getindex(cvars3d,'vp')
nrf3_dbz=getindex(cvars3d,'dbz')
nrf2_sst=getindex(cvars2d,'sst')
nrf2_gust=getindex(cvars2d,'gust')
nrf2_vis=getindex(cvars2d,'vis')
nrf2_pblh=getindex(cvars2d,'pblh')
nrf2_ps=getindex(cvars2d,'ps')
nrf2_wspd10m=getindex(cvars2d,'wspd10m')
nrf2_td2m=getindex(cvars2d,'td2m')
nrf2_mxtm=getindex(cvars2d,'mxtm')
nrf2_mitm=getindex(cvars2d,'mitm')
nrf2_pmsl=getindex(cvars2d,'pmsl')
nrf2_howv=getindex(cvars2d,'howv')
nrf3_sfwter =getindex(cvars3d,'sfwter')
nrf3_vpwter =getindex(cvars3d,'vpwter')
nrf2_tcamt=getindex(cvars2d,'tcamt')
nrf2_lcbas=getindex(cvars2d,'lcbas')
nrf2_cldch=getindex(cvars2d,'cldch')
nrf2_uwnd10m=getindex(cvars2d,'uwnd10m')
nrf2_vwnd10m=getindex(cvars2d,'vwnd10m')
! W, Cloud and hydrometeor fields
nrf3_ql =getindex(cvars3d,'ql')
nrf3_qi =getindex(cvars3d,'qi')
nrf3_qr =getindex(cvars3d,'qr')
nrf3_qs =getindex(cvars3d,'qs')
nrf3_qg =getindex(cvars3d,'qg')
nrf3_qnr =getindex(cvars3d,'qnr')
nrf3_w =getindex(cvars3d,'w')
if(nrf3_q>0 .and. qoption==2)then
do k=1,nsig
do j=1,mlat
varq(j,k)=min(max(real(corz(j,k,nrf3_q),r_kind),0.0015_r_kind),one)
corz(j,k,nrf3_q)=one
enddo
enddo
endif
if( nrf3_dbz>0 )then
if(.not. nrf3_t>0) then
write(6,*)'not as expect,stop'
stop
endif
corz(:,:,nrf3_dbz)=10.0_r_kind
hwll(:,:,nrf3_dbz)=hwll(:,:,nrf3_t)
vz(:,:,nrf3_dbz)=vz(:,:,nrf3_t)
endif
if (nrf3_oz>0) then
factoz = 0.0002_r_kind*r25
corz(:,:,nrf3_oz)=factoz
! hwll:,:,nrf3_oz)=400000.0_r_kind
vz(:,:,nrf3_oz)=vz_oz
end if
if(nrf3_cw > 0)then
if (cwcoveqqcov_ ) then
corz(:,:,nrf3_cw)=corz(:,:,nrf3_q)
hwll(:,:,nrf3_cw)=hwll(:,:,nrf3_q)
vz(:,:,nrf3_cw)=vz(:,:,nrf3_q)
else
corz(:,:,nrf3_cw)=zero
if (cwoption==2) then
do k=1,nsig
if (ges_prslavg(k)>15.0_r_kind) then
do j=1,mlat
varcw(j,k)=max(real(corz(j,k,nrf3_cw),r_kind),zero)
corz(j,k,nrf3_cw)=one
enddo
end if
enddo
end if
if (cwoption==1 .or. cwoption==3) then
do k=1,nsig
if (ges_prslavg(k)>15.0_r_kind) then
do j=1,mlat
corz(j,k,nrf3_cw)=one
end do
end if
end do
hwll(:,:,nrf3_cw)=0.5_r_kind*hwll(:,:,nrf3_q)
vz(:,:,nrf3_cw)=0.5_r_kind*vz(:,:,nrf3_q)
end if
end if
else if (icloud_cv .and. n_clouds_fwd>0 .and. cwoption==3) then
do n=1,size(cvars3d)
do ic=1,n_clouds_fwd
if(trim(cvars3d(n))==trim(cloud_names_fwd(ic))) then
ivar=n
do k=1,nsig
do i=1,mlat
corz(i,k,ivar)=one
end do
end do
hwll(:,:,ivar)=0.5_r_kind*hwll(:,:,nrf3_q)
vz (:,:,ivar)=0.5_r_kind*vz (:,:,nrf3_q)
exit
end if
end do
end do
end if ! n_clouds_fwd>0 .and. nrf3_cw<=0 .and. cwoption==3
if (nrf3_sfwter>0) then
if(mype==0) write(6,*)'Replace default with appropriate statistics for variable sfwter'
corz(1:mlat,1:nsig,nrf3_sfwter)=corz(1:mlat,1:nsig,nrf3_sf)
hwll(0:mlat+1,1:nsig,nrf3_sfwter)=hwll(0:mlat+1,1:nsig,nrf3_sf)
endif
if (nrf3_vpwter>0) then
if(mype==0) write(6,*)'Replace default with appropriate statistics for variable vpwter'
corz(1:mlat,1:nsig,nrf3_vpwter)=corz(1:mlat,1:nsig,nrf3_vp)
hwll(0:mlat+1,1:nsig,nrf3_vpwter)=hwll(0:mlat+1,1:nsig,nrf3_vp)
endif
! W and Cloud/hydrometeor variables share same error structure with humidity q by default
! (not with log transformed cloud variables)
if (nrf3_ql>0) then
if(mype==0) &
write(6,*)myname_,"@pe=",mype," set b_error stats of ql = qv(moisture) and nrf3_ql=",nrf3_ql
corz(:,:,nrf3_ql) = corz(:,:,nrf3_q)
hwll(:,:,nrf3_ql) = hwll(:,:,nrf3_q)
vz(:,:, nrf3_ql) = vz(:,:, nrf3_q)
end if
if (nrf3_qi>0) then
if(mype==0) &
write(6,*)myname_,"@pe=",mype," set b_error stats of qi = qv(moisture) and nrf3_qi=",nrf3_qi
corz(:,:,nrf3_qi) = corz(:,:,nrf3_q)
hwll(:,:,nrf3_qi) = hwll(:,:,nrf3_q)
vz(:,:, nrf3_qi) = vz(:,:, nrf3_q)
end if
if (nrf3_qr>0) then
if(mype==0) &
write(6,*)myname_,"@pe=",mype," set b_error stats of qr = qv(moisture) and nrf3_qr=",nrf3_qr
corz(:,:,nrf3_qr) = corz(:,:,nrf3_q)
hwll(:,:,nrf3_qr) = hwll(:,:,nrf3_q)
vz(:,:, nrf3_qr) = vz(:,:, nrf3_q)
end if
if (nrf3_qs>0) then
if(mype==0) &
write(6,*)myname_,"@pe=",mype," set b_error stats of qs = qv(moisture) and nrf3_qs=",nrf3_qs
corz(:,:,nrf3_qs) = corz(:,:,nrf3_q)
hwll(:,:,nrf3_qs) = hwll(:,:,nrf3_q)
vz(:,:, nrf3_qs) = vz(:,:, nrf3_q)
end if
if (nrf3_qg>0) then
if(mype==0) &
write(6,*)myname_,"@pe=",mype," set b_error stats of qg = qv(moisture) and nrf3_qg=",nrf3_qg
corz(:,:,nrf3_qg) = corz(:,:,nrf3_q)
hwll(:,:,nrf3_qg) = hwll(:,:,nrf3_q)
vz(:,:, nrf3_qg) = vz(:,:, nrf3_q)
end if
if (nrf3_qnr>0) then
if(mype==0) &
write(6,*)myname_,"@pe=",mype," set b_error stats of qnr = qv(moisture) and nrf3_nr=",nrf3_qnr
corz(:,:,nrf3_qnr) = corz(:,:,nrf3_q)
hwll(:,:,nrf3_qnr) = hwll(:,:,nrf3_q)
vz(:,:, nrf3_qnr) = vz(:,:, nrf3_q)
end if
if (nrf3_w>0) then
if(mype==0) &
write(6,*)myname_,"@pe=",mype," set b_error stats of w = qv(moisture) and nrf3_w=",nrf3_w
corz(:,:,nrf3_w) = corz(:,:,nrf3_q)
hwll(:,:,nrf3_w) = hwll(:,:,nrf3_q)
vz(:,:, nrf3_w) = vz(:,:, nrf3_q)
end if
! 2d variable
do n=1,nc2d
loc=nrf2_loc(n)
if(loc <= 0)cycle
! If we want to use the sst in global file uncomment the following if statement
if (n==nrf2_sst) then
! if(.not. usenewgfsberror)then
do i=1,mlat
corp(i,n)=zero_3
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_sf)
end do
! end if
else if (n==nrf2_gust) then
do i=1,mlat
corp(i,n)=three
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_q)
end do
else if (n==nrf2_vis) then
do i=1,mlat
corp(i,n)=3.0_r_kind
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_t)
end do
else if (n==nrf2_pblh) then
do i=1,mlat
corp(i,n)=500.0_r_kind
end do
do i=0,mlat+1
hwllp(i,n)=three*hwll(i,1,nrf3_t)
end do
else if (n==nrf2_wspd10m) then
do i=1,mlat
corp(i,n)=three
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_q)
end do
else if (n==nrf2_td2m) then
raux=maxval(corz(1:mlat,1,nrf3_q))
do i=1,mlat
corp(i,n)=(corz(i,1,nrf3_q)/raux)*three !tentatively
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_q) !tentatively
end do
else if (n==nrf2_mxtm) then
do i=1,mlat
corp(i,n)=corz(i,1,nrf3_t)
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_t)
end do
else if (n==nrf2_mitm) then
do i=1,mlat
corp(i,n)=corz(i,1,nrf3_t)
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_t)
end do
else if (n==nrf2_pmsl) then
do i=1,mlat
corp(i,n)=corp(i,nrf2_ps)
end do
do i=0,mlat+1
hwllp(i,n)=hwllp(i,nrf2_ps)
end do
else if (n==nrf2_howv) then
call read_howv_stats(mlat,1,2,cov_dum)
do i=1,mlat
corp(i,n)=cov_dum(i,1,1) !#ww3
hwllp(i,n) = cov_dum(i,1,2)
end do
hwllp(0,n) = hwllp(1,n)
hwllp(mlat+1,n) = hwllp(mlat,n)
if (mype==0) print*, 'corp(i,n) = ', corp(:,n)
if (mype==0) print*, ' hwllp(i,n) = ', hwllp(:,n)
! corp(:,n)=cov_dum(:,1)
!do i=1,mlat
! corp(i,n)=0.4_r_kind !#ww3
!end do
! do i=0,mlat+1
! hwllp(i,n)=hwll(i,1,nrf3_sf) !tentatively !#ww3 hwllp(i,n)=150000_r_kind !
! end do
else if (n==nrf2_tcamt) then
do i=1,mlat
corp(i,n)=50.0_r_kind
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_t)
end do
else if (n==nrf2_lcbas) then
do i=1,mlat
corp(i,n)=40000.0_r_kind
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_t)
end do
if(print_verbose)print*, 'm_berror_reg: maxhwllp_lcbas=',maxval(hwllp(:,n))
else if (n==nrf2_cldch) then
do i=1,mlat
corp(i,n)=3.0_r_kind
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_t)
end do
if(print_verbose)print*, 'm_berror_reg: maxhwllp_cldch=',maxval(hwllp(:,n))
else if (n==nrf2_uwnd10m) then
do i=1,mlat
corp(i,n)=three
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_q)
end do
else if (n==nrf2_vwnd10m) then
do i=1,mlat
corp(i,n)=three
end do
do i=0,mlat+1
hwllp(i,n)=hwll(i,1,nrf3_q)
end do
end if
enddo
! motley variable
do n=1,mvars
loc=nmotl_loc(n)
if(loc <=0)cycle
loc=n+nc2d
if (cvarsmd(n)=='sti' .or. cvarsmd(n)=='stl') then
do i=1,mlat
corp(i,loc)=corz(i,1,nrf3_sf)
end do
do i=0,mlat+1
hwllp(i,loc)=hwll(i,1,nrf3_sf)
end do
else if (cvarsmd(n)=='pswter') then
do i=1,mlat
corp(i,loc)=corp(i,nrf2_ps)
end do
do i=0,mlat+1
hwllp(i,loc)=hwllp(i,nrf2_ps)
end do
else if (cvarsmd(n)=='twter') then
do i=1,mlat
corp(i,loc)=corz(i,1,nrf3_t)
end do
do i=0,mlat+1
hwllp(i,loc)=hwll(i,1,nrf3_t)
end do
else if (cvarsmd(n)=='qwter') then
do i=1,mlat
corp(i,loc)=corz(i,1,nrf3_q)
end do
do i=0,mlat+1
hwllp(i,loc)=hwll(i,1,nrf3_q)
end do
else if (cvarsmd(n)=='gustwter') then
do i=1,mlat
corp(i,loc)=corp(i,nrf2_gust)
end do
do i=0,mlat+1
hwllp(i,loc)=hwllp(i,nrf2_gust)
end do
else if (cvarsmd(n)=='wspd10mwter') then
do i=1,mlat
corp(i,loc)=corp(i,nrf2_wspd10m)
end do
do i=0,mlat+1
hwllp(i,loc)=hwllp(i,nrf2_wspd10m)
end do
else if (cvarsmd(n)=='td2mwter') then
do i=1,mlat
corp(i,loc)=corp(i,nrf2_td2m)
end do
do i=0,mlat+1
hwllp(i,loc)=hwllp(i,nrf2_td2m)
end do
else if (cvarsmd(n)=='mxtmwter') then
do i=1,mlat
corp(i,loc)=corp(i,nrf2_mxtm)
end do
do i=0,mlat+1
hwllp(i,loc)=hwllp(i,nrf2_mxtm)
end do
else if (cvarsmd(n)=='mitmwter') then
do i=1,mlat
corp(i,loc)=corp(i,nrf2_mitm)
end do
do i=0,mlat+1
hwllp(i,loc)=hwllp(i,nrf2_mitm)
end do
else if (cvarsmd(n)=='uwnd10mwter') then
do i=1,mlat
corp(i,loc)=corp(i,nrf2_uwnd10m)
end do
do i=0,mlat+1
hwllp(i,loc)=hwllp(i,nrf2_uwnd10m)
end do
else if (cvarsmd(n)=='vwnd10mwter') then
do i=1,mlat
corp(i,loc)=corp(i,nrf2_vwnd10m)
end do
do i=0,mlat+1
hwllp(i,loc)=hwllp(i,nrf2_vwnd10m)
end do
! if not found use default
else
do i=1,mlat
corp(i,loc)=corz_default
end do
do i=0,mlat+1
hwllp(i,loc)=hwll_default
end do
end if
enddo
deallocate(nrf3_loc,nrf2_loc,nmotl_loc)
! Normalize vz with del sigmma and convert to vertical grid units!
if(.not. twodvar_regional)then
dlsig(1)=rlsig(1)-rlsig(2)
do k=2,nsig-1
dlsig(k)=half*(rlsig(k-1)-rlsig(k+1))
enddo
dlsig(nsig)=rlsig(nsig-1)-rlsig(nsig)
do n=1,nc3d
do j=0,mlat+1
do k=1,nsig
vz(k,j,n)=vz(k,j,n)*dlsig(k)
end do
end do
end do
end if
return
end subroutine berror_read_wgt_reg
!++++
subroutine read_howv_stats(nlat,nlon,npar,arrout)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_howv_stats
! prgmmr: stelios org: ??? date: 2016-08-03
!
! abstract: improrts the bkerror stats for howv (variance, correlation lengths,
! etc). Each quantity has been calculated externaly, it is interpolated at
! the grid of URMA and it is saved as binary file. The exact format
! is shown in the code.
! The imported arrays can be 2d or 1d, aka can be lon/lat dependent or
! two one of the two (most probably on lat)
!
! For cross reference the MATLAB code is provided
! Export ::
! fid = fopen(['URMA_variance_Q',num2str(i1),'.bin'],'w');
! dum =squeeze(variance(i1,:,:));
! fwrite(fid,dum,'double');
! fclose(fid);
!
! Import ::
! fid = fopen(['URMA_variance_Q',num2str(i1),'1.bin'],'r');
! dum = fread(fid, 'double');
! A(:,:,i1) = reshape(dum,1597,2345);
! fclose(fid);
!
! To use it: 1. the data files (eg variance, correlation lengths have
! to be located at the running path. So the driving scripts have to
! copying the file of interest to the tmp dir used for e.g.:
!
! cp [...]/fix/urma2p5/URMA_variance_lat.bin URMA_variance1d.bin
!
! program history log:
! 2016-08-03 stelios
! 2016-08-26 stelios : Compatible with GSI.
! input argument list:
! filename - The name of the file
! output argument list:
! arr_out - One or Two dimensional field of quantity of interest
!
! attributes:
! language: f90
! machine: theia/gyre
!
!$$$ end documentation block
!
use kinds,only : r_kind, i_kind
!
implicit none
! Declare passed variables
integer(i_kind), intent(in )::nlat,nlon,npar
real(r_kind), dimension(nlat ,nlon, npar), intent( out)::arrout
! Declare local variables
integer(i_kind) :: reclength,i,j,i_npar
character(256) , parameter :: myname='read_how_stats : '
logical :: file_exists
integer(i_kind), parameter :: lun34=34
character(len=256),dimension(npar) :: filename
integer(i_kind), parameter :: dp1 = kind(1D0)
!
filename(1) = 'howv_var_berr.bin'
filename(2) = 'howv_lng_berr.bin'
!
arrout(:,:,1)=0.42_r_kind
arrout(:,:,2)=50000.0_r_kind
reclength=nlat*r_kind
!
do i_npar = 1,npar
inquire(file=trim(filename(i_npar)), exist=file_exists)
if (file_exists)then
open (unit=lun34 ,file=trim(filename(i_npar)) ,status='old' &
,form='unformatted' ,access='direct' ,recl=reclength )
do j = 1,nlon
read(unit=lun34 ,rec=j) (arrout(i,j,i_npar), i=1,nlat)
enddo
close(unit=lun34)
else
print*,myname, trim(filename(i_npar)) // ' does not exist'
end if
end do
end subroutine read_howv_stats
end module m_berror_stats_reg