program grad_obs ! !********************************************************************** ! * . . . ! * PROGRAM: grad_obs ! * PRGMMR: M. ZUPANSKI ORG: CIRA/CSU DATE: 2003-09-03 ! * ! * ABSTRACT: GRADIENT IN ENSEMBLE SPACE ! * Calculate observation cost function gradient with respect to cntrl vrbl ! * in ensemble subspace. ! * --- Calculate inner product for current obs time ! * --- Sum all to create the total gradient ! * ! * PROGRAM LOG: ! * ! * 09/03/2003 ..... M. ZUPANSKI: ! * ! ********************************************************************** ! -1/2 ! INPUT: x(i) = R [KM(x+delx)-KM(x)] ! -1/2 ! INPUT: z(i) = R [KM(x)-y_obs] ! T ! OUTPUT: g(i) = x z gradient in ensemble space ! !-------------------------------------------------------------------- #ifdef MPI_USE include "mpif.h" #endif !----- integer,parameter::idim0=20 ! res cntrl file integer,parameter::idim=21 ! res file integer,parameter::iptrb=22 ! res (all ensemble incr at obs location) integer,parameter::gradout=51 ! cost-function gradient !----- !--------- non-dimensional obs increments ----------------------- real,allocatable::res(:) real,allocatable::res_j(:) real,allocatable::res_cntrl(:) !--------- inner-products (temporary) --------- real,allocatable::sumx(:) real,allocatable::sum_loc(:) !--------- gradient ---------------------- real,allocatable::g(:) integer,allocatable::jlen(:),jdisp(:) integer IERR integer NENS integer NENS_START integer NALL integer mobs 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 !================================================ !-- read ensemble size REWIND 15 READ(15,ENSEMBLE_SIZE) if(MPIRANK.eq.0) then write(*,*) "grad_obs: NENS_START, NENS=",NENS_START,NENS end if 1110 FORMAT('res_',I4.4) !- Total number of degrees of freedom NALL=NENS-NENS_START+1 !== !== 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 CALL para_range(NENS_START,NENS,NPROC,MPIRANK,jsta,jend) #ifdef MPI_USE CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #endif !== N_LOC=jlen(MPIRANK) if(MPIRANK.eq.0) then do NN=0,NPROC-1 write(*,*) MPIRANK," jlen=",jlen(NN)," jdisp=",jdisp(NN) end do end if !==================================================================== !! MPI is used for inner-products ! !!-- g = total gradient in ensemble space (sum over all obs times) !!-- sumx = partial gradient (at given obs time) if(MPIRANK.eq.0) then allocate(g(NENS_START:NENS)) g=0.0 end if allocate(sumx(1:NALL)) ! ! read in the obs incr dimensions for each obs time ! and (control) innovation vectors ! CLOSE(idim0) OPEN(idim0,file='res_cntrl',status='unknown',form='unformatted') read(idim0) mobs allocate(res_cntrl(max(1,mobs))) read(idim0) (res_cntrl(k),k=1,mobs) CLOSE(idim0,STATUS='KEEP') if(MPIRANK.eq.0) then write(*,*) MPIRANK," res_cntrl: Obs incr dimension mobs=",mobs write(*,*) " min,max res_cntrl=",minval(res_cntrl),maxval(res_cntrl) end if ! ! read in the obs incr dimensions of importance ! allocate(res(max(1,mobs))) CLOSE(idim) OPEN(idim,file='res',status='unknown',form='unformatted') read(idim) read(idim) (res(k),k=1,mobs) CLOSE(idim,STATUS='KEEP') !------------------------------------------------------ allocate(res_j(1:max(1,mobs))) allocate(sum_loc(1:N_LOC)) res_j(:)=0.0 #ifdef MPI_USE CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) #endif nzero=0 nnonzero=0 kk=0 do JJ=jsta,jend kk=kk+1 ! ! read in the ensemble perturbations ! write(*,*) MPIRANK," read res_j: mobs=",mobs," JJ=",JJ WRITE(FILNAMEJ,1110) JJ CLOSE(iptrb) OPEN(iptrb,file=FILNAMEJ,status='unknown',form='unformatted') read(iptrb) read(iptrb) (res_j(k),k=1,mobs) CLOSE(iptrb,STATUS='KEEP') !---- inner product matrix -------- sum_loc(kk)=0.0 mnonzero=0 mnonzero_ctl=0 do k=1,mobs xx=(res_cntrl(k)-res_j(k))*res(k) if(res_cntrl(k).ne.0.0) mnonzero_ctl=mnonzero_ctl+1 sum_loc(kk)=sum_loc(kk)+(res_cntrl(k)-res_j(k))*res(k) if(xx.eq.0.0) then nzero=nzero+1 else nnonzero=nnonzero+1 mnonzero=mnonzero+1 !!! write(*,*)" kk=",kk," k=",k," res_cntrl,res_j=",res_cntrl(k),res_j(k) end if end do write(*,*) " kk=",kk," mnonzero=",mnonzero," mnonzero/mobs=",float(mnonzero)/mobs end do nntotal=nzero+nnonzero write(*,*) " kk=",kk," nzero=",nzero," nnonzero=",nnonzero," nntotal=",nntotal write(*,*) " mnonzero_ctl=",mnonzero_ctl,float(mnonzero_ctl)/mobs if(nntotal.gt.0) then write(*,*) " nzero=",float(nzero)/nntotal," nnonzero=",float(nnonzero)/nntotal endif deallocate(res_j) deallocate(res_cntrl) deallocate(res) ! 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 !==================================================================== ! GRADIENT g IN ENSEMBLE SPACE ! {negative sign due to {R**(-1/2)}*[H(x+dx)-H(x)]*R**(-1/2)*(y-Hx)} !==================================================================== if(MPIRANK.eq.0) then Kindex=0 do jj=NENS_START,NENS Kindex=Kindex+1 g(jj)=g(jj)-sumx(Kindex) end do end if !==================================================================== if(MPIRANK.eq.0) then do jj=NENS_START,min(20,NENS) write(*,101) jj,g(jj) end do ! write out the gradient (g) rewind gradout write(gradout) (g(K),K=NENS_START,NENS) deallocate(g) end if 101 format('g(',I4,')=',E17.5) !==================================================================== #ifdef MPI_USE CALL MPI_FINALIZE(IERR) #endif stop end