program innov_norml !********************************************************************** ! * . . . ! * PROGRAM: INNOV_NORML ! * PRGMMR: M. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-09 ! * ! * ABSTRACT: CALCULATE NORMALIZED INNOVATIONS (R+HPH**T)**(-1/2)*[y-Hx] ! * ! * PROGRAM LOG: ! * ! * 09/09/2003 ..... M. ZUPANSKI: Original 'transform.F' ! * 10/17/2003 ..... M. ZUPANSKI: Innovation normalization ! * ! ********************************************************************** !---------------------------------- ! INPUT: ! sqrt(Pf) - fcst_xxxx ! sqrt(I+A)**(-1/2) - sqrtIPA_inv ! z - control variable in ensemble space ! ! OUTPUT: ! x - 1-d control variable in model space ! !---------------------------------- #ifdef MPI_USE include "mpif.h" #endif !----- integer,parameter::ires_cntrl=20 ! res_0000 file # integer,parameter::ires=21 ! res_0000 file # integer,parameter::ieigen=22 ! A_eigenfile file # integer,parameter::iV_file=23 ! V eigenvector file # integer,parameter::incr=50 ! Obs increment file # integer,parameter::innov_file=51 ! Normalized innovation file # !----- real,parameter::Weps=1.e-04 ! threshold for eigenvalue (inverse calculation) !----- integer :: NENS integer :: N_obs integer :: N_LOC,jsta,jend integer :: NPROC,MPIRANK character*20 :: name_res character*20 :: name_eigen character*20 :: name_ens character*20 :: name_incr character*30 :: fileout character*30 :: FILNAMEJ character*30 :: FILE_INCR integer,allocatable :: jlen(:),jdisp(:) real,dimension(:),allocatable::sum_loc,sumx real,dimension(:),allocatable::x,x_b real,dimension(:),allocatable::W,W_h real,dimension(:,:),allocatable::Z,zz real,dimension(:),allocatable::res,res_cntrl,res_incr real,dimension(:),allocatable::V real,dimension(:),allocatable::z_1 real,dimension(:),allocatable::x_innov !----- ! ! DECLARE NAMELIST ! NAMELIST /ENSEMBLE_SIZE/ NENS !==============start calculation=================== #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 !================================================ REWIND 15 READ(15,ENSEMBLE_SIZE) if(MPIRANK.eq.0) then write(*,*) "innov_norml: number of ensemble members=",NENS end if 100 format(A,i4.4) !-- input filename write(name_res,101) 101 format('res.post_') !-- output filename write(name_eigen,102) 102 format('eigenV_') !-- ensemble matrix name write(name_ens,103) 103 format('A_eigenfile') !-- increment name write(name_incr,104) 104 format('incr_') !-- A_eigenfile -------------------- CLOSE(ieigen) OPEN(UNIT=ieigen,FILE=name_ens,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) ieigen,' OPEN UNIT ERROR IER=',IER READ(ieigen) NN if(NN.ne.NENS) then write(*,*) "Innov_norml DIMENSION PROBLEM: NN=",NN," NENS=",NENS stop end if allocate(W_h(1:NENS)) allocate(Z(1:NENS,1:NENS)) READ(ieigen) (W_h(i),i=1,NENS) READ(ieigen) ((Z(i,j),i=1,NENS),j=1,NENS) CLOSE(ieigen,STATUS='KEEP') write(*,*) "--------------------" write(*,*) "Eigenvalues (I+W)**(-1/2) and eigenvectors U " do i=1,8 write(*,111) i,W_h(i),(Z(i,j),j=1,8) end do 111 format(I3,E10.3,8(E10.3)) write(*,*) "--------------------" !== from (I+Lambda)**(-1/2)=W_h, calculate Lambda*(-1/2)=W ========= !== W=W_h/sqrt(1-W_h**2) allocate(W(1:NENS)) do N=1,NENS diff=1./W_h(N)**2-1. write(*,*) N," Lambda=",diff if(diff.ge.Weps) then W(N)=1./sqrt(diff) else W(N)=0.0 end if end do if(MPIRANK.eq.0) then do ii=1,NENS write(*,*) ii," Lambda**(-1/2)=",W(ii)," (I+Lamb)**(-1/2)=",W_h(ii) end do end if !-- create zz=Lambda**(-1/2)*U ------------- allocate(zz(1:NENS,1:NENS)) do j=1,NENS zz(i,j)=0.0 do i=1,NENS zz(i,j)=zz(i,j)+W(j)*Z(i,j) end do end do deallocate(W) ; deallocate(Z) !==================================================================== !== !== define local dimensions and indexes !== allocate(jdisp(0:NPROC-1)) ; allocate(jlen(0:NPROC-1)) do irank=0,NPROC-1 CALL para_range(1,NENS,NPROC,irank,jsta,jend) jlen(irank)=jend-jsta+1 jdisp(irank)=jsta-1 end do #ifdef MPI_USE CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #endif CALL para_range(1,NENS,NPROC,MPIRANK,jsta,jend) !== N_LOC=jlen(MPIRANK) !== ! write(*,*) "NENS,NPROC,MPIRANK,jsta,jend=",NENS,NPROC,MPIRANK,jsta,jend if(MPIRANK.eq.0) then do nn=0,NPROC-1 write(*,*) "nn=",nn," jdisp=",jdisp(nn)," jlen=",jlen(nn) end do end if !==== calculate and write-out Z-matrix (ensemble increments at obs location) === ! ! read in the control obs residual ! OPEN(UNIT=ires_cntrl,FILE='res.post_0000',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) ires_cntrl,' OPEN UNIT ERROR IER=',IER READ(ires_cntrl) N_obs allocate(res_cntrl(1:N_obs)) READ(ires_cntrl) res_cntrl CLOSE(ires_cntrl,STATUS='KEEP') allocate(res_incr(1:N_obs)) allocate(res(1:N_obs)) do jj=jsta,jend WRITE(FILNAMEJ,100) trim(name_res),jj OPEN(UNIT=ires,FILE=FILNAMEJ,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) ires,' OPEN UNIT ERROR IER=',IER READ(ires) NN if(NN.ne.N_obs) then write(*,*) "Innov_norml RES DIMENSION PROBLEM: NN=",NN," N_obs=",N_obs stop end if READ(ires) res CLOSE(ires,STATUS='KEEP') !-- increments do ii=1,N_obs res_incr(ii)=res_cntrl(ii)-res(ii) end do WRITE(FILE_INCR,100) trim(name_incr),jj OPEN(UNIT=incr,FILE=FILE_INCR,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) incr,' OPEN UNIT ERROR IER=',IER WRITE(incr) N_obs WRITE(incr) res_incr CLOSE(incr,STATUS='KEEP') end do deallocate(res_incr) ; deallocate(res) !====================== matrix-matrix product (get V eigenvectots) ================= call MATRIX_MATRIX & (incr,iV_file,name_incr,name_eigen,NENS,MPIRANK,NPROC,jsta,N_LOC,zz) deallocate(zz) !==================================================================== ! (V**T)*res_0000 allocate(sum_loc(1:N_LOC)) sum_loc(:)=0.0 ii=0 do N=jsta,jend ii=ii+1 write(fileout,100) trim(name_eigen),N CLOSE(iV_file) OPEN(UNIT=iV_file,FILE=fileout,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) iV_file,' OPEN UNIT ERROR IER=',IER READ(iV_file) NNx if(NNx.ne.N_obs) then write(*,*) "Innov_norml DIMENSION PROBLEM: NNx=",NNx," N_obs=",N_obs stop end if allocate(V(1:N_obs)) READ(iV_file) V CLOSE(iV_file,STATUS='KEEP') do i=1,N_obs sum_loc(ii)=sum_loc(ii)+res_cntrl(i)*V(i) end do end do allocate(sumx(1:NENS)) #ifdef MPI_USE CALL MPI_ALLGATHERV & (sum_loc,N_LOC,MPI_REAL4,sumx,jlen,jdisp,MPI_REAL4,MPI_COMM_WORLD,IERR) #else sumx=sum_loc #endif deallocate(sum_loc) deallocate(res_cntrl) !==================================================================== ! (I+Lambda)**(-1/2)*(V**T)*res allocate(z_1(1:NENS)) do i=1,NENS z_1(i)=sumx(i)*W_h(i) end do deallocate(sumx) ; deallocate(W_h) !==================================================================== allocate(sum_loc(1:N_obs)) allocate(sumx(1:N_obs)) call MATRIX_VECTOR & (iV_file,name_eigen,NENS,MPIRANK,NPROC,jsta,jend,N_obs,z_1,sum_loc) #ifdef MPI_USE CALL MPI_ALLREDUCE & (sum_loc,sumx,N_obs,MPI_REAL4,MPI_SUM,MPI_COMM_WORLD,IERR) #else sumx=sum_loc #endif deallocate(sum_loc) deallocate(z_1) deallocate(jlen) ; deallocate(jdisp) !========== write 1-d control variable ============== !---- read background (first guess) vector if(MPIRANK.eq.0) then allocate(x_innov(1:N_obs)) x_innov(:)=sumx(:) CLOSE(innov_file) OPEN(UNIT=innov_file,FILE='innov',FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(*,*) ' imodelb OPEN UNIT ERROR IER=',IER write(innov_file) N_obs write(innov_file) x_innov CLOSE(innov_file,status='KEEP') deallocate(x_innov) end if deallocate(sumx) if(MPIRANK.eq.0) then write(*,*) "end innov_norml" end if #ifdef MPI_USE CALL MPI_FINALIZE(IERR) #endif stop end