!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! 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: izero,ione,zero,one use gridmod, only: nsig implicit none private ! except ! reconfigurable parameters, via NAMELIST/setup/ public :: berror_stats ! reconfigurable filename ! interfaces to file berror_stats. 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 !EOP ___________________________________________________________________ character(len=*),parameter :: myname='m_berror_stats_reg' ! Reconfigurable parameters, vai NAMELISt/setup/ character(len=256),save :: berror_stats = "berror_stats" ! filename integer(i_kind),parameter :: default_unit_ = 22_i_kind integer(i_kind),parameter :: ERRCODE=2_i_kind 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_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 gridmod,only : nlat,nlon use guess_grids, only: ges_psfcavg,ges_prslavg implicit none integer(i_kind) ,intent(in ) :: msig,mlat real(r_kind),dimension(0:mlat+ione,nsig,nsig),intent( out) :: agvi real(r_kind),dimension(0:mlat+ione,nsig) ,intent( out) :: wgvi real(r_kind),dimension(0:mlat+ione,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,n,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==izero) 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) ) allocate ( agv_avn(0:mlat+ione,1:msig,1:msig) ) allocate ( bv_avn(0:mlat+ione,1:msig),wgv_avn(0:mlat+ione,1:msig) ) ! 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=ione m1=2_i_kind lsig(k)=ione coef1(k)=one coef2(k)=zero else if(rlsig(k)<=rlsigo(msig))then m=msig-ione m1=msig lsig(k)=msig-ione coef1(k)=zero coef2(k)=one else m_loop: do m=1,msig-ione m1=m+ione 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-ione coef2(k)=one coef1(k)=zero endif endif end do ! Load agv wgv bv do k=1,nsig m=lsig(k) m1=m+ione 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+ione 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 agvi(0,:,:)=agvi(1,:,:) wgvi(0,:)=wgvi(1,:) bvi(0,:)=bvi(1,:) agvi(mlat+ione,:,:)=agvi(mlat,:,:) wgvi(mlat+ione,:)=wgvi(mlat,:) bvi(mlat+ione,:)=bvi(mlat,:) deallocate (agv_avn,bv_avn,wgv_avn,clat_avn,sigma_avn,rlsigo) 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,mype,unit) use kinds,only : r_single,r_kind use gridmod,only : nlat,nlon,nsig use control_vectors,only: nrf,nrf2,nrf3,nrf_var,nrf2_loc,nrf3_loc, & nrf3_oz,nrf3_q,nrf3_cw,nrf3_sf,nrf2_sst,nvars use jfunc,only: varq,qoption use guess_grids, only: ges_psfcavg,ges_prslavg use constants, only: ione,zero,one implicit none 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(1:mlat,nsig,nrf3),intent(inout):: corz real(r_kind),dimension(1:mlat,nrf2),intent(out):: corp real(r_kind),dimension(0:mlat+1,nsig,nrf3),intent(out):: hwll real(r_kind),dimension(0:mlat+1,nvars-nrf3),intent(out):: hwllp real(r_kind),dimension(nsig,0:mlat+1,1:nrf3),intent(out):: vz 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 !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:: clat_avn,sigma_avn real(r_single),dimension(:,:),allocatable:: bv_avn,wgv_avn,corqq_avn real(r_single),dimension(:,:,:),allocatable:: agv_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 var logical,dimension(nrf):: nrf_err integer(i_kind) :: inerr,istat integer(i_kind) :: nsigstat,nlatstat,isig integer(i_kind) :: loc,nn,m1,m,i,n,j,k real(r_kind) :: corq2x real(r_kind) :: factoz allocate ( clat_avn(mlat) ) allocate ( sigma_avn(1:msig) ) allocate ( agv_avn(0:mlat+ione,1:msig,1:msig) ) allocate ( bv_avn(0:mlat+ione,1:msig),wgv_avn(0:mlat+ione,1:msig) ) ! 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 ! 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 ! compute vertical(pressure) interpolation index and weight do k=1,nsig rlsig(k)=log(ges_prslavg(k)/ges_psfcavg) enddo if(mype==izero) then write(6,*) myname_,'(PREWGT_REG): read error amplitudes ', & '"',trim(berror_stats),'". ', & 'mype,nsigstat,nlatstat =', & mype,nsigstat,nlatstat end if ! Read amplitudes nrf_err=.false. read: do read(inerr,iostat=istat) var, isig if (istat /= 0) exit allocate ( corz_avn(1:mlat,1:isig) ) allocate ( hwll_avn(0:mlat+1,1:isig) ) allocate ( vztdq_avn(1:isig,0:mlat+1) ) if (var/='q') then read(inerr) corz_avn else allocate ( corqq_avn(1:mlat,1:isig) ) read(inerr) corz_avn,corqq_avn end if read(inerr) hwll_avn if (isig>1) then read(inerr) vztdq_avn end if ! load the variances do n=1,nrf if (var==nrf_var(n)) then nrf_err(n)=.true. loc=n exit end if end do if (isig==msig) then do n=1,nrf3 if (nrf3_loc(n)==loc) then if (var=='q' .and. qoption==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 do i=0,mlat+ione hwll_tmp(i,k,n)=hwll_avn(i,k) vz_tmp(k,i,n)=vztdq_avn(k,i) end do end do exit end if end do end if if (isig==1) then do n=1,nrf2 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 hwllp(0,n)=hwll_avn(0,1) hwllp(mlat+1,n)=hwll_avn(mlat+1,1) exit end if end do end if deallocate ( corz_avn ) deallocate ( hwll_avn ) deallocate ( vztdq_avn ) if (var=='q') deallocate ( corqq_avn ) enddo read close(inerr) deallocate(clat_avn,sigma_avn) deallocate(agv_avn,bv_avn,wgv_avn) ! 3d variable do n=1,nrf3 loc=nrf3_loc(n) 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+ione 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 end if enddo 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 do n=1,nrf3 loc=nrf3_loc(n) if (nrf_err(loc)) cycle if (n==nrf3_oz) then factoz = 0.0002_r_kind*r25 corz(:,:,nrf3_oz)=factoz ! hwll(:,:,nrf3_oz)=400000.0_r_kind vz(:,:,nrf3_oz)=vz_oz nrf_err(loc)=.true. else if (n==nrf3_cw) then corz(:,:,nrf3_cw)=corz(:,:,nrf3_q) hwll(:,:,nrf3_cw)=hwll(:,:,nrf3_q) vz(:,:,nrf3_cw)=vz(:,:,nrf3_q) nrf_err(loc)=.true. end if end do ! 2d variable do n=1,nrf2 loc=nrf2_loc(n) if (nrf_err(loc)) cycle if (n==nrf2_sst) 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) hwllp(i,nrf2+1)=hwll(i,1,nrf3_sf) hwllp(i,nrf2+2)=hwll(i,1,nrf3_sf) end do nrf_err(loc)=.true. end if enddo ! Do final check to make sure that background errors have been loaded for all variables if(mype==0) then do n=1,nrf if (.not. nrf_err(n)) then write(6,*) 'READ_WGT_REG: ***ERROR*** fail to load error variance for ', nrf_var(n) call stop2(333) end if end do end if return end subroutine berror_read_wgt_reg end module m_berror_stats_reg