program A_matrix_elements !********************************************************************** ! * . . . ! * PROGRAM: A_matrix_elements ! * PRGMMR: M. ZUPANSKI ORG: CIRA/CSU DATE: 2003-08-28 ! * ! * ABSTRACT: SQUARE-ROOT ANALYSIS ERROR COVARIANCE CALCULATION ! * Calculate elements of the matrix A, as ensemble inner-products ! * in observation space. ! * ! * PROGRAM LOG: ! * ! * 08/28/2003 ..... M. ZUPANSKI: ! * ! ********************************************************************** ! -1/2 ! INPUT: x(i) = R [KM(x+delx)-KM(x)] ! -1/2 ! OUTPUT: a(i) = A - matrix (in upper-packed storage mode for EVD) ! !-------------------------------------------------------------------- #ifdef MPI_USE include "mpif.h" #endif integer,parameter::idim=21 ! res cntrl file integer,parameter::iptrb=22 ! res (ensemble ptrb at obs location) file integer,parameter::matrixa=51 ! A matrix file integer,parameter::matrixa_2d=52 ! A matrix file !----- !--------- non-dimensional obs perturbations (residuals) -------- real,allocatable::res_i(:) real,allocatable::res_j(:) real,allocatable::res_cntrl(:) !--------- inner-product A matrix (temporary) --------- real,allocatable::sumx(:) real,allocatable::sum_loc(:) ! !--------- symmetric matrix A elements (in upper-packed storage mode) ----------- real,allocatable::a(:) real,allocatable::a2d(:,:) integer,allocatable::jlen(:),jdisp(:) integer IERR integer NALL integer NENS integer NENS_NSTART CHARACTER*8 FILNAMEI,FILNAMEJ CHARACTER*4 ENVIR ! ! 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 !================================================ !------------------------------- if(MPIRANK.eq.0) then write(*,*) "A_matrix: start reading input 15" end if REWIND 15 READ(15,ENSEMBLE_SIZE) if(MPIRANK.eq.0) then write(*,*) "A_matrix: NENS_START, NENS=",NENS_START,NENS end if 1110 FORMAT('res_',I4.4) !==================================================================== !=== create A matrix == !==================================================================== NALL=NENS-NENS_START+1 NSIZE=NALL*(NALL+1)/2 write(*,*) "NSIZE=",NSIZE !== !== define global indexes for the local sub-domain (jsta,jend) !== allocate(jdisp(0:NPROC-1)) allocate(jlen(0:NPROC-1)) do irank=0,NPROC-1 CALL para_range(1,NSIZE,NPROC,irank,jsta,jend) jlen(irank)=jend-jsta+1 jdisp(irank)=jsta-1 end do CALL para_range(1,NSIZE,NPROC,MPIRANK,jsta,jend) #ifdef MPI_USE CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #endif !== N_LOC=jlen(MPIRANK) !== ! ! define ii,jj indexes (start and end of loop for each task) ! K=0 do JJ=NENS_START,NENS do II=NENS_START,JJ K=K+1 if(K.eq.jdisp(MPIRANK)+1) then islocal=ii jslocal=jj end if if(K.eq.jdisp(MPIRANK)+N_LOC) then ielocal=ii jelocal=jj end if end do end do if(MPIRANK.eq.0) then do nn=0,NPROC-1 write(*,*) "nn=",nn," jdisp=",jdisp(nn)," jlen=",jlen(nn) end do end if !==================================================================== !! MPI is used for inner-products ! !!-- a = total A matrix (sum of all obs times) !!-- sumx = partial A matrix (for given obs time) if(MPIRANK.eq.0) then allocate(a(1:NSIZE)) a=0.0 end if allocate(sumx(1:NSIZE)) ! ! read in the control obs residual dimensions ! CLOSE(idim) OPEN(idim,file='res_cntrl',status='unknown',form='unformatted') read(idim) mtot allocate(res_cntrl(1:max(1,mtot))) read(idim) (res_cntrl(k),k=1,mtot) CLOSE(idim,status='KEEP') if(MPIRANK.eq.0) then write(*,*) "res_cntrl min,max=",minval(res_cntrl),maxval(res_cntrl) write(*,*) "Obs residual dimension mtot=",mtot end if allocate(res_i(1:max(1,mtot))) allocate(res_j(1:max(1,mtot))) allocate(sum_loc(1:N_LOC)) res_i(:)=0.0 res_j(:)=0.0 #ifdef MPI_USE CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #endif kk=0 do JJ=jslocal,jelocal ! ! read in the ensemble perturbations ! WRITE(FILNAMEJ,1110) JJ CLOSE(iptrb) OPEN(iptrb,file=FILNAMEJ,status='unknown',form='unformatted') read(iptrb) read(iptrb) (res_j(k),k=1,mtot) CLOSE(iptrb,STATUS='KEEP') !! write(*,*) jj," res_j min,max=",minval(res_j),maxval(res_j) if(JJ.eq.jslocal) then ii_start=islocal else ii_start=NENS_START end if if(JJ.eq.jelocal) then ii_end=ielocal else ii_end=JJ end if do II=ii_start,ii_end kk=kk+1 WRITE(FILNAMEI,1110) II CLOSE(iptrb) OPEN(iptrb,file=FILNAMEI,status='unknown',form='unformatted') read(iptrb) read(iptrb) (res_i(k),k=1,mtot) CLOSE(iptrb,STATUS='KEEP') !! write(*,*) ii," res_i min,max=",minval(res_i),maxval(res_i) !---- inner product matrix -------- !-NOTE: R**(-1/2)*[H(x+dx)-H(x)] = res_cntrl - res sum_loc(kk)=0.0 do k=1,mtot sum_loc(kk)=sum_loc(kk)+(res_cntrl(k)-res_i(k))*(res_cntrl(k)-res_j(k)) end do end do end do deallocate(res_i) deallocate(res_j) deallocate(res_cntrl) ! !! write(*,*) MPIRANK," N_LOC=",N_LOC," kk=",kk ! sumx=0.0 ! #ifdef MPI_USE CALL MPI_GATHERV & (sum_loc,N_LOC,MPI_REAL4,sumx,jlen,jdisp,MPI_REAL4,0,MPI_COMM_WORLD,IERR) #else sumx=sum_loc #endif ! deallocate(sum_loc) !! END MPI !==================================================================== ! CREATE MATRIX A IN UPPER-PACKED STORAGE MODE !==================================================================== if(MPIRANK.eq.0) then ! Kindex=0 do jj=NENS_START,NENS do ii=NENS_START,jj Kindex=Kindex+1 a(Kindex)=a(Kindex)+sumx(Kindex) end do end do ! end if !==================================================================== if(MPIRANK.eq.0) then allocate(a2d(NENS_START:NENS,NENS_START:NENS)) a2d(:,:)=0.0 Kindex=0 do jj=NENS_START,NENS do ii=NENS_START,jj Kindex=Kindex+1 if(Kindex.lt.30) write(*,101) Kindex,ii,jj,a(Kindex) a2d(ii,jj)=a(Kindex) a2d(jj,ii)=a2d(ii,jj) end do end do write(*,*) "========" write(*,*) "A-MATRIX" write(*,*) "========" do ii=NENS_START,min(NENS,10) write(*,102) ii,(a2d(ii,jj),jj=NENS_START,min(8+NENS_START-1,NENS)) end do CLOSE(matrixa_2d) OPEN(matrixa_2d,file='matrix_A_2d',status='unknown', & form='unformatted') write(matrixa_2d) a2d CLOSE(matrixa_2d,status='KEEP') ! write out the A-matrix elements ! (all in upper-packed storage mode) CLOSE(matrixa) OPEN(matrixa,file='matrix_A',status='unknown',form='unformatted') write(matrixa) (a(K),K=1,NSIZE) CLOSE(matrixa,status='KEEP') deallocate(a2d) end if 101 format(I4,' A(',I3,',',I3,')=',E17.5) 102 format(I5,8E10.3) !==================================================================== #ifdef MPI_USE CALL MPI_FINALIZE(IERR) #endif stop end