program cost_filt_grad ! !!========================================================== !!! THIS IS NOT UPDATED!!! ONLY AN OLD VERSION FROM 4DVAR. !!========================================================== !====== COST-FILT GRADIENT CALCULATION USING ENSEMBLES ================= !------- calculate digital filter inner product --------- ! ! INPUT: XENS = M(x+delx)-FM(x+delx) ! INPUT: XCNTRL = M(x)-FM(x) ! ! F - (Nonlinear) Digital Filter operator ! T -1 ! OUTPUT: g(i) = (x-y) P y - gradient in ensemble space ! ! OUTPUT: FPEN - digital filter cost function ! !-------------------------------------------------------------------- include 'mpif.h' INCLUDE 'ntype.h' !!! NTYPE= No. of cntrl vrbl types !----- integer,parameter::INCR_C=21 ! control fcst increment (y) integer,parameter::INCR_E=22 ! ensemble fcst incr (x) integer,parameter::gradout=51 ! digital filter gradient integer,parameter::IPEN=52 ! digital filter cost-function !----- !--------- non-dimensional obs perturbations i and j ----------- real,allocatable::XCNTRL(:,:,:) real,allocatable::XENS(:,:,:) real,allocatable::FPEN(:) !--------- inner-products (temporary) --------- real,allocatable::sumx(:) real,allocatable::sum_loc(:) !--------- gradient ---------------------- real,allocatable::g(:) integer,allocatable::jlen(:),jdisp(:) integer IERR integer NENS,NTOTAL integer NENS_START CHARACTER*50 FILNAMEI,FILNAMEJ CHARACTER*50 FILNAME CHARACTER*4 ENVIR INCLUDE 'CNTRL4D.comm' INCLUDE 'namelist.h' !!! define namelists (all) ! ! DECLARE NAMELIST ! NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS !==============start calculation=================== ! ! START MPI ! CALL MPI_INIT(IERR) IF (IERR .NE. MPI_SUCCESS) STOP 'FAILED TO INIT MPI' ! Note: NPROC is same as NPES in 4dvar ! Note: MPIRANK is same as MYPE in 4dvar 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) !--- read input files 'cntrl_index.name' (33) and 'cntrl_specs.name' (34) INCLUDE 'read_namelist' ! read all namelists !---------------------------------------- !--- define model space indexes (imeta,jmeta)--------- !---------------------------------------- imeta=IM jmeta=JM !================================================ NEIG=15 if(MPIRANK.eq.0) then write(*,*) "RANDOM PTRB: start reading NEIG=",NEIG end if READ(NEIG,ENSEMBLE_SIZE) NTOTAL=NENS if(MPIRANK.eq.0) then write(*,*) "RANDOM_PTRB: number of ensemble members=",NTOTAL end if CALL INITCTRL !==================================================================== !=== create gradient component for each obs time =================== !==================================================================== NSIZE=NTOTAL !== !== define local dimensions and indexes !== allocate(jdisp(0:NPROC-1)) allocate(jlen(0:NPROC-1)) do irank=0,NPROC-1 CALL para_range(1,NSIZE,NPROC,irank,jsta,jend) !!! write(*,*) irank," jsta=",jsta," jend=",jend jlen(irank)=jend-jsta+1 jdisp(irank)=jsta-1 end do CALL para_range(1,NSIZE,NPROC,MPIRANK,jsta,jend) CALL MPI_BARRIER(MPI_COMM_WORLD,IERR) !== N_LOC=jlen(MPIRANK) write(*,*) MPIRANK," jsta=",jsta," jend=",jend," N_LOC=",N_LOC if(MPIRANK.eq.0) then do NN=0,NPROC-1 write(*,*) MPIRANK," jlen=",jlen(NN)," jdisp=",jdisp(NN) end do end if !----------------------------------------------------------- if(MPIRANK.eq.0) then !!-- g = total gradient in ensemble space (sum over all obs times) allocate(g(NSIZE)) g=0.0 !---- define cost-function ---------------- ALLOCATE(FPEN(NLEV)) FPEN=0.0 FUNTOT=0.0 end if !----------------------------------------------------------- IF(EPSL1.GT.0) THEN !! start calculation !----------------------------------------------------------- allocate(XCNTRL(1:imeta,1:jmeta,1:NLEV)) allocate(XENS(1:imeta,1:jmeta,1:NLEV)) ! !--- read input vector = x ! CLOSE(INCR_C) OPEN(INCR_C,file='oldfilt_incr',status='unknown',form='unformatted') READ(INCR_C) (((XCNTRL(i,j,k),i=1,imeta),j=1,jmeta),k=1,NLEV) CLOSE(INCR_C,STATUS='KEEP') !! write(*,*)MPIRANK," READ: XCNTRL min,max=",minval(XCNTRL),maxval(XCNTRL) CALL MAT_VEC (imeta,jmeta,NLEV,XCNTRL) !! write(*,*)MPIRANK," MATVEC: XCNTRL min,max=",minval(XCNTRL),maxval(XCNTRL) !==================================================================== !! MPI is used for inner-products ! !!-- g = total gradient in ensemble space (sum over all obs times) !!-- sumx = partial gradient (at given obs time) allocate(sumx(NSIZE)) allocate(sum_loc(N_LOC)) kk=0 do JJ=jsta,jend kk=kk+1 write(*,*) "current ensemble ",JJ ! ! read in the ensemble perturbations ! WRITE(FILNAMEJ,1110) JJ 1110 format('oldfilt_incr.',i3.3) CLOSE(INCR_E) OPEN(INCR_E,file=FILNAMEJ,status='unknown',form='unformatted') READ(INCR_E) (((XENS(i,j,k),i=1,imeta),j=1,jmeta),k=1,NLEV) CLOSE(INCR_E,STATUS='KEEP') !! write(*,*)MPIRANK,JJ," READ: XENS min,max=",minval(XENS),maxval(XENS) CALL MAT_VEC (imeta,jmeta,NLEV,XENS) !! write(*,*)MPIRANK,JJ," MATVEC: XENS min,max=",minval(XENS),maxval(XENS) !---- inner product matrix -------- sum_loc(kk)=0.0 do k=1,NLEV do j=1,jmeta do i=1,imeta sum_loc(kk)=sum_loc(kk)+(XENS(i,j,k)-XCNTRL(i,j,k))*XCNTRL(i,j,k) end do end do end do end do !! end JJ (ensemble) loop ! !! write(*,*) MPIRANK," N_LOC=",N_LOC," kk=",kk ! sumx=0.0 ! CALL MPI_GATHERV & (sum_loc,N_LOC,MPI_REAL,sumx,jlen,jdisp,MPI_REAL,0,MPI_COMM_WORLD,IERR) ! deallocate(sum_loc) deallocate(XENS) !! END MPI !================================================================= !--- CALCULATE VIRTUAL PENALTY COST FUNCTION if(MPIRANK.EQ.0) THEN write(3,*) "========================" write(3,*) "GW PENALTY COST FUNCTION" write(3,*) "========================" write(*,*) "========================" write(*,*) "GW PENALTY COST FUNCTION" write(*,*) "========================" !------------------------------------------------------------ !---- define FPEN as a cost-function ---------------- do k=1,NLEV FPEN(k)=0.0 do j=1,jmeta do i=1,imeta FPEN(k)=FPEN(k) + 0.5*XCNTRL(i,j,k)*XCNTRL(i,j,k) end do end do end do FUNCTOT=0.0 do NL=1,NLEV FUNCTOT=FUNCTOT+FPEN(NL) end do end if !!! end rank=0 deallocate(XCNTRL) !==================================================================== ! GRADIENT g IN ENSEMBLE SPACE !==================================================================== if(MPIRANK.eq.0) then ! Kindex=0 do jj=1,NTOTAL Kindex=Kindex+1 g(Kindex)=g(Kindex)+sumx(Kindex) end do ! end if !==================================================================== deallocate(sumx) !----------------------------------------------------------- ENDIF !! end calculation !----------------------------------------------------------- !--- WRITE OUT THE PENALTY COST FUNCTION AND GRADIENT -------- if(MPIRANK.eq.0) then write(*,*) "DIGITAL FILTER GRADIENT: EPSL1=",EPSL1 write(3,*) "DIGITAL FILTER GRADIENT: EPSL1=",EPSL1 do jj=1,10 write(*,101) jj,g(jj) write(3,101) jj,g(jj) end do 101 format('g(',I3,')=',E17.5) write(*,*) "DIGITAL FILTER COST FUNCTION: EPSL1=",EPSL1 WRITE(*,*) " " WRITE(*,*) "DIGITAL FILTER COST-FUNCTION" WRITE(*,1015) FUNCTOT write(3,*) "DIGITAL FILTER COST FUNCTION: EPSL1=",EPSL1 WRITE(3,*) " " WRITE(3,*) "DIGITAL FILTER COST-FUNCTION" WRITE(3,1015) FUNCTOT 1015 FORMAT('FPENTOT=',E15.5) CLOSE(IPEN) OPEN(IPEN,file='costp',status='unknown',form='unformatted') WRITE(IPEN) FUNCTOT,FPEN CLOSE(IPEN,STATUS='KEEP') close(gradout) open(gradout,file='gradnew_filt',status='unknown',form='unformatted') write(gradout) (g(K),K=1,NTOTAL) close(gradout,STATUS='KEEP') deallocate(FPEN) deallocate(g) endif !==================================================================== CALL MPI_FINALIZE(IERR) stop end