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: ! U,Lambda - EVD of matrix A=(Z**T)*Z [ (I+lambda)**(-1/2) = W_h ] ! res_cntrl- innovation vector: R**(-1/2)*[y-H(x)] ! Z - observation matrix: R**(-1/2)*[H(x+dx)-H(x)] ! ! OUTPUT: ! x_innov - normalized innovation vector: [I+Z*U*sigma*(U**T)*(Z**T)]*res_cntrl !---------------------------------- #ifdef MPI_USE include "mpif.h" #endif !----- integer,parameter::ires_cntrl=20 ! res_cntrl file # integer,parameter::ires=21 ! res_iiii file # integer,parameter::ieigen=22 ! A_eigenfile file # integer,parameter::incr=50 ! Obs increment file # integer,parameter::innov_file=51 ! Normalized innovation file # !----- integer :: NENS integer :: NENS_START integer :: NALL integer :: N_obs integer :: N_LOC,jsta,jend integer :: NPROC,MPIRANK character*20 :: name_res 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::res,res_cntrl,res_incr real,dimension(:),allocatable::z_1,z_2,z_3 real,dimension(:),allocatable::x_innov real,dimension(:),allocatable::xlambda real,dimension(:),allocatable::beta0 real,dimension(:),allocatable::gammak_old,gammak real,dimension(:),allocatable::sigmak_old,sigmak real,dimension(:),allocatable::W_h real,dimension(:,:),allocatable::U !----- ! ! DECLARE NAMELIST ! NAMELIST /ENSEMBLE_SIZE/ NENS_START,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: NENS_START, NENS=",NENS_START,NENS end if !- Total number of degrees of freedom NALL=NENS-NENS_START+1 100 format(A,i4.4) !-- ensemble matrix name write(name_ens,101) 101 format('A_eigenfile') !-- input filename write(name_res,102) 102 format('res_') !-- normalized innovation file name write(name_incr,103) 103 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.NALL) then write(*,*) "Innov_norml DIMENSION PROBLEM: NN=",NN," NALL=",NALL stop end if allocate(W_h(NENS_START:NENS)) allocate(U(NENS_START:NENS,NENS_START:NENS)) READ(ieigen) (W_h(i),i=NENS_START,NENS) READ(ieigen) ((U(i,j),i=NENS_START,NENS),j=NENS_START,NENS) CLOSE(ieigen,STATUS='KEEP') write(*,*) "--------------------" write(*,*) "Eigenvalues (I+W)**(-1/2) and eigenvectors U " do i=NENS_START,8 write(*,111) i,W_h(i),(U(i,j),j=NENS_START,8) end do 111 format(I3,E10.3,8(E10.3)) write(*,*) "--------------------" !== from (I+Lambda)**(-1/2)=W_h, calculate Lambda =================== allocate(xlambda(NENS_START:NENS)) do N=NENS_START,NENS diff=1./W_h(N)**2-1. xlambda(N)=max(diff,0.0) end do allocate(beta0(NENS_START:NENS)) allocate(sigmak_old(NENS_START:NENS)) allocate(gammak_old(NENS_START:NENS)) allocate(sigmak(NENS_START:NENS)) allocate(gammak(NENS_START:NENS)) sigmak_old(:)=0.0 gammak_old(:)=0.0 sigmak(:)=0.0 gammak(:)=0.0 beta0(:)=-1.0*W_h(:)**2 Nrec=10 do krec=1,Nrec do i=NENS_START,NENS sigmak(i)=0.5*(sigmak_old(i)+beta0(i)-gammak_old(i)-beta0(i)*xlambda(i)*gammak_old(i)) xinv=1./(1.+xlambda(i)*sigmak(i)) gammak(i)=sigmak(i)*xinv sigmak_old(i)=sigmak(i) gammak_old(i)=gammak(i) end do if(MPIRANK.eq.0) then do ii=NENS_START,NENS write(*,200) krec,xlambda(ii),sigmak(ii),gammak(ii) end do end if end do 200 format(I3,' xlambda=',E12.5,' sigma=',E12.5,' gamma=',E12.5) deallocate(xlambda) deallocate(gammak_old) ; deallocate(gammak) deallocate(sigmak_old) deallocate(W_h) !==================================================================== !== !== 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) !== ! write(*,*) "NALL,NPROC,MPIRANK,jsta,jend=",NALL,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_cntrl',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)) !=================================================================== ! Create Z=R**(-1/2)*(H(x+dx)-H(x)) ! (Z**T)*res_cntrl allocate(sum_loc(1:N_LOC)) sum_loc(:)=0.0 ii=0 do jj=jsta,jend ii=ii+1 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 (columns of Z-matrix) do kk=1,N_obs res_incr(kk)=res_cntrl(kk)-res(kk) end do !-- save for later 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') !-- (Z**T)*res_cntrl do i=1,N_obs sum_loc(ii)=sum_loc(ii)+res_cntrl(i)*res_incr(i) end do end do deallocate(res_incr) ; deallocate(res) !--------------------------------------------------- allocate(sumx(1:NALL)) #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) !==================================================================== ! U*(Z**T)*res allocate(z_1(NENS_START:NENS)) do i=NENS_START,NENS z_1(i)=0.0 Kindex=0 do j=NENS_START,NENS Kindex=Kindex+1 z_1(i)=z_1(i)+U(j,i)*sumx(Kindex) end do end do deallocate(sumx) !==================================================================== ! Sigma*(U**T)*(Z**T)*res allocate(z_2(NENS_START:NENS)) do i=NENS_START,NENS z_2(i)=sigmak(i)*z_1(i) end do deallocate(z_1) ; deallocate(sigmak) !==================================================================== ! U*Sigma*(U**T)*(Z**T)*res allocate(z_3(NENS_START:NENS)) do i=NENS_START,NENS z_3(i)=0.0 do j=NENS_START,NENS z_3(i)=z_3(i)+U(i,j)*z_2(j) end do end do deallocate(z_2) ; deallocate(U) !==================================================================== ! Z*U*Sigma*(U**T)*(Z**T)*res allocate(sum_loc(1:N_obs)) allocate(sumx(1:N_obs)) call MATRIX_VECTOR & (incr,name_incr,NENS_START,NENS,MPIRANK,NPROC,jsta,jend,N_obs,z_3,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_3) deallocate(jlen) ; deallocate(jdisp) !---- read background (first guess) vector if(MPIRANK.eq.0) then allocate(x_innov(1:N_obs)) do i=1,N_obs x_innov(i)=res_cntrl(i)+sumx(i) end do 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(res_cntrl) ; deallocate(sumx) if(MPIRANK.eq.0) then write(*,*) "end innov_norml" end if #ifdef MPI_USE CALL MPI_FINALIZE(IERR) #endif stop end