var/convertor/ens_obs_var.f90

   1
   2 !xlf90 -g -qrealsize=8 ens_obs_var.f90 -o ens_obs_var.exe
   3
   4 program ens_obs_var
   5
   6      parameter (ne=4)
   7      character(len=40) :: platform, date, name, id
   8
   9      CHARACTER (LEN = 80)          :: filename1, filename2
  10      CHARACTER (LEN = 120)         :: fmt_info, &
  11                                       fmt_srfc, &
  12                                       fmt_srfc2, &
  13                                       fmt_each, fmt_each_iv
  14      character (len = 160)         :: header_string, eheader_string
  15
  16      integer                       :: n, iunit(ne),ens_size
  17      real                          :: data, error
  18      real                          :: missing_v, missing_qc, missing_e
  19      integer                       :: iqc, iost
  20      real                :: xlon,xlat
  21      real                :: data1,error1, data2,error2, &
  22                             data3,error3, data4,error4, data5,error5, &
  23                             data6,error6, data7,error7, data8,error8, &
  24                             data9,error9
  25      integer             :: iqc1, iqc2, iqc3, iqc4, iqc5, iqc6, iqc7, iqc8, iqc9
  26      real, dimension(ne) :: edata1,eerror1, edata2,eerror2, &
  27                             edata3,eerror3, edata4,eerror4, edata5,eerror5, &
  28                             edata6,eerror6, edata7,eerror7, edata8,eerror8, &
  29                             edata9,eerror9, euinv,evinv,etinv,eqinv
  30      integer,dimension(ne):: eiqc1, eiqc2, eiqc3, eiqc4, eiqc5, eiqc6, eiqc7, eiqc8, eiqc9
  31      real                :: utotsp, vtotsp, ttotsp, qtotsp
  32
  33      fmt_info = '(A12,1X,A19,1X,A40,1X,I6,3(F12.3,11X),6X,A5)'
  34      fmt_srfc = '(F12.3,I4,F7.2,F12.3,I4,F7.3)'
  35      fmt_srfc2 = '(7(:,f12.3,i4,f7.2))'
  36      fmt_each = '(3(F12.3,I4,F7.2),11X,3(F12.3,I4,F7.2),11X,1(F12.3,I4,F7.2)))'
  37      fmt_each_iv = '(3(F12.3,I4,F7.2),11X,3(F12.3,I4,F7.2),11X,1(F12.3,I4,F7.2),4f17.7)'
  38
  39      missing_v  = -888888.000
  40      missing_qc = -88.0
  41      missing_e  = -888.0
  42
  43      open(12,file='filtered_obs',status='old')
  44      open(13,file='filtered_obs.var')
  45
  46      do n=1,ne
  47         iunit(n)=100+n
  48         write(filename1,'(a,i2.2)') 'filtered_obs.e',n
  49         open(iunit(n),file=filename1,status='old')
  50      enddo
  51
  52      do
  53        read(unit=12,fmt='(a)',iostat=iost) header_string
  54        do n=1,ne
  55           read(unit=iunit(n),fmt='(a)',iostat=iost) eheader_string
  56        enddo
  57        if (iost /= 0) then
  58           print*, 'Problem reading obs header, exit'
  59           stop
  60        endif
  61        write(unit=13,fmt='(a)') header_string
  62        if (header_string(1:6) == '#-----' .or. header_string(1:21) == 'MODIFIED FILTERED OBS') exit
  63      enddo
  64
  65      istamax=250000
  66
  67      do ista=1,istamax
  68         read(12,fmt=trim(fmt_info),end=999) platform,date,name,ilev,xlat,xlon,elevation,id
  69         write(13,fmt=trim(fmt_info))        platform,date,name,ilev,xlat,xlon,elevation,id
  70
  71         ! print*,'ista=',ista
  72         ! print*,platform,date,name,ilev,xlat,xlon,elevation,id
  73         read(12,fmt=trim(fmt_srfc),end=999) data1,iqc1,error1,data2,iqc2,error2
  74         write(13,fmt=trim(fmt_srfc2))        data1,iqc1,error1,data2,iqc2,error2
  75
  76         do n=1,ne
  77            read(iunit(n),fmt=trim(fmt_info),end=999)platform,date,name,ilev,xlat,xlon,elevation,id
  78            read(iunit(n),fmt=trim(fmt_srfc),end=999)edata1(n),eiqc1(n),eerror1(n),edata2(n),eiqc2(n),eerror2(n)
  79            if(edata1(n).ne.data1 .or. eiqc1(n).ne.iqc1 .or. eerror1(n).ne.error1) then
  80              print*,'Observation mismatch!, ensemble#',n
  81              print*,platform,date,name,ilev,xlat,xlon,elevation,id
  82              print*,data1,iqc1,error1,data2,iqc2,error2
  83              print*,edata1(n),eiqc1(n),eerror1(n),edata2(n),eiqc2(n),eerror2(n)
  84              stop
  85            endif
  86         enddo
  87
  88         do k=1,ilev
  89            read(12,fmt=trim(fmt_each),end=999) data3,iqc3,error3, data4,iqc4,error4, data5,iqc5,error5, &
  90                                                data6,iqc6,error6, data7,iqc7,error7, data8,iqc8,error8, &
  91                                                data9,iqc9,error9
  92            do n=1,ne
  93               read(iunit(n),fmt=trim(fmt_each_iv),end=999)   &
  94                   edata3(n),eiqc3(n),eerror3(n), edata4(n),eiqc4(n),eerror4(n), edata5(n),eiqc5(n),eerror5(n), &
  95                   edata6(n),eiqc6(n),eerror6(n), edata7(n),eiqc7(n),eerror7(n), edata8(n),eiqc8(n),eerror8(n), &
  96                   edata9(n),eiqc9(n),eerror9(n), euinv(n),evinv(n),etinv(n),eqinv(n)
  97            enddo
  98            call compute_total_spread(edata4,eiqc4,eerror4,euinv,ne,utotsp)
  99            call compute_total_spread(edata5,eiqc5,eerror5,evinv,ne,vtotsp)
 100            call compute_total_spread(edata7,eiqc7,eerror7,etinv,ne,ttotsp)
 101            call compute_total_spread(edata9,eiqc9,eerror9,eqinv,ne,qtotsp)
 102
 103            if(qtotsp .ge. 0) qtotsp=max(0.01,qtotsp)
 104            write(13,fmt=trim(fmt_each))data3,iqc3,error3, data4,iqc4,utotsp, data5,iqc5,vtotsp, &
 105                                        data6,iqc6,error6, data7,iqc7,ttotsp, data8,iqc8,error8, &
 106                                        data9,iqc9,qtotsp
 107         enddo
 108      enddo
 109 999  print*, ista-1
 110      if(ista-1.eq.istamax)then
 111      print*,'you should give larger istamax'
 112
 113      endif
 114
 115 contains
 116
 117 subroutine compute_total_spread(edata,eqc,eerr,einv,esize,etotsp)
 118     integer :: esize, i
 119     real    :: edata(esize), eerr(esize), einv(esize), etotsp
 120     integer :: eqc(esize)
 121     real    :: emean, evar
 122
 123     do i=1,esize
 124        if (edata(i)==missing_r .or. eqc(i)==missing_qc .or. eerr(i)==missing_e) then
 125           etotsp=missing_e
 126           return
 127        endif
 128        if (eerr(i).ne.eerr(1)) then
 129           print*,'Ensemble obs error mismatch: ',eerr(1), eerr(i)
 130           stop
 131        endif
 132     enddo
 133     call compute_mean_var(einv,esize,emean,evar)
 134     etotsp=sqrt(eerr(1)**2+evar)
 135     return
 136 end subroutine compute_total_spread
 137
 138 subroutine compute_mean_var(evalue,esize,emean,evar)
 139     integer :: esize
 140     real    :: evalue(esize), emean, evar
 141     emean=sum(evalue)/esize
 142     evar=sum((evalue-emean)**2)/(esize-1)
 143 end subroutine compute_mean_var
 144
 145 end program
 146
 147