program anlsis_group_obs ! ********************************************************************** ! * . . . ! * PROGRAM: anlsis_group_obs ! * ! * PRGMMR: D. ZUPANSKI ORG: CIRA/CSU DATE: 2004-09-28 ! * ! * ABSTRACT: THIS PROGRAM READS GLOBAL OBSERVATION FILES AND CREATES ! * FILES WITH SMALL NUMBER OF OBSERVATIONS ! * ! * INPUT FILES: ! * obs_${OBSTYPE}.${r3ddate} ! * address_${OBSTYPE}.${r3ddate} ! * first_guess_list.h (include file) ! * first_guess ! ! * OUTPUT FILES: ! * obs_${OBSTYPE}.${r3ddate}_ncateg ! * address_${OBSTYPE}.${r3ddate}_ncateg ! * ! * PROGRAM LOG: ! * ! * 09/28/2004 ..... D. ZUPANSKI ! * 10/12/2004 ..... M. ZUPANSKI: include more options for categories ! * ! ********************************************************************** integer,parameter :: i_loc=22 ! obs address file unit integer,parameter :: i_obs=23 ! obs file unit integer,parameter :: i_loc_new=105 ! new obs address file unit integer,parameter :: i_obs_new=53 ! new obs file unit logical,parameter :: option_0=.true. ! Obs Categories - 0th option (ALL obs in ONE categ) logical,parameter :: option_1=.false. ! Obs Categories - 1st option (just by order) logical,parameter :: option_2=.false. ! Obs Categories - 2nd option (every 5th, 10th obs) logical,parameter :: option_3=.false. ! Obs Categories - 3rd option (x-y grid) logical,parameter :: option_4=.false. ! Obs Categories - 4th option (lon-lat grid) !-------------------------------------------------------------------------------------- real,parameter :: deg_incr=4.50 !!! average grid-point increment (degrees) !---------------------------------------------------------------- !! real,parameter :: deg_incr=2.25 !!! average grid-point increment (degrees) !!!! (Could relate this to the resolution from script: SWMresol )-------------- !-------------------------------------------------------------------------------------- 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*40 :: file_obs character*40 :: file_loc character*40 :: obs_form1,obs_form2,obs_form3,obs_form4 character*40 :: obs_form5,obs_form6,obs_form7,obs_form8 character*40 :: loc_form1,loc_form2,loc_form3,loc_form4 character*40 :: loc_form5,loc_form6,loc_form7,loc_form8 character*40 :: obs_form character*40 :: loc_form character*9 :: fg_name,obs_name character*7 :: obstype character*14 :: rdate real, dimension(:), allocatable :: obs,obs_err real, dimension(:), allocatable :: x_lon,y_lat real, dimension(:,:), allocatable :: categ_obs,categ_obs_err real, dimension(:,:), allocatable :: categ_x,categ_y integer, dimension(:), allocatable :: kloc,kobs integer :: obsjj,locjj integer, dimension(:), allocatable :: i1_loc,i2_loc,i3_loc,i4_loc integer, dimension(:), allocatable :: icateg integer, dimension(:), allocatable :: icateg_start,icateg_end integer, dimension(:), allocatable :: jcateg_start,jcateg_end integer, dimension(:,:), allocatable :: j1_loc,j2_loc,j3_loc,j4_loc integer NENS,NENS_NSTART integer::ncateg,N_degree,N_min integer,parameter::I_degree=6 !------------------------------ NAMELIST /MODEL_DIMENSION/ NNXP,NNYP,NNZP,NNSOIL,NNSNOW NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS !------------------------------ write(*,*) " start group_obs" !--- read ensemble size ------- REWIND 15 READ(15,ENSEMBLE_SIZE) !------------------------------ !--- 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) imax=NNXP(1) jmax=NNYP(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)) do i=1,fgvar_max READ(105,*) fg_name fgvar(i)%name=fg_name write(*,*) fgvar(i)%name end do CLOSE(105,STATUS='KEEP') !-- used for naming the otput files ----------- call getenv("OBSTYPE",obstype) call getenv("rdate",rdate) write(*,*) "obstype=",obstype," rdate=",rdate obs_form1="('obs_',A6,'.',A14,'_',I1)" obs_form2="('obs_',A6,'.',A14,'_',I2)" obs_form3="('obs_',A6,'.',A14,'_',I3)" obs_form4="('obs_',A6,'.',A14,'_',I4)" obs_form5="('obs_',A6,'.',A14,'_',I5)" obs_form6="('obs_',A6,'.',A14,'_',I6)" obs_form7="('obs_',A6,'.',A14,'_',I7)" obs_form8="('obs_',A6,'.',A14,'_',I8)" loc_form1="('address_',A6,'.',A14,'_',I1)" loc_form2="('address_',A6,'.',A14,'_',I2)" loc_form3="('address_',A6,'.',A14,'_',I3)" loc_form4="('address_',A6,'.',A14,'_',I4)" loc_form5="('address_',A6,'.',A14,'_',I5)" loc_form6="('address_',A6,'.',A14,'_',I6)" loc_form7="('address_',A6,'.',A14,'_',I7)" loc_form8="('address_',A6,'.',A14,'_',I8)" !---------------------------------------- !-- read ALL obs ------------------------ rewind i_obs read(i_obs) kk_res allocate(obs(1:kk_res)) allocate(obs_err(1:kk_res)) allocate(x_lon(1:kk_res)) allocate(y_lat(1:kk_res)) !!!! read(i_obs) obs,obs_err read(i_obs) obs,obs_err,x_lon,y_lat write(*,*) " min,max obs=",minval(obs),maxval(obs) write(*,*) " min,max obs_err=",minval(obs_err),maxval(obs_err) write(*,*) " min,max x_lon=",minval(x_lon),maxval(x_lon) write(*,*) " min,max y_lat=",minval(y_lat),maxval(y_lat) write(*,*) "No. of observations=",kk_res !-------------------------------------------------------------------------------------- !--- define ncateg based on the number of ensembles and the number of observations ---- allocate(icateg(1:kk_res)) !------------------ if(option_0) then !------------------ ncateg=1 icateg(1:kk_res)=1 write(*,*) "option=0, ncateg=",ncateg !------------------ elseif(option_1) then !------------------ N_degree=(NENS-NENS_START+1)*I_degree iwork1=int(kk_res/N_degree) iwork2=MOD(kk_res,N_degree) if(iwork2.eq.0) then ncateg=iwork1 else ncateg=iwork1+1 end if write(*,*) "I_degree=",I_degree," Nobs=",kk_res," No. categories=",ncateg allocate(jcateg_start(1:ncateg)) allocate(jcateg_end(1:ncateg)) !-------------------------------------------------------------------------------- !--- This is NOT for MPI !!! ---------------------------------------------------- !--- NORMALLY, para_range is used to distribute PEs ----------------------------- !--- HERE, para_range is used to distribute categories ------------------------- !-------------------------------------------------------------------------------- do irank=1,ncateg call para_range(1,kk_res,ncateg,irank-1,jsta,jend) jcateg_start(irank)=jsta jcateg_end(irank)=jend write(*,*) "irank=",irank," jsta,jend=",jsta,jend end do do irank=1,ncateg icateg(jcateg_start(irank):jcateg_end(irank))=irank end do write(*,*) "option=1, ncateg=",ncateg !------------------ elseif(option_2) then !------------------ N_degree=(NENS-NENS_START+1)*I_degree iwork1=int(kk_res/N_degree) iwork2=MOD(kk_res,N_degree) if(iwork2.eq.0) then ncateg=iwork1 else ncateg=iwork1+1 end if ncateg=2 !!!in case I would like to set-up ncateg do irank=1,ncateg do iobs=irank,kk_res,ncateg icateg(iobs)=irank end do end do write(*,*) "option=2, ncateg=",ncateg !------------------ elseif(option_3) then !------------------ Nx_block=20 !! x-size of a block (grid-point units) Ny_block=20 !! y-size of a block (grid-point units) iwork1=int(imax/Nx_block) iwork2=MOD(imax,Nx_block) if(iwork2.eq.0) then ncateg_x=iwork1 else ncateg_x=iwork1+1 end if iwork1=int(jmax/Ny_block) iwork2=MOD(jmax,Ny_block) if(iwork2.eq.0) then ncateg_y=iwork1 else ncateg_y=iwork1+1 end if ncateg=ncateg_x*ncateg_y write(*,*) "Ncateg_x=",ncateg_x," Ncateg_y=",ncateg_y, & " No. categories=",ncateg !--- get start and end indexes for each category (box) --------------------- do irank=1,ncateg_x call para_range(1,imax,ncateg_x,irank-1,ista,iend) icateg_start(irank)=ista icateg_end(irank)=iend write(*,*) "irank=",irank," ista,iend=",ista,iend end do do irank=1,ncateg_y call para_range(1,jmax,ncateg_y,irank-1,jsta,jend) jcateg_start(irank)=jsta jcateg_end(irank)=jend write(*,*) "irank=",irank," jsta,jend=",jsta,jend end do !----- get icateg(iobs) ------------------------ do iobs=1,kk_res irank=0 do ir=1,ncateg_x if(x_index(iobs).ge.icateg_start(ir).and.x_index(iobs).le.icateg_end(ir)) then do jr=1,ncateg_y if(y_index(iobs).ge.jcateg_start(jr).and.y_index(iobs).le.jcateg_end(jr)) then irank=irank+1 icateg(iobs)=irank exit end if end do end if end do end do !! obs loop write(*,*) "option=3, ncateg=",ncateg !------------------ elseif(option_4) then !------------------ glat_block=20.0 glon_block=20.0 imax=INT(360./deg_incr + 0.001) jmax=INT(180./deg_incr + 0.001) + 1 Nx_block=INT(glon_block/deg_incr) !! x-size of a block (grid-point units) Ny_block=INT(glat_block/deg_incr) !! y-size of a block (grid-point units) iwork1=int(imax/Nx_block) iwork2=MOD(imax,Nx_block) if(iwork2.eq.0) then ncateg_x=iwork1 else ncateg_x=iwork1+1 end if iwork1=int(jmax/Ny_block) iwork2=MOD(jmax,Ny_block) if(iwork2.eq.0) then ncateg_y=iwork1 else ncateg_y=iwork1+1 end if ncateg=ncateg_x*ncateg_y write(*,*) "Ncateg_x=",ncateg_x," Ncateg_y=",ncateg_y," No. categories=",ncateg !--- get start and end indexes for each category (box) --------------------- do irank=1,ncateg_x call para_range(1,imax,ncateg_x,irank-1,ista,iend) icateg_start(irank)=ista icateg_end(irank)=iend write(*,*) "irank=",irank," ista,iend=",ista,iend end do do irank=1,ncateg_y call para_range(1,jmax,ncateg_y,irank-1,jsta,jend) jcateg_start(irank)=jsta jcateg_end(irank)=jend write(*,*) "irank=",irank," jsta,jend=",jsta,jend end do !----- get icateg(iobs) ------------------------ do iobs=1,kk_res irank=0 do ir=1,ncateg_x if(x_index(iobs).ge.icateg_start(ir).and.x_index(iobs).le.icateg_end(ir)) then do jr=1,ncateg_y if(y_index(iobs).ge.jcateg_start(jr).and.y_index(iobs).le.jcateg_end(jr)) then irank=irank+1 icateg(iobs)=irank exit end if end do end if end do end do !! obs loop write(*,*) "option=4, ncateg=",ncateg !------------------ end if !------------------ rewind i_loc !-- allocate for obs categories --------------------------- allocate(categ_obs(1:ncateg,1:kk_res)) allocate(categ_obs_err(1:ncateg,1:kk_res)) allocate(categ_x(1:ncateg,1:kk_res)) !! x-index or lon-coordinate allocate(categ_y(1:ncateg,1:kk_res)) !! y-index or lat-coordinate allocate(kloc(1:ncateg)) allocate(kobs(1:ncateg)) do jj=1,ncateg if(jj.lt.10) then obs_form=obs_form1 loc_form=loc_form1 else if(jj.lt.100) then obs_form=obs_form2 loc_form=loc_form2 else if(jj.lt.1000) then obs_form=obs_form3 loc_form=loc_form3 else if(jj.lt.10000) then obs_form=obs_form4 loc_form=loc_form4 else if(jj.lt.100000) then obs_form=obs_form5 loc_form=loc_form5 else if(jj.lt.1000000) then obs_form=obs_form6 loc_form=loc_form6 else if(jj.lt.10000000) then obs_form=obs_form7 loc_form=loc_form7 else if(jj.lt.100000000) then obs_form=obs_form8 loc_form=loc_form8 else if(jj.ge.1000000000) then write(*,*) " PROBLEM: categ=",jj," Maximum: 100000000" stop endif write(file_loc,fmt=loc_form) obstype,rdate,jj write(*,*) "file_loc=",file_loc locjj=i_loc_new+jj OPEN(locjj,FILE=file_loc,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*) locjj,' OPEN UNIT ERROR IER=',IER end do kobs(:)=0 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 read(i_loc) obs_name 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 !--------------------------------------------------------------------------- allocate(j1_loc(1:ncateg,1:nobs)) allocate(j2_loc(1:ncateg,1:nobs)) allocate(j3_loc(1:ncateg,1:nobs)) allocate(j4_loc(1:ncateg,1:nobs)) !-- kloc(:)=0 do k=1,nobs iobs=iobs+1 kloc(icateg(iobs))=kloc(icateg(iobs))+1 kobs(icateg(iobs))=kobs(icateg(iobs))+1 ii=icateg(iobs) categ_obs(ii,kobs(ii))=obs(iobs) categ_obs_err(ii,kobs(ii))=obs_err(iobs) categ_x(ii,kobs(ii))=x_lon(iobs) categ_y(ii,kobs(ii))=y_lat(iobs) j1_loc(ii,kloc(ii))=i1_loc(k) j2_loc(ii,kloc(ii))=i2_loc(k) j3_loc(ii,kloc(ii))=i3_loc(k) j4_loc(ii,kloc(ii))=i4_loc(k) end do !--- write out for each obs category ---------------------------- do jj=1,ncateg if(jj.lt.10) then obs_form=obs_form1 loc_form=loc_form1 else if(jj.lt.100) then obs_form=obs_form2 loc_form=loc_form2 else if(jj.lt.1000) then obs_form=obs_form3 loc_form=loc_form3 else if(jj.lt.10000) then obs_form=obs_form4 loc_form=loc_form4 else if(jj.lt.100000) then obs_form=obs_form5 loc_form=loc_form5 else if(jj.lt.1000000) then obs_form=obs_form6 loc_form=loc_form6 else if(jj.lt.10000000) then obs_form=obs_form7 loc_form=loc_form7 else if(jj.lt.100000000) then obs_form=obs_form8 loc_form=loc_form8 else if(jj.ge.1000000000) then write(*,*) " PROBLEM: categ=",jj," Maximum: 100000000" stop endif kk_categ=kloc(jj) locjj=i_loc_new+jj write(locjj) obs_name write(locjj) kk_categ write(locjj) j1_loc(jj,1:kk_categ),j2_loc(jj,1:kk_categ), & j3_loc(jj,1:kk_categ),j4_loc(jj,1:kk_categ) end do !! jj-loop, over obs categories !-- deallocate(i1_loc) deallocate(i2_loc) deallocate(i3_loc) deallocate(i4_loc) deallocate(j1_loc) deallocate(j2_loc) deallocate(j3_loc) deallocate(j4_loc) !----------------------------------------------------------- else if(fgvar_list(j)%ndim.eq.3) then read(i_loc) obs_name 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 !--------------------------------------------------------------------------- allocate(j1_loc(1:ncateg,1:nobs)) allocate(j2_loc(1:ncateg,1:nobs)) allocate(j3_loc(1:ncateg,1:nobs)) !-- kloc(:)=0 do k=1,nobs iobs=iobs+1 kloc(icateg(iobs))=kloc(icateg(iobs))+1 kobs(icateg(iobs))=kobs(icateg(iobs))+1 ii=icateg(iobs) categ_obs(ii,kobs(ii))=obs(iobs) categ_obs_err(ii,kobs(ii))=obs_err(iobs) categ_x(ii,kobs(ii))=x_lon(iobs) categ_y(ii,kobs(ii))=y_lat(iobs) j1_loc(ii,kloc(ii))=i1_loc(k) j2_loc(ii,kloc(ii))=i2_loc(k) j3_loc(ii,kloc(ii))=i3_loc(k) end do !--- write out for each obs category ---------------------------- do jj=1,ncateg if(jj.lt.10) then obs_form=obs_form1 loc_form=loc_form1 else if(jj.lt.100) then obs_form=obs_form2 loc_form=loc_form2 else if(jj.lt.1000) then obs_form=obs_form3 loc_form=loc_form3 else if(jj.lt.10000) then obs_form=obs_form4 loc_form=loc_form4 else if(jj.lt.100000) then obs_form=obs_form5 loc_form=loc_form5 else if(jj.lt.1000000) then obs_form=obs_form6 loc_form=loc_form6 else if(jj.lt.10000000) then obs_form=obs_form7 loc_form=loc_form7 else if(jj.lt.100000000) then obs_form=obs_form8 loc_form=loc_form8 else if(jj.ge.1000000000) then write(*,*) " PROBLEM: categ=",jj," Maximum: 100000000" stop endif kk_categ=kloc(jj) locjj=i_loc_new+jj write(locjj) obs_name write(locjj) kk_categ write(locjj) j1_loc(jj,1:kk_categ),j2_loc(jj,1:kk_categ), & j3_loc(jj,1:kk_categ) end do !! jj-loop, over obs categories !-- deallocate(i1_loc) deallocate(i2_loc) deallocate(i3_loc) deallocate(j1_loc) deallocate(j2_loc) deallocate(j3_loc) !--------------------------------------------------------------------------- else if(fgvar_list(j)%ndim.eq.2) then read(i_loc) obs_name read(i_loc) nobs allocate(i1_loc(1:nobs)) allocate(i2_loc(1:nobs)) read(i_loc) i1_loc,i2_loc !--------------------------------------------------------------------------- allocate(j1_loc(1:ncateg,1:nobs)) allocate(j2_loc(1:ncateg,1:nobs)) !-- kloc(:)=0 do k=1,nobs iobs=iobs+1 kloc(icateg(iobs))=kloc(icateg(iobs))+1 kobs(icateg(iobs))=kobs(icateg(iobs))+1 ii=icateg(iobs) categ_obs(ii,kobs(ii))=obs(iobs) categ_obs_err(ii,kobs(ii))=obs_err(iobs) categ_x(ii,kobs(ii))=x_lon(iobs) categ_y(ii,kobs(ii))=y_lat(iobs) j1_loc(ii,kloc(ii))=i1_loc(k) j2_loc(ii,kloc(ii))=i2_loc(k) end do !--- write out for each obs category ---------------------------- do jj=1,ncateg if(jj.lt.10) then obs_form=obs_form1 loc_form=loc_form1 else if(jj.lt.100) then obs_form=obs_form2 loc_form=loc_form2 else if(jj.lt.1000) then obs_form=obs_form3 loc_form=loc_form3 else if(jj.lt.10000) then obs_form=obs_form4 loc_form=loc_form4 else if(jj.lt.100000) then obs_form=obs_form5 loc_form=loc_form5 else if(jj.lt.1000000) then obs_form=obs_form6 loc_form=loc_form6 else if(jj.lt.10000000) then obs_form=obs_form7 loc_form=loc_form7 else if(jj.lt.100000000) then obs_form=obs_form8 loc_form=loc_form8 else if(jj.ge.1000000000) then write(*,*) " PROBLEM: categ=",jj," Maximum: 100000000" stop endif kk_categ=kloc(jj) locjj=i_loc_new+jj write(locjj) obs_name write(locjj) kk_categ write(locjj) j1_loc(jj,1:kk_categ),j2_loc(jj,1:kk_categ) end do !! jj-loop, over obs categories !-- deallocate(i1_loc) deallocate(i2_loc) deallocate(j1_loc) deallocate(j2_loc) !--------------------------------------------------------------------------- else if(fgvar_list(j)%ndim.eq.1) then read(i_loc) obs_name read(i_loc) nobs allocate(i1_loc(1:nobs)) read(i_loc) i1_loc !--------------------------------------------------------------------------- allocate(j1_loc(1:ncateg,1:nobs)) !-- kloc(:)=0 do k=1,nobs iobs=iobs+1 kloc(icateg(iobs))=kloc(icateg(iobs))+1 kobs(icateg(iobs))=kobs(icateg(iobs))+1 ii=icateg(iobs) categ_obs(ii,kobs(ii))=obs(iobs) categ_obs_err(ii,kobs(ii))=obs_err(iobs) categ_x(ii,kobs(ii))=x_lon(iobs) categ_y(ii,kobs(ii))=y_lat(iobs) j1_loc(ii,kloc(ii))=i1_loc(k) end do !--- write out for each obs category ---------------------------- do jj=1,ncateg if(jj.lt.10) then obs_form=obs_form1 loc_form=loc_form1 else if(jj.lt.100) then obs_form=obs_form2 loc_form=loc_form2 else if(jj.lt.1000) then obs_form=obs_form3 loc_form=loc_form3 else if(jj.lt.10000) then obs_form=obs_form4 loc_form=loc_form4 else if(jj.lt.100000) then obs_form=obs_form5 loc_form=loc_form5 else if(jj.lt.1000000) then obs_form=obs_form6 loc_form=loc_form6 else if(jj.lt.10000000) then obs_form=obs_form7 loc_form=loc_form7 else if(jj.lt.100000000) then obs_form=obs_form8 loc_form=loc_form8 else if(jj.ge.1000000000) then write(*,*) " PROBLEM: categ=",jj," Maximum: 100000000" stop endif kk_categ=kloc(jj) locjj=i_loc_new+jj write(locjj) obs_name write(locjj) kk_categ write(locjj) j1_loc(jj,1:kk_categ) end do !! jj-loop, over obs categories !-- deallocate(i1_loc) deallocate(j1_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 deallocate(icateg) !---------------------------------------------- !-- close category address files --------------- do jj=1,ncateg locjj=i_loc_new+jj CLOSE(locjj,status='keep') end do !--- write obs categories ---------------------------- do jj=1,ncateg if(jj.lt.10) then obs_form=obs_form1 loc_form=loc_form1 else if(jj.lt.100) then obs_form=obs_form2 loc_form=loc_form2 else if(jj.lt.1000) then obs_form=obs_form3 loc_form=loc_form3 else if(jj.lt.10000) then obs_form=obs_form4 loc_form=loc_form4 else if(jj.lt.100000) then obs_form=obs_form5 loc_form=loc_form5 else if(jj.lt.1000000) then obs_form=obs_form6 loc_form=loc_form6 else if(jj.lt.10000000) then obs_form=obs_form7 loc_form=loc_form7 else if(jj.lt.100000000) then obs_form=obs_form8 loc_form=loc_form8 else if(jj.ge.1000000000) then write(*,*) " PROBLEM: categ=",jj," Maximum: 100000000" stop endif write(file_obs,fmt=obs_form) obstype,rdate,jj write(*,*) "file_obs=",file_obs OPEN(i_obs_new,FILE=file_obs,FORM='UNFORMATTED',IOSTAT=IER) IF(IER.NE.0)WRITE(*,*) obsjj,' OPEN UNIT ERROR IER=',IER write(i_obs_new) kobs(jj) write(i_obs_new) categ_obs(jj,1:kobs(jj)),categ_obs_err(jj,1:kobs(jj)) !!! write(i_obs_new) categ_obs(jj,1:kobs(jj)),categ_obs_err(jj,1:kobs(jj)), & !!! categ_x (jj,1:kobs(jj)),categ_y (jj,1:kobs(jj)) CLOSE(i_obs_new,status='keep') write(*,*) jj," nobs=",kobs(jj) & ," obs,obs_err min,max=",minval(categ_obs(jj,1:kobs(jj))) & ,maxval(categ_obs(jj,1:kobs(jj))) & ,minval(categ_obs_err(jj,1:kobs(jj))) & ,maxval(categ_obs_err(jj,1:kobs(jj))) end do !! jj - obs categories !---------------------------------------------- rewind 55 write(55,*) ncateg write(*,*) " end group_obs" end program anlsis_group_obs