program anlsis_operator ! ********************************************************************** ! * . . . ! * PROGRAM: anlsis_operator ! * ! * PRGMMR: D. ZUPANSKI, M. Zupanski ORG: CIRA/CSU DATE: 2003-10-01 ! * ! * ABSTRACT: THIS PROGRAM READS FIRST GUESS VARIABLES AND ! * OBSERVATIONS, AND CALCULATES THE RESIDUAL [R**(-1/2)*(y-Hx)] ! * ! * INPUT FILES: ! * fguess ! * obs_${OBSTYPE}.${r3ddate} ! * first_guess_list.h (include file) ! * first_guess ! ! * OUTPUT FILES: ! * outres ! * ! * PROGRAM LOG: ! * ! * 10/01/2003 ..... D. ZUPANSKI ! * 10/07/2003 ..... D. ZUPANSKI, M. Zupanski ! * ! ********************************************************************** integer,parameter :: i_fguess=21 ! first guess file unit integer,parameter :: i_loc=22 ! obs address file unit integer,parameter :: i_obs=23 ! obs file unit integer,parameter :: i_res=51 ! residual file unit integer,parameter :: out_rms1=52 ! rmsobs file unit integer,parameter :: out_rms2=53 ! rmsobs file unit integer,parameter:: infl_unit=55 ! inflation file # INCLUDE 'define_var_list_type.h' INCLUDE 'define_fgvar_type.h' TYPE (var_list_type), dimension(:), allocatable :: fgvar_list TYPE (fgvar_type), dimension(:), allocatable :: fgvar integer, dimension(1) :: NNXP,NNYP,NNZP,NNSOIL,NNSNOW integer :: fgvar_max character*20 :: filename character*9 :: fg_name,obs_name real :: obs_err_inv real :: xx real, dimension(:), allocatable :: res real, dimension(:), allocatable :: obs,obs_err real, dimension(:), allocatable :: rms character*9,dimension(:) :: rms_name real, dimension(:), allocatable :: fg1d real, dimension(:,:), allocatable :: fg2d real, dimension(:,:,:), allocatable :: fg3d real, dimension(:,:,:,:), allocatable :: fg4d character*3 H_obs integer, dimension(:), allocatable :: i1_loc,i2_loc,i3_loc,i4_loc !------------------------------ NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW NAMELIST /HOBS/ H_obs !------------------------------ write(*,*) " start anlsis_operator" !--- get observation operator type CLOSE(201) OPEN(UNIT=201,FILE='H_obs.parm',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(3,*)' 201 OPEN UNIT ERROR IER=',IER READ(201,HOBS) CLOSE(201,STATUS='KEEP') !------------------------------ !--- get model dimensions write(filename,200) 200 format('model_dimension.name') CLOSE(500) OPEN(UNIT=500,FILE=filename,FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0) WRITE(3,*)' 500 OPEN UNIT ERROR IER=',IER READ(500,MODEL_DIMENSION) CLOSE(500,STATUS='KEEP') write(*,*) "NNXP=",NNXP(1)," NNYP=",NNYP(1)," NNZP=",NNZP(1) INCLUDE 'first_guess_list.h' write(*,*) "max_num_of_first_guess=",max_num_of_first_guess do nfg=1,max_num_of_first_guess write(*,*) "ndim ",fgvar_list(nfg)%ndim write(*,*) "start1 ",fgvar_list(nfg)%start_index(1) write(*,*) "start2 ",fgvar_list(nfg)%start_index(2) write(*,*) "start3 ",fgvar_list(nfg)%start_index(3) write(*,*) "start4 ",fgvar_list(nfg)%start_index(4) write(*,*) "end1 ",fgvar_list(nfg)%end_index(1) write(*,*) "end2 ",fgvar_list(nfg)%end_index(2) write(*,*) "end3 ",fgvar_list(nfg)%end_index(3) write(*,*) "end4 ",fgvar_list(nfg)%end_index(4) write(*,*) "name ",fgvar_list(nfg)%name write(*,*) "description ",fgvar_list(nfg)%description end do CLOSE(105) OPEN(105,FILE='first_guess_anlsis',FORM='FORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*)' 105 OPEN UNIT ERROR IER=',IER READ(105,*) fgvar_max allocate(fgvar(1:fgvar_max)) allocate(rms(1:fgvar_max)) rms(:)=0.0 allocate(rms_name(1:fgvar_max)) do i=1,fgvar_max READ(105,*) fg_name fgvar(i)%name=fg_name write(*,*) fgvar(i)%name end do CLOSE(105,STATUS='KEEP') !---------------------------------------------- rewind i_obs read(i_obs) kk_res !---------------------------------------------- !---------------------------------------------- allocate(res(1:kk_res)) allocate(obs(1:kk_res)) allocate(obs_err(1:kk_res)) read(i_obs) obs,obs_err ! alfa=1.0 ! CLOSE(infl_unit) ! OPEN(UNIT=infl_unit,FILE='infl_file',FORM='FORMATTED',IOSTAT=IER) ! read(infl_unit,*,end=230) alfa ! CLOSE(infl_unit) !230 obs_err(:)=obs_err*alfa write(*,*) " min,max obs=",minval(obs),maxval(obs) write(*,*) " min,max obs_err=",minval(obs_err),maxval(obs_err) rewind i_fguess rewind i_loc iobs=0 nvar=0 icount=0 do i=1,fgvar_max write(*,*) "do i=",i do j=1,max_num_of_first_guess write(*,*) "do j=",j if(fgvar(i)%name.eq.fgvar_list(j)%name) then nvar=nvar+1 i1s=fgvar_list(j)%start_index(1) i2s=fgvar_list(j)%start_index(2) i3s=fgvar_list(j)%start_index(3) i4s=fgvar_list(j)%start_index(4) i1e=fgvar_list(j)%end_index(1) i2e=fgvar_list(j)%end_index(2) i3e=fgvar_list(j)%end_index(3) i4e=fgvar_list(j)%end_index(4) !-- read first guess on unit i_fguess if(fgvar_list(j)%ndim.eq.4) then allocate(fg4d(i1s:i1e,i2s:i2e,i3s:i3e,i4s:i4e)) read(i_fguess) fg_name read(i_fguess) fg4d(i1s:i1e,i2s:i2e,i3s:i3e,i4s:i4e) read(i_loc,end=220) obs_name if(obs_name.ne.fg_name) then backspace i_loc exit end if read(i_loc) nobs allocate(i1_loc(1:nobs)) allocate(i2_loc(1:nobs)) allocate(i3_loc(1:nobs)) allocate(i4_loc(1:nobs)) read(i_loc) i1_loc,i2_loc,i3_loc,i4_loc do k=1,nobs iobs=iobs+1 obs_err_inv=1./obs_err(iobs) i1=i1_loc(k) i2=i2_loc(k) i3=i3_loc(k) i4=i4_loc(k) call Hx(H_obs,fg4d(i1,i2,i3,i4),xx) res(iobs)=obs_err_inv*(obs(iobs)-xx) rms(nvar)=rms(nvar)+(obs(iobs)-xx)**2 end do if(nobs.gt.0) then rms(nvar)=sqrt(rms(nvar)/nobs) endif rms_name(nvar)=obs_name deallocate(fg4d) deallocate(i1_loc) deallocate(i2_loc) deallocate(i3_loc) deallocate(i4_loc) else if(fgvar_list(j)%ndim.eq.3) then allocate(fg3d(i1s:i1e,i2s:i2e,i3s:i3e)) read(i_fguess) fg_name read(i_fguess) fg3d(i1s:i1e,i2s:i2e,i3s:i3e) write(*,*)" fg_name=",fg_name write(*,*)" min,max fg3d=",minval(fg3d),maxval(fg3d) read(i_loc,end=220) obs_name write(*,*)" 1 obs_name=",obs_name if(obs_name.ne.fg_name) then backspace i_loc exit end if read(i_loc) nobs allocate(i1_loc(1:nobs)) allocate(i2_loc(1:nobs)) allocate(i3_loc(1:nobs)) read(i_loc) i1_loc,i2_loc,i3_loc do k=1,nobs iobs=iobs+1 obs_err_inv=1./obs_err(iobs) i1=i1_loc(k) i2=i2_loc(k) i3=i3_loc(k) call Hx(H_obs,fg3d(i1,i2,i3),xx) res(iobs)=obs_err_inv*(obs(iobs)-xx) rms(nvar)=rms(nvar)+(obs(iobs)-xx)**2 end do if(nobs.gt.0) then rms(nvar)=sqrt(rms(nvar)/nobs) endif rms_name(nvar)=obs_name write(*,*) " 2 obs_name=",obs_name," nobs=",nobs," iobs=",iobs deallocate(fg3d) deallocate(i1_loc) deallocate(i2_loc) deallocate(i3_loc) else if(fgvar_list(j)%ndim.eq.2) then allocate(fg2d(i1s:i1e,i2s:i2e)) read(i_fguess) fg_name read(i_fguess) fg2d(i1s:i1e,i2s:i2e) read(i_loc,end=220) obs_name if(obs_name.ne.fg_name) then backspace i_loc exit end if read(i_loc) nobs allocate(i1_loc(1:nobs)) allocate(i2_loc(1:nobs)) read(i_loc) i1_loc,i2_loc do k=1,nobs iobs=iobs+1 obs_err_inv=1./obs_err(iobs) i1=i1_loc(k) i2=i2_loc(k) call Hx(H_obs,fg2d(i1,i2),xx) res(iobs)=obs_err_inv*(obs(iobs)-xx) rms(nvar)=rms(nvar)+(obs(iobs)-xx)**2 end do if(nobs.gt.0) then rms(nvar)=sqrt(rms(nvar)/nobs) endif rms_name(nvar)=obs_name deallocate(fg2d) deallocate(i1_loc) deallocate(i2_loc) else if(fgvar_list(j)%ndim.eq.1) then allocate(fg1d(i1s:i1e)) read(i_fguess) fg_name read(i_fguess) fg1d(i1s:i1e) read(i_loc,end=220) obs_name if(obs_name.ne.fg_name) then backspace i_loc exit end if read(i_loc) nobs allocate(i1_loc(1:nobs)) read(i_loc) i1_loc do k=1,nobs iobs=iobs+1 obs_err_inv=1./obs_err(iobs) i1=i1_loc(k) call Hx(H_obs,fg1d(i1),xx) res(iobs)=obs_err_inv*(obs(iobs)-xx) rms(nvar)=rms(nvar)+(obs(iobs)-xx)**2 write(*,*) "fg1d,xx=",fg1d(i1),xx," obs,err=",obs(iobs),obs_err(iobs) end do if(nobs.gt.0) then rms(nvar)=sqrt(rms(nvar)/nobs) endif rms_name(nvar)=obs_name deallocate(fg1d) deallocate(i1_loc) else write(*,*) "FGUESS DOES NOT HAVE APPROPRIATE ndim=",fgvar_list(j)%ndim end if !---------------------------------------------- end if !!!if(fgvar(i)%name.eq.fgvar_list(j)%name) then end do end do 220 continue !---------------------------------------------- !-- write out residuals ! do k=1,kk_res ! write(*,*) k,res(k) ! end do write(*,*) " kk_res=",kk_res," min,max res=",minval(res),maxval(res) if(iobs.ne.kk_res) then write(*,*) "PROBLEM iobs=",iobs," kk_res=",kk_res else rewind i_res WRITE(i_res) kk_res WRITE(i_res) res end if 100 format(6E20.10) write(*,*)" rms_name=",rms_name k_start=1 k_end=min(3,fgvar_max) write(out_rms1,100) (rms(kk),kk=k_start,k_end) if (fgvar_max.gt.k_end) then k_start=k_end+1 k_end=fgvar_max write(out_rms2,100) (rms(kk),kk=k_start,k_end) end if write(*,*) " end anlsis_operator" end program anlsis_operator