subroutine penal(xhat) !$$$ subprogram documentation block ! . . . . ! subprogram: penal oberror tuning ! prgmmr: wu org: np23 date: 2005-08-26 ! ! abstract: randomized estimation of Tr(KH) and Tr(HK) and ! adaptive tuning ! ! ! program history log: ! 2005-08-15 wu - oberror tuning ! 2008-03-24 wu - use convinfo ikx as index for oberr tune ! 2008-05-27 safford - rm unused vars ! 2008-12-03 todling - update in light of state vector and obs binning ! ! usage: intt(st,rt) ! input argument list: ! xhat - increment in grid space ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_single,r_kind,i_kind use mpimod, only: ierror,mpi_comm_world,mpi_sum,mpi_rtype,mype use constants, only: izero,ione,zero,one use gsi_4dvar, only: nobs_bins use obsmod, only: qhead,qptr,thead,tptr,whead,wptr,pshead,psptr use converr, only:etabl use jfunc, only: jiterstart,jiter use convinfo, only:ictype,nconvtype,ioctype use state_vectors implicit none ! Declare passed variables type(state_vector),intent(in ) :: xhat ! Declare passed variables real(r_kind),save,dimension(33,200) :: penalty,trace ! Declare local variables real(r_kind) err2 integer(i_kind) i,n,k,l,m,ibin real(r_kind) tpenalty(33,nconvtype),ttrace(33,nconvtype) real(r_kind) valu,valv,val,so(33,nconvtype),cat_num(33,nconvtype),sosum,tcat_num(33,nconvtype) integer(i_kind) itype,ncat,k1 ncat=nconvtype*33 if(jiter==jiterstart)then trace=zero penalty=zero do ibin=1,nobs_bins ! Moisture qptr => qhead(ibin)%head do while (associated(qptr)) n=qptr%kx itype=ictype(n) if(itype==120_i_kind)then k1=qptr%k1 else k1=ione endif err2=qptr%raterr2*qptr%err2 ! err=sqrt(err2) ! Forward model val= qptr%wij(1)* xhat%q(qptr%ij(1))+qptr%wij(2)* xhat%q(qptr%ij(2))& +qptr%wij(3)* xhat%q(qptr%ij(3))+qptr%wij(4)* xhat%q(qptr%ij(4))& +qptr%wij(5)* xhat%q(qptr%ij(5))+qptr%wij(6)* xhat%q(qptr%ij(6))& +qptr%wij(7)* xhat%q(qptr%ij(7))+qptr%wij(8)* xhat%q(qptr%ij(8)) trace(k1,n)=trace(k1,n)-qptr%qpertb*val*err2 penalty(k1,n)=penalty(k1,n)+(val-qptr%res)**2*err2 qptr => qptr%llpoint end do ! if(mype==29_i_kind)write(0,*)'q2 trace,pen=',trace(k1,n),penalty(k1,n),k1,n ! Temperature tptr => thead(ibin)%head do while (associated(tptr)) n=tptr%kx itype=ictype(n) if(itype==120_i_kind)then k1=tptr%k1 else k1=ione endif err2=tptr%raterr2*tptr%err2 ! err=sqrt(err2) ! Forward model val= tptr%wij(1)* xhat%t(tptr%ij(1))+tptr%wij(2)* xhat%t(tptr%ij(2))& +tptr%wij(3)* xhat%t(tptr%ij(3))+tptr%wij(4)* xhat%t(tptr%ij(4))& +tptr%wij(5)* xhat%t(tptr%ij(5))+tptr%wij(6)* xhat%t(tptr%ij(6))& +tptr%wij(7)* xhat%t(tptr%ij(7))+tptr%wij(8)* xhat%t(tptr%ij(8)) trace(k1,n)=trace(k1,n)-tptr%tpertb*val*err2 penalty(k1,n)=penalty(k1,n)+(val-tptr%res)**2*err2 tptr => tptr%llpoint end do ! Surface pressure psptr => pshead(ibin)%head do while (associated(psptr)) n=psptr%kx itype=ictype(n) k1=ione err2=psptr%raterr2*psptr%err2 ! err=sqrt(err2) ! Forward model val= psptr%wij(1)* xhat%p(psptr%ij(1))+psptr%wij(2)* xhat%p(psptr%ij(2))& +psptr%wij(3)* xhat%p(psptr%ij(3))+psptr%wij(4)* xhat%p(psptr%ij(4)) trace(k1,n)=trace(k1,n)-psptr%ppertb*val*err2 penalty(k1,n)=penalty(k1,n)+(val-psptr%res)**2*err2 psptr => psptr%llpoint end do ! Winds wptr => whead(ibin)%head do while (associated(wptr)) n=wptr%kx itype=ictype(n) if(itype==220_i_kind .or. itype==223_i_kind .or. itype==233_i_kind .or. itype==245_i_kind)then k1=wptr%k1 else k1=ione endif err2=wptr%raterr2*wptr%err2 ! err=sqrt(err2) ! Forward model valu= wptr%wij(1)* xhat%u(wptr%ij(1))+wptr%wij(2)* xhat%u(wptr%ij(2))& +wptr%wij(3)* xhat%u(wptr%ij(3))+wptr%wij(4)* xhat%u(wptr%ij(4))& +wptr%wij(5)* xhat%u(wptr%ij(5))+wptr%wij(6)* xhat%u(wptr%ij(6))& +wptr%wij(7)* xhat%u(wptr%ij(7))+wptr%wij(8)* xhat%u(wptr%ij(8)) valv= wptr%wij(1)* xhat%v(wptr%ij(1))+wptr%wij(2)* xhat%v(wptr%ij(2))& +wptr%wij(3)* xhat%v(wptr%ij(3))+wptr%wij(4)* xhat%v(wptr%ij(4))& +wptr%wij(5)* xhat%v(wptr%ij(5))+wptr%wij(6)* xhat%v(wptr%ij(6))& +wptr%wij(7)* xhat%v(wptr%ij(7))+wptr%wij(8)* xhat%v(wptr%ij(8)) trace(k1,n)=trace(k1,n)-(wptr%upertb*valu+wptr%vpertb*valv)*err2 penalty(k1,n)=penalty(k1,n)+((valu-wptr%ures)**2+(valv-wptr%vres)**2)*err2 wptr => wptr%llpoint end do end do ! ibin else ! jiter cat_num=zero do ibin=1,nobs_bins ! Moisture ! ratiomin=one qptr => qhead(ibin)%head do while (associated(qptr)) n=qptr%kx itype=ictype(n) if(itype==120_i_kind)then k1=qptr%k1 else k1=ione endif err2=qptr%raterr2*qptr%err2 ! err=sqrt(err2) ! Forward model val= qptr%wij(1)* xhat%q(qptr%ij(1))+qptr%wij(2)* xhat%q(qptr%ij(2))& +qptr%wij(3)* xhat%q(qptr%ij(3))+qptr%wij(4)* xhat%q(qptr%ij(4))& +qptr%wij(5)* xhat%q(qptr%ij(5))+qptr%wij(6)* xhat%q(qptr%ij(6))& +qptr%wij(7)* xhat%q(qptr%ij(7))+qptr%wij(8)* xhat%q(qptr%ij(8)) cat_num(k1,n)=cat_num(k1,n)+one trace(k1,n)=trace(k1,n)+qptr%qpertb*val*err2 qptr => qptr%llpoint end do ! if(mype==29_i_kind)write(0,*)'q2 trace,pen=',trace(k1,n),cat_num(k1,n),k1,n ! Temperature tptr => thead(ibin)%head do while (associated(tptr)) n=tptr%kx itype=ictype(n) if(itype==120_i_kind)then k1=tptr%k1 else k1=ione endif err2=tptr%raterr2*tptr%err2 ! err=sqrt(err2) ! Forward model val= tptr%wij(1)* xhat%t(tptr%ij(1))+tptr%wij(2)* xhat%t(tptr%ij(2))& +tptr%wij(3)* xhat%t(tptr%ij(3))+tptr%wij(4)* xhat%t(tptr%ij(4))& +tptr%wij(5)* xhat%t(tptr%ij(5))+tptr%wij(6)* xhat%t(tptr%ij(6))& +tptr%wij(7)* xhat%t(tptr%ij(7))+tptr%wij(8)* xhat%t(tptr%ij(8)) cat_num(k1,n)=cat_num(k1,n)+one trace(k1,n)=trace(k1,n)+tptr%tpertb*val*err2 tptr => tptr%llpoint end do ! Surface pressure psptr => pshead(ibin)%head do while (associated(psptr)) n=psptr%kx itype=ictype(n) k1=ione err2=psptr%raterr2*psptr%err2 ! err=sqrt(err2) ! Forward model val= psptr%wij(1)* xhat%p(psptr%ij(1))+psptr%wij(2)* xhat%p(psptr%ij(2))& +psptr%wij(3)* xhat%p(psptr%ij(3))+psptr%wij(4)* xhat%p(psptr%ij(4)) cat_num(k1,n)=cat_num(k1,n)+one trace(k1,n)=trace(k1,n)+psptr%ppertb*val*err2 psptr => psptr%llpoint end do ! Winds wptr => whead(ibin)%head do while (associated(wptr)) n=wptr%kx itype=ictype(n) if(itype==220_i_kind .or. itype==223_i_kind .or. itype==233_i_kind .or. itype==245_i_kind)then k1=wptr%k1 else k1=ione endif err2=wptr%raterr2*wptr%err2 ! err=sqrt(err2) ! Forward model valu= wptr%wij(1)* xhat%u(wptr%ij(1))+wptr%wij(2)* xhat%u(wptr%ij(2))& +wptr%wij(3)* xhat%u(wptr%ij(3))+wptr%wij(4)* xhat%u(wptr%ij(4))& +wptr%wij(5)* xhat%u(wptr%ij(5))+wptr%wij(6)* xhat%u(wptr%ij(6))& +wptr%wij(7)* xhat%u(wptr%ij(7))+wptr%wij(8)* xhat%u(wptr%ij(8)) valv= wptr%wij(1)* xhat%v(wptr%ij(1))+wptr%wij(2)* xhat%v(wptr%ij(2))& +wptr%wij(3)* xhat%v(wptr%ij(3))+wptr%wij(4)* xhat%v(wptr%ij(4))& +wptr%wij(5)* xhat%v(wptr%ij(5))+wptr%wij(6)* xhat%v(wptr%ij(6))& +wptr%wij(7)* xhat%v(wptr%ij(7))+wptr%wij(8)* xhat%v(wptr%ij(8)) cat_num(k1,n)=cat_num(k1,n)+one trace(k1,n)=trace(k1,n)+(wptr%upertb*valu+wptr%vpertb*valv)*err2 wptr => wptr%llpoint end do do n=1,nconvtype do k=1,33 trace(k,n)=cat_num(k,n)-trace(k,n) enddo enddo end do ! ibin call mpi_allreduce(trace,ttrace,ncat,mpi_rtype,mpi_sum, & mpi_comm_world,ierror) call mpi_allreduce(penalty,tpenalty,ncat,mpi_rtype,mpi_sum, & mpi_comm_world,ierror) call mpi_allreduce(cat_num,tcat_num,ncat,mpi_rtype,mpi_sum, & mpi_comm_world,ierror) if(mype==izero)then do n=1,nconvtype write(233,*)'obs type=',ictype(n),trim(ioctype(n)) do k=1,33 if(tcat_num(k,n)>zero .and. tcat_num(k,n)<10._r_kind)write(223,*)k,n,tcat_num(k,n) write(233,*)k,n,tpenalty(k,n),ttrace(k,n),tcat_num(k,n) enddo enddo so=one do n=1,nconvtype do k=1,33 if(ttrace(k,n) /= zero .and. tcat_num(k,n)>10._r_kind) then so(k,n)=tpenalty(k,n)/ttrace(k,n) write(234,*)k,n,ictype(n),trim(ioctype(n)),so(k,n) endif if(so(k,n) >= zero) then so(k,n)=sqrt(so(k,n)) else so(k,n)=one endif enddo enddo sosum=zero do i=1,ncat sosum=sosum+(so(i,1)-one)**2 enddo write(235,*)'sosum=',sosum ! Update etabl do i=1,nconvtype if(trim(ioctype(i))=='t')then m=2_i_kind elseif(trim(ioctype(i))=='q')then m=3_i_kind elseif(trim(ioctype(i))=='uv')then m=4_i_kind elseif(trim(ioctype(i))=='ps')then m=5_i_kind else cycle endif l=ictype(i) ! Enough obs to define the vertical profile if((l==120_i_kind.and.m/=5_i_kind) .or. l==220_i_kind .or. l==223_i_kind .or. l==233_i_kind .or. l==245_i_kind)then write(235,*)l,trim(ioctype(i)),'33' do k=1,33 if( etabl(l,k,m) < 1.e8_r_single) etabl(l,k,m)=etabl(l,k,m)*so(k,i) end do else write(235,*)l,trim(ioctype(i)),'1' do k=1,33 if( etabl(l,k,m) < 1.e8_r_single) etabl(l,k,m)=etabl(l,k,m)*so(1,i) end do endif enddo ! Write out err table open(59,file='errtable_out',form='formatted') rewind 59 do l=100,299 if(etabl(l,1,1)==1100._r_single)then write(59,100)l 100 format(1x,i3,' OBSERVATION TYPE') do k=1,33 write(59,110)(etabl(l,k,i),i=1,6) 110 format(1x,6e12.5) end do endif ! etable1=1100 end do close(59) endif ! mype==0 call mpi_finalize(ierror) stop endif ! jiter return end subroutine penal