program chi_square_bound ! !---------------------------------------------------------------------------------------------- !--- Scale the sqrt analysis error covariance if out of preassigned bounds for chi-square test -- ! If NOT, do not change sqrtPa ------------------------------------------------------------- !---------------------------------------------------------------------------------------------- ! ! INPUT: x_1 = normalized innovations ! x = control state (for dimensions only) ! ! INPUT,OUTPUT: sqrtPx = sqrt analysis error covariance ! !-------------------------------------------------------------------- #ifdef MPI_USE include "mpif.h" #endif real,parameter::chi_square_max=1.10 ! the upper bound for chi_square ! (may be adjusted by the ensemble size as well) integer,parameter::chi=20 ! current chi_square integer,parameter::icntrl=21 ! control x integer,parameter::ieigen=22 ! A_eigenfile file integer,parameter::isqrtPx=50 ! ensemble perturbation (sqrtPa) integer idim !----- real,dimension(:),allocatable :: W_h,xlambda real,dimension(:),allocatable :: sqrtPx integer,allocatable :: jlen(:),jdisp(:) real :: chi_square,chi_rel_err real :: eps real :: trace_chi,trace_Paobs integer :: N_LOC,jsta,jend integer :: NPROC,MPIRANK integer :: NENS integer :: NENS_START integer :: NALL character*11 Pxname NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS !==================================================================== #ifdef MPI_USE ! ! START MPI ! CALL MPI_INIT(IERR) IF (IERR .NE. MPI_SUCCESS) STOP 'FAILED TO INIT MPI' CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPROC, IERR) CALL MPI_COMM_RANK(MPI_COMM_WORLD, MPIRANK, IERR) CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #else write(*,*) "This is a NO MPI run" MPIRANK=0 NPROC=1 #endif if(MPIRANK.eq.0) then write(*,*) "start chi_square_bound" end if !-- read ensemble size REWIND 15 READ(15,ENSEMBLE_SIZE) if(MPIRANK.eq.0) then write(*,*) "chi_square_bound: NENS_START, NENS=",NENS_START,NENS end if ! Total number of degrees of freedom NALL=NENS-NENS_START+1 !==================================================================== !-- read chi_square --------------- rewind chi read(chi,500) N_dim,chi_square 500 format(I15,E20.10) write(*,*) "INPUT chi_square=",chi_square," N_obs=",N_dim !==============start calculation=================== !============ if(chi_square.gt.chi_square_max) then !============ chi_rel_err=chi_square_max/chi_square write(*,*) "sqrtPa inflation NEEDED: chi_rel_err=",chi_rel_err !-- A_eigenfile -------------------- CLOSE(ieigen) OPEN(UNIT=ieigen,FILE='A_eigenfile',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) ieigen,' OPEN UNIT ERROR IER=',IER READ(ieigen) NN if(NN.ne.NALL) then write(*,*) MPIRANK," chi_square_bound DIMENSION PROBLEM: NN=",NN," NALL=",NALL stop end if allocate(W_h(NENS_START:NENS)) READ(ieigen) (W_h(i),i=NENS_START,NENS) READ(ieigen) CLOSE(ieigen,STATUS='KEEP') !== from (I+Lambda)**(-1/2)=W_h, calculate Lambda =================== allocate(xlambda(NENS_START:NENS)) do N=NENS_START,NENS W_h(N)=W_h(N)*sqrt(float(N_dim)) diff=1./W_h(N)**2-1. xlambda(N)=max(diff,0.0) end do !--- calculate trace ------------------- ! trace_chi = sum [1./(1+xlambda)] ! trace_Paobs = sum [xlambda/(1+xlambda)] trace_chi=0.0 trace_Paobs=0.0 DO k=1,NALL trace_chi=trace_chi+W_h(k)*W_h(k) trace_Paobs=trace_Paobs+xlambda(k)*W_h(k)*W_h(k) END DO deallocate(W_h) ; deallocate(xlambda) ! if(MPIRANK.eq.0) then ! write(*,*) "------------ " ! write(*,*) " trace_chi=",trace_chi," trace_Paobs=",trace_Paobs ! write(*,*) "------------ " ! end if !! eps=sqrt((NALL-chi_rel_err*trace_chi)/trace_Paobs) CLOSE(54) OPEN(54,file='trace',status='unknown',form='formatted') read(54,*) Rodgers_inf_cont,trace_inv CLOSE(54) write(*,*) " trace_inv=",trace_inv," Rodgers_inf_cont=",Rodgers_inf_cont xx=(N_dim-chi_rel_err*trace_inv)/Rodgers_inf_cont if(xx.gt.0.0) then eps=sqrt(xx) else eps=1.0 end if ! eps=sqrt((N_dim-chi_rel_err*trace_chi)/trace_Paobs) if(eps.gt.1.) then if(MPIRANK.eq.0) then write(*,*) "chi_square_bound is proceeding NORMALLY" write(*,*) "chi_rel_err=",chi_rel_err write(*,*) "eps=",eps end if !==================================================================== !== !== define local dimensions and indexes !== allocate(jdisp(0:NPROC-1)) ; allocate(jlen(0:NPROC-1)) do irank=0,NPROC-1 CALL para_range(NENS_START,NENS,NPROC,irank,jsta,jend) jlen(irank)=jend-jsta+1 jdisp(irank)=jsta-NENS_START end do #ifdef MPI_USE CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #endif CALL para_range(NENS_START,NENS,NPROC,MPIRANK,jsta,jend) !== N_LOC=jlen(MPIRANK) !== ! if(MPIRANK.eq.0) then do nn=0,NPROC-1 write(*,*) "nn=",nn," jdisp=",jdisp(nn)," jlen=",jlen(nn) end do end if !==================================================================== !-- read control variable, just to get dimensions CLOSE(icntrl) OPEN(icntrl,file='x1d_cntrl',status='unknown',form='unformatted') READ(icntrl) idim write(*,*) "input: idim=",idim READ(icntrl) CLOSE(icntrl,status='KEEP') !==================================================================== allocate(sqrtPx(1:idim)) 2000 format('sqrtPa_',i4.4) DO N=jsta,jend WRITE(Pxname,2000) N CLOSE(isqrtPx) OPEN(isqrtPx,file=Pxname,status='unknown',form='unformatted') READ(isqrtPx) READ(isqrtPx) sqrtPx !--- input statistics -------------------- sumP=0.0 do i=1,idim sumP=sumP+sqrtPx(i)**2 end do sumP=sqrt(sumP/float(idim)) write(*,100) Pxname,minval(sqrtPx),maxval(sqrtPx),sumP do i=1,idim sqrtPx(i) = eps * sqrtPx(i) end do !-- write perturbation (sqrtPx) BACKSPACE isqrtPx WRITE(isqrtPx) sqrtPx CLOSE(isqrtPx,status='KEEP') !--- output statistics -------------------- sumP=0.0 do i=1,idim sumP=sumP+sqrtPx(i)**2 end do sumP=sqrt(sumP/float(idim)) write(*,200) Pxname,minval(sqrtPx),maxval(sqrtPx),sumP !----------------------------------- END DO !-------------------------------------------- 100 format('INPUT: ',a11,' min,max,avg=',3E15.5) 200 format('OUTPUT: ',a11,' min,max,avg=',3E15.5) deallocate(sqrtPx) else if(MPIRANK.eq.0) then write(*,*) "PROBLEM: Although inflation is anticipated, the values are NOT GOOD!" write(*,*) "chi_rel_err=",chi_rel_err," eps (should be greater than 1)=",eps write(*,*) "ERROR EXIT" end if end if !============ ELSE !============ chi_rel_err=1. write(*,*) "sqrtPa inflation NOT needed, NORMAL EXIT" !============ END IF !============ if(MPIRANK.eq.0) then write(*,*) "end chi_square_bound" end if !==================================================================== #ifdef MPI_USE CALL MPI_FINALIZE(IERR) #endif stop end program chi_square_bound