!-----------------------------------------------------------------------------! ! Group on Data Assimilation Development - GDAD/CPTEC/INPE ! !-----------------------------------------------------------------------------! !BOP ! ! ! TITLE: GSI Conventional Diagnostic Documentation \\ Version 0.01 ! ! !AUTHORS: Jo\~a Gerd Z. de Mattos ! ! !AFFILIATION: Group on Data Assimilation Development, CPTEC/INPE, Cachoeira Paulista - SP ! ! !DATE: May 15, 2012 ! ! !INTRODUCTION: Package Overview ! ! diag_GSI_conv is a Fortran 90 collection of routines/functions for ! . ! . ! . ! !EOP !--------------------------------------------------------------------- !BOP !--------------------------------------------------------------------- ! !ROUTINE: diag_gsi_cov.f90 !--------------------------------------------------------------------- ! INPE/CPTEC Data Assimilation Group !--------------------------------------------------------------------- ! ! !REVISION HISTORY: ! 15 May 2012 - J. G. de Mattos - Initial Version ! ! ! !--------------------------------------------------------------------- ! PROGRAM diag_GSI_conv ! ! !USES: ! use kinds, only: r_kind,r_single,i_kind use m_utils use m_string IMPLICIT NONE ! ! Auxiliary variables for diagnostic file reading ! character(len=3) :: var character(8),allocatable,dimension(:) :: cdiagbuf real(r_single),allocatable,dimension(:,:) :: rdiagbuf real(r_single),allocatable,dimension(:,:) :: diag integer(i_kind) :: idate integer(i_kind) :: nchar,ninfo,nobs,mype ! ! misc. ! real(r_single) :: zero_single real(r_single) :: tiny_single real :: small_num character(8) :: stationID character :: ch character(len=1024) :: fname, ctlfname, mapfname integer :: lu integer :: idx integer :: iflg integer :: i,j,k, P, V, ic,ttt integer :: HSN(2) integer :: ios character(len=50) :: rowfmt character(len=3) :: DataType character(len=10) :: TimeFile character(len=1024) :: tmp INTEGER :: nymd INTEGER :: nhms real :: Tim integer :: Lev integer :: Flag real, allocatable :: rtmp(:) logical :: OK ! ! Variaveis ! type(data), dimension(nvars+1) :: conv ! !EOP !------------------------------------------------------------------- !BOC ! if ( iargc() .lt. 2) then write(*,'(a)')'error: diag_gsi.x [ YYYYMMDDHH ] [ 01/02 ]' write(*,'(a)')' 01 - Process Diag FGS files' write(*,'(a)')' 02 - Process Diag ANL files' call exit() endif call getarg(1,TimeFile) ! Data YYYYMMDDHH read(TimeFile,'(I8.8,I8.8)') nymd,nhms IF(nhms.LT.100) nhms = nhms * 10000 IF(nhms.GT.100.AND.nhms.LT.10000) nhms = nhms * 100 call getarg(2,DataType) ! 01- Guess, 02- Analysis ! year = INT ( ymd / 10000 ) ! month = MOD ( ymd, 10000 ) / 100 ! day = MOD ( ymd, 100 ) ! ! select case(hms) ! case( : 99 ) ! hour = hms ! minute = 0 ! second = 0 ! case(100:9999) ! hour = INT ( hms / 100 ) ! minute = MOD ( hms, 10 ) ! second = 0 ! case default ! hour = INT ( hms / 10000 ) ! minute = MOD ( hms, 10000 ) / 100 ! second = MOD ( hms, 100 ) ! end select call readcard() ! ! lim_qm: Valores baseados na descricao contida no manual do GSI. ! Estes valores são utilizados para controlar os dados aceitos ou ! rejeitados. ! The parameter lim_qm is a threshold set in GSI read in procedure. The current values of ! lim_qm in GSI are listed in the following table: ! ! +----------------------+---------------+---------------+ ! |The value of namelist | lim_qm for Ps | lim_qm others | ! |option noiqc | | | ! +----------------------+---------------+---------------+ ! |True (without OI QC) | 7 | 8 | ! +----------------------+---------------+---------------+ ! |False (with OI QC) | 4 | 4 | ! +----------------------+---------------+---------------+ ! do i=1,nvars conv(i)%name = vname(i) if (noiqc) then conv(i)%lim_qm = 8 if(vname(i).eq." ps")conv(i)%lim_qm = 7 else conv(i)%lim_qm = 4 endif enddo ! ! ! tiny_single = tiny(zero_single) small_num = 1.e3 * tiny_single ! ! Abrindo e lendo arquivos de diagnostico gerados pelo GSI ! Se o GSI foi executado em paralelo são abertos mais do que um arquivo. Isto ! é controlado pela variaveis nproc. lu = 10 print*, '>>>>>>>> NPROC: ',nproc DO I= 0,nproc-1 FName = TRIM(DirIN)//'/'//TRIM(FInp) CALL str_template(FName, nymd,nhms) call replace_(FName, "%e", int2str(I,'(I4.4)')) call replace_(FName, "%type", TRIM(DataType)) #ifdef SUMMARY print*,trim(fname) #endif OPEN ( UNIT = lu+I, & FILE = trim(fname), & STATUS = 'OLD', & IOSTAT = ios, & ACCESS = 'SEQUENTIAL', & FORM = 'UNFORMATTED', CONVERT='SWAP' ) if(ios > 0 ) then write(*,*) ' file is unavailabe: ', trim(fname) stop 123 endif ENDDO ! ! Reading observation diagnostics ! ! read(lu, ERR=997) idate ! time of observation analisys ! ! ! zero fill to counters ! do i=0,nlev conv(:)%use ( i ) = 0 conv(:)%nuse ( i ) = 0 conv(:)%rej ( i ) = 0 conv(:)%mon ( i ) = 0 Do j=1,2 conv(:)%rmse ( i,j ) = 0 conv(:)%vies ( i,j ) = 0 enddo enddo conv(:)%count(1)=0 conv(:)%count(2)=0 conv(:)%nobs = 0 ttt=0 idx=1 var='tt' mype=0 FILES : DO P=0,nproc-1 ! if(P.eq.93.or.& ! P.eq.99.or.& ! P.eq.100.or.& ! P.ge.101)cycle #ifdef DEBUG write(*,)'=======================',p,'=======================' #endif VARIABLES : DO ! print*,conv(idx)%nobs,var,mype read(lu+P, iostat=ios,end=110) var, nchar,ninfo,nobs,mype ! if(ios/=0)read(lu+P, iostat=ios,end=110) ! if(ios/=0)then print*,var, nchar,ninfo,nobs,mype,ios ! cycle ! endif idx = maxloc(index(vname(:),var),1) #ifdef DEBUG write(*,'(2I,A,5I)')P,conv(idx)%nobs,var, nchar,ninfo,nobs,mype, ios #endif ! if (idx.eq.5)conv(idx)%nobs=ttt ! if(idx.eq.5)print*,conv(idx)%nobs, nobs, conv(idx)%nobs+nobs if (nobs > 0) then print*, "1",nobs allocate(cdiagbuf(nobs),rdiagbuf(ninfo,nobs)) print*, "2",nobs read(lu+P,ERR=999,end=109) cdiagbuf, rdiagbuf STATION : do i=1,nobs j = conv(idx)%nobs + i k = minloc(rdiagbuf(6,i)-levs,mask=(rdiagbuf(6,i)-levs).GE.0,DIM=1) !if (mype.eq.94)print*,i,j, ii, conv(idx)%nobs, idx, conv(:)%nobs, conv(5)%nobs if (var .eq. " ps") k = 1 conv(idx)%kx (j) = nint(rdiagbuf(1,i)) ! observation type conv(idx)%lat (j) = rdiagbuf(3,i) ! observation latitude (degrees) conv(idx)%lon (j) = rdiagbuf(4,i) ! observation longitude (degrees) conv(idx)%lev (j) = levs(k) ! observation level reference conv(idx)%prs (j) = rdiagbuf(6,i) ! observation pressure (hPa) conv(idx)%time (j) = rdiagbuf(8,i) * 60 ! obs time (minutes relative to analysis time) conv(idx)%pbqc (j) = rdiagbuf(9,i) ! input prepbufr qc or event mark conv(idx)%iuse (j) = int(rdiagbuf(12,i)) ! analysis usage flag (1=use, -1=monitoring ) conv(idx)%iusev(j) = int(rdiagbuf(11,i)) ! analysis usage flag ( value ) conv(idx)%robs (j) = rdiagbuf(17,i) ! observation conv(idx)%diff (j) = rdiagbuf(18,i) ! obs-ges used in analysis (K) conv(idx)%rmod (j) = conv(idx)%robs(j)-conv(idx)%diff(j) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! Este passo causa falha de segmentacao quando rdiagbuf muito proximo a zero !! Necessita Investigacao ! ! conv(idx)%sigo (j) = udef ! ! if (rdiagbuf(16,i) > small_num) then! final inverse observation error (K**-1) ! conv(idx)%sigo(j) = 1.0/rdiagbuf(16,i) ! end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! get station ID ! stationID = cdiagbuf(i) iflg = 0 do ic=8,1,-1 ch = stationID(ic:ic) if (ch > ' ' .and. ch <= 'z') then iflg = 1 else stationID(ic:ic) = ' ' end if if (ch == ' ' .and. iflg == 1) then stationID(ic:ic) = '_' endif enddo conv(idx)%ID(j) = stationID ! ! Counting number of observations accepted, reject and monitored ! if (conv(idx)%iuse(j).eq. 1) conv(idx)%use(k) = conv(idx)%use(k) + 1 if (conv(idx)%iuse(j).eq.-1) conv(idx)%nuse(k) = conv(idx)%nuse(k) + 1 ! ! The QC process creates a number indicating the data quality for each observation. ! These numbers are called QC markers in PrepBUFR files and are important as parts of the ! observation information. GSI uses QC markers to decide how to use the data. A brief ! summary of the meaning of the QC markers is as follows: ! ! +-----------------+-----------------------------------------------------------+ ! | QC markes range | Data Process in GSI | ! +-----------------+-----------------------------------------------------------+ ! | > 15 or |GSI skips these observations during reading procedure. That| ! | <= 0 |means these observations are tossed | ! +-----------------+-----------------------------------------------------------+ ! | >= lim_qm |These observations will be in monitoring status. That means| ! | and |these observations will be read in and be processed through| ! | < = 15 |GSI QC process (gross check) and innovation calculation | ! | |stage but will not be used in inner iteration. | ! +-----------------+-----------------------------------------------------------+ ! | > 0 |Observations will be used in further gross check (failure | ! | and |observation will be list in rejection), innovation | ! | < lim_qm |caalculation, and the analysis (inner iteration). | ! +-----------------+-----------------------------------------------------------+ ! if (conv(idx)%iuse(j).eq.-1 .and.& (conv(idx)%pbqc(j) > 15.0 .or.& conv(idx)%pbqc(j) <= 0.0)) conv(idx)%rej(k) = conv(idx)%rej(k) + 1 if (conv(idx)%iuse(j).eq.-1 .and.& (conv(idx)%pbqc(j) > 0.0 .and.& conv(idx)%pbqc(j) < conv(idx)%lim_qm)) conv(idx)%rej(k) = conv(idx)%rej(k) + 1 if (conv(idx)%iuse(j).eq.-1.and.& (conv(idx)%pbqc(j) >= conv(idx)%lim_qm .and.& conv(idx)%pbqc(j) <= 15.0)) conv(idx)%mon(k) = conv(idx)%mon(k) + 1 ! ! Data adjustments ! if (var == "gps")then conv(idx)%robs (j) = rdiagbuf(17,i) ! observation (data(igps,i)) conv(idx)%incref (j) = rdiagbuf(5,i) ! incremental refractivity (x100 %)---- (data(igps,i)-nrefges)/data(igps,i) conv(idx)%diff (j) = rdiagbuf(17,i)*rdiagbuf(5,i) ! obs-ges used in analysis conv(idx)%diffpt (j) = conv(idx)%diff(j)*100 ! obs-ges used in analysis *100% conv(idx)%rmod (j) = conv(idx)%robs(j)-conv(idx)%diff(j) ! ges refrativity (nrefges) conv(idx)%elevgps (j) = rdiagbuf(7,i) ! height in meters/ impact height in meters conv(idx)%kmelev (j) = conv(idx)%elevgps(j)/1000 ! height in km/ impact height in km conv(idx)%hobl (j) = rdiagbuf(19,i) ! model vertical grid (midpoint) conv(idx)%errf (j) = rdiagbuf(16,i) ! final inverse observation error due to superob factor (N**-1) conv(idx)%normdiff (j) = conv(idx)%diff(j)*conv(idx)%errf(j) ! normalized bending angles differences endif ! ! When the data is q, unit convert kg/kg -> g/kg **/ ! if (var == " q") then conv(idx)%robs(j) = conv(idx)%robs(j) * 1000.0 conv(idx)%rmod(j) = conv(idx)%rmod(j) * 1000.0 conv(idx)%diff(j) = conv(idx)%diff(j) * 1000.0 ! conv(idx)%sigo(j) = conv(idx)%sigo(j) * 1000.0 end if ! ! When the data is pw, replase the rprs to udef **/ ! if (var .eq. " pw") conv(idx)%time(j) = udef ! ! ! ! ! When the data is uv U return on uv idx and V return nvars + 1 ! if ( var .eq. " uv")then ! ! UWIND ! conv(idx)%robs(j) = rdiagbuf(17,i) conv(idx)%diff(j) = rdiagbuf(18,i) conv(idx)%rmod(j) = conv(idx)%robs(j)-conv(idx)%diff(j) ! ! VWIND ! conv(nvars+1)%robs(j) = rdiagbuf(20,i) conv(nvars+1)%diff(j) = rdiagbuf(21,i) conv(nvars+1)%rmod(j) = conv(nvars+1)%robs(j)-conv(nvars+1)%diff(j) ! ! SPD ! conv(idx)%robs(j) = sqrt(conv(idx)%robs(j)**2 + conv(nvars+1)%robs(j)**2) conv(idx)%rmod(j) = sqrt(conv(idx)%rmod(j)**2 + conv(nvars+1)%rmod(j)**2) conv(idx)%diff(j) = conv(idx)%robs(j)-conv(idx)%rmod(j) endif if(conv(idx)%robs(j) .gt. 1.0e8) then conv(idx)%robs(j) = udef conv(idx)%diff(j) = udef endif ! ! Statistical calculations ! if (conv(idx)%diff(j).ne.udef)then if (conv(idx)%lat(j) .gt. 0.0)then conv(idx)%vies(k,1) = conv(idx)%vies(k,1) + conv(idx)%diff(j) conv(idx)%rmse(k,1) = conv(idx)%rmse(k,1) + conv(idx)%diff(j)*conv(idx)%diff(j) conv(idx)%count(1)=conv(idx)%count(1)+1 elseif(conv(idx)%lat(j) .lt. 0.0)then conv(idx)%vies(k,2) = conv(idx)%vies(k,2) + conv(idx)%diff(j) conv(idx)%rmse(k,2) = conv(idx)%rmse(k,2) + conv(idx)%diff(j)*conv(idx)%diff(j) conv(idx)%count(2)=conv(idx)%count(2)+1 endif endif ! ! ! enddo STATION 109 continue deallocate(cdiagbuf,rdiagbuf) else read(lu+P) endif conv(idx)%nobs = conv(idx)%nobs + nobs END DO VARIABLES 110 continue CLOSE(lu+P) END DO FILES ! ! Adding observations at all levels ! do i=1,nvars conv(i)%use (0) = sum(conv(i)%use (1:nlev)) conv(i)%nuse(0) = sum(conv(i)%nuse(1:nlev)) conv(i)%rej (0) = sum(conv(i)%rej (1:nlev)) conv(i)%mon (0) = sum(conv(i)%mon (1:nlev)) do j=1,2 conv(i)%vies(0,j) = sum(conv(i)%vies(1:nlev,j)) conv(i)%rmse(0,j) = sum(conv(i)%rmse(1:nlev,j)) conv(i)%vies( 0,j) = conv(i)%vies( 0,j) / (conv(i)%count(j)*nlev) conv(i)%vies(1:nlev,j) = conv(i)%vies(1:nlev,j) / conv(i)%count(j) conv(i)%rmse( 0,j) = sqrt( conv(i)%rmse( 0,j) / (conv(i)%count(j)*nlev) ) conv(i)%rmse(1:nlev,j) = sqrt( conv(i)%rmse(1:nlev,j) / conv(i)%count(j) ) enddo enddo ! ! writing output to plot on GMT and gnuplot ! ! ! Line and bar Graphics ! ! ! write to file ! IF(DataType.EQ."01") DataType="fgs" ! It'll work till 6 outerloops IF(DataType.EQ."02" .OR. DataType.EQ."03" .OR. DataType.EQ."04" .OR. DataType.EQ."05" .OR. DataType.EQ."06") DataType="anl" ! ! Spatial Output ! DO v=1,nvars FName = TRIm(DirOu)//'/'//TRIM(FDatSSO) call replace_(FName, "%obstype", trim(adjustl(vname(v)))) call replace_(FName, "%type", TRIM(DataType)) call str_template(FName,nymd,nhms) Open(Unit = 10, & File = trim(FName),& Status = 'unknown',& Form = 'Unformatted',& access = 'sequential'& ! recordtype = 'stream'& ) ! print*,trim(adjustl(vname(v))), conv(v)%nobs Tim = 0.0 ! ! O Lev = 2 identifica que havera niveis na vertical ! Lev = 2 Flag = 1 DO I=1,conv(v)%nobs write(StationID,'(I8)')I WRITE(10)& StationID ,& conv(v)%lat(I),& conv(v)%lon(I),& Tim, & Lev, & Flag ! print*, "LAT: ", conv(v)%lat(I) ! print*, "LON: ", conv(v)%lon(I) ! ! A primeira variavel escrita corresponde ao nivel da observacao e não conta no ! ctl. ! WRITE(10)& conv(v)%lev (I),& REAL(I) ,& conv(v)%kx (I),& conv(v)%lat (I),& conv(v)%lon (I),& conv(v)%lev (I),& conv(v)%prs (I),& conv(v)%time (I),& conv(v)%pbqc (I),& conv(v)%iuse (I),& conv(v)%iusev(I),& conv(v)%robs (I),& conv(v)%rmod (I),& conv(v)%diff (I) ENDDO WRITE(10)' ENDFILE',0.0,0.0,0.0,0,1 ! Lev = 0 ! WRITE(10)& ! conv(v)%ID(I) ,& ! conv(v)%lat(I),& ! conv(v)%lon(I),& ! Tim, & ! Lev, & ! Flag Close(10) FName = TRIM(FDatSSO) call replace_(FName, "%obstype", trim(adjustl(vname(v)))) call replace_(FName, "%type", TRIM(DataType)) call str_template(FName,nymd,nhms) MapFName = TRIM(FMAP) call replace_(MapFName, "%obstype", trim(adjustl(vname(v)))) call replace_(MapFName, "%type", TRIM(DataType)) call str_template(MapFName,nymd,nhms) ctlFName = TRIm(DirOu)//'/'//TRIM(FCTL) call replace_(CtlFName, "%obstype", trim(adjustl(vname(v)))) call replace_(CtlFName, "%type", TRIM(DataType)) call str_template(CtlFName,nymd,nhms) tmp='%h2Z%d2%MC%y4' call str_template(tmp,nymd,nhms) call writectl(FName,ctlFname,MapFName,tmp) ! A execucao do stnmap foi incluida no shell script run_stat.sh !write(tmp,'(2A)')'stnmap -i ',TRIM(ctlFName) !call system(tmp) ENDDO !dset ps_spatial_summary_obs_01.bin !DTYPE station !STNMAP psspatial_summary_obs_01.map !UNDEF -999.0 !TITLE Station Data Sample !TDEF 1 linear 12z18jan1992 12hr !VARS 12 !ID 1 99 Station ID !kx 1 99 observation type !rlat 1 99 observation latitude (degrees) !rlon 1 99 observation longitude (degrees) !rlev 1 99 observation level reference !prs 1 99 observation pressure (hPa) !rtime 1 99 obs time (minutes relative to analysis time) !pbqc 1 99 input prepbufr qc or event mark !iuse 1 99 analysis usage flag (1=use, -1=monitoring ) !iusev 1 99 analysis usage flag ( value ) !robs 1 99 observation !rmod 1 99 model !diff 1 99 OMF !ENDVARS ! ! Spatial Output ! allocate(rtmp(nlev)) DO v=1,nvars ! WRITE(FName,'(A,A,A)')'spatial_summary_obs_count_',trim(adjustl(vname(v))),'.dat' ! Open(Unit=10, File=Trim(FName)) FName = TRIm(DirOu)//'/'//TRIM(FTxtSSO) call replace_(FName, "%obstype", trim(adjustl(vname(v)))) call replace_(FName, "%type", TRIM(DataType)) call str_template(FName,nymd,nhms) Open(Unit=10,File=trim(FName)) DO I=1,conv(v)%nobs ! if(conv(v)%iuse(I).ne.1)cycle WRITE(10,'(13(F15.6,2x))')& conv(v)%kx (I),& conv(v)%iuse (I),& conv(v)%robs (I),& conv(v)%lev (I),& conv(v)%lat (I),& conv(v)%lon (I),& conv(v)%diff (I),& conv(v)%diffpt (I),& conv(v)%rmod (I),& conv(v)%kmelev (I),& conv(v)%incref (I),& conv(v)%hobl (I),& conv(v)%normdiff (I) ! conv(v)%errf (I),& !conv(v)%time (I),& ENDDO Close(10) ! FName = TRIm(DirOu)//'/KX_%obstype_%type.txt' ! call replace_(FName, "%obstype", trim(adjustl(vname(v)))) ! call replace_(FName, "%type", TRIM(DataType)) ! call str_template(FName,nymd,nhms) ! Open(Unit=10,File=trim(FName),form='formatted') ! DO I=1,conv(v)%nobs ! if(conv(v)%iuse(I).ne.1)cycle ! WRITE(10,'(F9.3)')conv(v)%kx (I) ! ENDDO ! Close(10) ! ! salvando valor médio sobre o globo de OMF/OMF para os dados assimilados ! DO I=1,nkx Flag = count(conv(v)%kx.eq.kx(I).and.conv(v)%iuse.eq.1) if (Flag.eq.0)cycle Open(Unit=11,File=trim(FName),position='append') DO J=1,nlev rtmp(J) = udef Flag = count(conv(v)%lev.eq.levs(J)) if (Flag.eq.0)cycle rtmp(J) = sum(conv(v)%diff,& mask=conv(v)%iuse.eq.1.and.conv(v)%lev.eq.levs(J).and.nint(conv(v)%kx).eq.kx(I))/& count(conv(v)%iuse.eq.1.and.conv(v)%lev.eq.levs(J).and.nint(conv(v)%kx).eq.kx(I)) if (isnan(rtmp(J))) rtmp(J) = udef ! funcao do gfortran e ifort (fortran padrao) ENDDO FName = TRIM(FTxtOMF) call replace_(FName, "%obstype", trim(adjustl(vname(v)))) call replace_(FName, "%type", TRIM(DataType)) call replace_(FName, "%RepType", int2str(kx(I),'(I4.4)')) call str_template(FName,nymd,nhms) Open(Unit=11,File=trim(FName),position='append') write(tmp,'(A,I3.2,A)')'(A,1x',nlev,'F9.3)' WRITE(11,tmp)Trim(TimeFile),(rtmp(J),j=1,nlev) close(11) ENDDO ! close(11) ! ! salvando valor médio sobre o globo de OMF/OMF dos dados monitorados ! Apenas para os dados 153 GPSIPW e 126 RASSDA vento sensor acustico adicionar outros ! DO I=1,nkx Flag = count((((kx(I).eq.153).and.vname(v).eq.' pw').or.((kx(I).eq.126).and.vname(v).eq.' t')).and.conv(v)%iuse.eq.-1) if (Flag.eq.0)cycle Open(Unit=11,File=trim(FName),position='append') DO J=1,nlev rtmp(J) = udef Flag = count(conv(v)%lev.eq.levs(J)) if (Flag.eq.0)cycle rtmp(J) = sum(conv(v)%diff,& mask=conv(v)%iuse.eq.-1.and.conv(v)%lev.eq.levs(J).and.nint(conv(v)%kx).eq.kx(I))/& count(conv(v)%iuse.eq.-1.and.conv(v)%lev.eq.levs(J).and.nint(conv(v)%kx).eq.kx(I)) if (isnan(rtmp(J))) rtmp(J) = udef ! funcao do gfortran e ifort (fortran padrao) ENDDO FName = TRIM(FTxtOMF) call replace_(FName, "%obstype", "%obstype_mon") call replace_(FName, "%obstype", trim(adjustl(vname(v)))) call replace_(FName, "%type", TRIM(DataType)) call replace_(FName, "%RepType", int2str(kx(I),'(I4.4)')) call str_template(FName,nymd,nhms) Open(Unit=11,File=trim(FName),position='append') write(tmp,'(A,I3.2,A)')'(A,1x',nlev,'F9.3)' WRITE(11,tmp)Trim(TimeFile),(rtmp(J),j=1,nlev) close(11) ENDDO ENDDO ! ! Count Summary ! FName = TRIm(DirOu)//'/'//TRIM(FSOC) call replace_(FName, "%type", TRIM(DataType)) call str_template(FName,nymd,nhms) Open(Unit=10, File=trim(FName)) WRITE(10,'( 6A13)')'Obs Name','# of OBS','Used','Not Used','Rejeited','Monitored' DO I=1,nvars WRITE(10,'(A12,5I12)') vname(i), conv(i)%nobs, conv(i)%use(0), conv(i)%nuse(0), conv(i)%rej(0), conv(i)%mon(0) ENDDO Close(10) ! ! used ! FName = TRIm(DirOu)//'/'//TRIM(FSOCL) call replace_(FName, "%stat", "tused") call replace_(FName, "%type", TRIM(DataType)) call str_template(FName,nymd,nhms) Open(Unit=11, File=trim(FName)) write(rowfmt,'(A,I4,A)') '(A9,',nvars,'A9)' write(11,rowfmt)'Levels',vname(1:nvars) write(rowfmt,'(A,I4,A)') '(F9.1,',nvars,'I9)' Do i=1,nlev write(11,rowfmt)levs(i),(conv(j)%use(i),j=1,nvars) Enddo ! ! Not Used ! FName = TRIm(DirOu)//'/'//TRIM(FSOCL) call replace_(FName, "%stat", "nused") call replace_(FName, "%type", TRIM(DataType)) call str_template(FName,nymd,nhms) Open(Unit=11, File=trim(FName)) write(rowfmt,'(A,I4,A)') '(A9,',nvars,'A9)' write(11,rowfmt)'Levels',vname(1:nvars) write(rowfmt,'(A,I4,A)') '(F9.1,',nvars,'I9)' Do i=1,nlev write(11,rowfmt)levs(i),(conv(j)%nuse(i),j=1,nvars) Enddo ! ! Rejeited ! FName = TRIm(DirOu)//'/'//TRIM(FSOCL) call replace_(FName, "%stat", "rejec") call replace_(FName, "%type", TRIM(DataType)) call str_template(FName,nymd,nhms) Open(Unit=11, File=trim(FName)) write(rowfmt,'(A,I4,A)') '(A9,',nvars,'A9)' write(11,rowfmt)'Levels',vname(1:nvars) write(rowfmt,'(A,I4,A)') '(F9.1,',nvars,'I9)' Do i=1,nlev write(11,rowfmt)levs(i),(conv(j)%rej(i),j=1,nvars) Enddo ! ! Monitored ! FName = TRIm(DirOu)//'/'//TRIM(FSOCL) call replace_(FName, "%stat", "monit") call replace_(FName, "%type", TRIM(DataType)) call str_template(FName,nymd,nhms) Open(Unit=11, File=trim(FName)) write(rowfmt,'(A,I4,A)') '(A9,',nvars,'A9)' write(11,rowfmt)'Levels',vname(1:nvars) write(rowfmt,'(A,I4,A)') '(F9.1,',nvars,'I9)' Do i=1,nlev write(11,rowfmt)levs(i),(conv(j)%mon(i),j=1,nvars) Enddo ! ! Vies and RMS Summary ! ! ! HN ! ! Open(Unit=10, File=trim(dirou)//'/RMSE_HN_'//bufr//'.txt') ! write(rowfmt,'(A,I4,A)') '(A9,',nlev+1,'F15.9)' ! DO I=1,nvars ! WRITE(10,rowfmt) vname(i), (conv(i)%rmse(J,1),J=0,nlev) ! ENDDO ! Close(10) ! Open(Unit=10, File=trim(dirou)//'/VIES_HN_'//bufr//'.txt') ! write(rowfmt,'(A,I4,A)') '(A9,',nlev+1,'F15.9)' ! DO I=1,nvars ! WRITE(10,rowfmt) vname(i), (conv(i)%vies(J,1),J=0,nlev) ! ENDDO ! Close(10) ! ! HS ! ! Open(Unit=10, File=trim(dirou)//'/RMSE_HS_'//bufr//'.txt') ! write(rowfmt,'(A,I4,A)') '(A9,',nlev+1,'F15.9)' ! DO I=1,nvars ! WRITE(10,rowfmt) vname(i), (conv(i)%rmse(J,2),J=0,nlev) ! ENDDO ! Close(10) ! Open(Unit=10, File=trim(dirou)//'/VIES_HS_'//bufr//'.txt') ! write(rowfmt,'(A,I4,A)') '(A9,',nlev+1,'F15.9)' ! DO I=1,nvars ! WRITE(10,rowfmt) vname(i), (conv(i)%vies(J,2),J=0,nlev) ! ENDDO ! Close(10) #ifdef SUMMARY ! ! Writing summary of count ! WRITE(6,'(1x,A75)')'---------------------------------------------------------------------------' WRITE(6,'(1x,A75)')' S U M M A R Y ' WRITE(6,'(1x,A75)')'---------------------------------------------------------------------------' WRITE(6,'( 6A12)')'Obs Name','# of OBS','Used','Not Used','Rejected','Monitored' WRITE(6,'(1X,A75)')'---------------------------------------------------------------------------' DO I=1,nvars WRITE(6,'(A12,5I12)') vname(i), conv(i)%nobs, conv(i)%use(0), conv(i)%nuse(0), conv(i)%rej(0), conv(i)%mon(0) ENDDO WRITE(6,'(1X,A75)')'---------------------------------------------------------------------------' WRITE(6,'(A12,5I12)')'Total',sum(conv(:)%nobs),sum(conv(:)%use(0)),sum(conv(:)%nuse(0)), sum(conv(:)%rej(0)), sum(conv(:)%mon(0)) WRITE(6,'(1X,A75)')'---------------------------------------------------------------------------' write(rowfmt,'(A,I4,A)') '(A6,',nlev,'F9.1)' WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------' write(6,'(1x,A120)')' U S E D B Y L E V E L ' write(6,rowfmt)'Top',(levs(i),i=1,nlev) write(6,rowfmt)'Bot',1200.0,(levs(i)-0.1,i=1,nlev-1) WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------' write(rowfmt,'(A,I4,A)') '(A6,',nlev,'I9)' do i=1,nvars write(6,rowfmt)vname(i), conv(i)%use(1:nlev) enddo write(rowfmt,'(A,I4,A)') '(A6,',nlev,'F9.1)' WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------' write(6,'(1x,A120)')' N O T U S E D B Y L E V E L ' write(6,rowfmt)'Top',(levs(i),i=1,nlev) write(6,rowfmt)'Bot',1200.0,(levs(i)-0.1,i=1,nlev-1) WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------' write(rowfmt,'(A,I4,A)') '(A6,',nlev,'I9)' do i=1,nvars write(6,rowfmt)vname(i), conv(i)%nuse(1:nlev) enddo write(rowfmt,'(A,I4,A)') '(A6,',nlev,'F9.1)' WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------' write(6,'(1x,A120)')' R E J E I T E D B Y L E V E L ' write(6,rowfmt)'Top',(levs(i),i=1,nlev) write(6,rowfmt)'Bot',1200.0,(levs(i)-0.1,i=1,nlev-1) WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------' write(rowfmt,'(A,I4,A)') '(A6,',nlev,'I9)' do i=1,nvars write(6,rowfmt)vname(i), conv(i)%rej(1:nlev) enddo write(rowfmt,'(A,I4,A)') '(A6,',nlev,'F9.1)' WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------' write(6,'(1x,A120)')' M O N I T O R E D B Y L E V E L ' write(6,rowfmt)'Top',(levs(i),i=1,nlev) write(6,rowfmt)'Bot',1200.0,(levs(i)-0.1,i=1,nlev-1) WRITE(6,'(1X,A120)')'------------------------------------------------------------------------------------------------------------------------' write(rowfmt,'(A,I4,A)') '(A6,',nlev,'I9)' do i=1,nvars write(6,rowfmt)vname(i), conv(i)%mon(1:nlev) enddo #endif STOP 'ENDS NORMALLY' 997 PRINT *,'error read in diag file 997' stop 1 998 PRINT *,'error read in diag file 998' stop 2 999 PRINT *,'error read in diag file 999' stop 3 END PROGRAM diag_GSI_conv ! !EOCdiag_gsi_conv.f90 !---------------------------------------------------------------------