var/da/da_radar/da_write_oa_radar_ascii.inc

   1 subroutine da_write_oa_radar_ascii ( ob, iv, re, it )
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: write out OMB and OMA vector structure for radar data.
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type (y_type),     intent(in)  :: ob       ! Observation structure.
  10    type (iv_type),    intent(in)  :: iv       ! O-B structure.
  11    type (y_type),     intent(in)  :: re       ! O-A structure.
  12    integer,           intent(in)  :: it       ! external iteration counter
  13
  14    character(len=filename_len)               :: filename , outfname
  15    character(len=filename_len), allocatable  :: filename1(:)
  16    integer :: num_obs, num, levels, ios
  17    integer :: n, i, k, m, kk, l
  18    integer :: oma_radar_unit, omb_radar_unit, omb_radar_iter_unit
  19    real    :: lat, lon, height, zk
  20    real    :: rv_obs, rv_inv, rv_error, rv_inc
  21    real    :: rf_obs, rf_inv, rf_error, rf_inc
  22    real    :: rn_obs, rn_inv, rn_error, rn_inc
  23    real    :: sn_obs, sn_inv, sn_error, sn_inc
  24    real    :: gr_obs, gr_inv, gr_error, gr_inc
  25    real    :: qv_obs, qv_inv, qv_error, qv_inc
  26    integer :: rv_qc, rf_qc
  27    integer :: rn_qc, sn_qc, gr_qc, qv_qc
  28    character(len=5) :: stn_id
  29
  30    if (trace_use) call da_trace_entry("da_write_oa_radar_ascii")
  31
  32    write(unit=message(1),fmt='(A)') 'Writing radar OMA ascii file'
  33    call da_message(message(1:1))
  34
  35    ! myproc is set in var/da/da_main/da_wrfvar_init1.inc
  36    ! it is 0 for serial mode
  37    write(unit=filename, fmt='(a,i2.2,a,i4.4)') 'radar_omb_oma_',it,'.', myproc
  38
  39    call da_get_unit(oma_radar_unit)
  40    open(unit=oma_radar_unit,file=trim(filename),form='formatted',iostat=ios)
  41    if (ios /= 0) Then
  42        call da_error(__FILE__,__LINE__, &
  43          (/"Cannot open oma radar file "//trim(filename)/))
  44    endif
  45
  46    if (iv % info(radar)%nlocal  >0 ) then
  47       num_obs = 0
  48       do n = 1, iv% info(radar)%nlocal
  49          if (iv%info(radar)%proc_domain(1,n)) num_obs=num_obs+1
  50       end do
  51       if (num_obs > 0) then
  52          write(oma_radar_unit,'(a20,i8)')'radar', num_obs
  53          num_obs = 0
  54          do n = 1, iv % info(radar)%nlocal
  55             if (iv%info(radar)%proc_domain(1,n)) then
  56                num_obs = num_obs + 1
  57                write(oma_radar_unit,'(i8)') iv % info(radar) % levels(n)
  58                do k = 1, iv % info(radar) % levels(n)
  59                   if ( .not. use_radar_rhv ) then
  60                      write(oma_radar_unit,'(2i8,a5,2f9.2,f17.7,2(2f17.7,i8,2f17.7),f17.7)')&
  61                         num_obs , k, 'RADAR', &
  62                         iv % info (radar) % lat(1,n), &  ! Latitude
  63                         iv % info (radar) % lon(1,n), &  ! Longitude
  64                         iv % radar(n) % height(k),    &  ! Obs height in m
  65                         ob%radar(n)%rv(k), &
  66                         iv%radar(n)%rv(k)%inv,iv%radar(n)%rv(k)%qc,iv%radar(n)%rv(k)%error,&
  67                         re%radar(n)%rv(k), & ! O, O-B, qc, err, O-A rv
  68                         ob%radar(n)%rf(k), &
  69                         iv%radar(n)%rf(k)%inv,iv%radar(n)%rf(k)%qc,iv%radar(n)%rf(k)%error,&
  70                         re%radar(n)%rf(k),iv%info(radar)%zk(k,n)   ! O, O-B, qc, err, O-A rf
  71                   else ! use_radar_rhv
  72                      if ( use_radar_rqv ) then
  73                         write(oma_radar_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')&
  74                            num_obs , k, 'RADAR', &
  75                            iv % info (radar) % lat(1,n), &  ! Latitude
  76                            iv % info (radar) % lon(1,n), &  ! Longitude
  77                            iv % radar(n) % height(k),    &  ! Obs height in m
  78                            ob%radar(n)%rv(k),   &  ! radial velocity
  79                            iv%radar(n)%rv(k)%inv,iv%radar(n)%rv(k)%qc,iv%radar(n)%rv(k)%error,&
  80                            re%radar(n)%rv(k),   &  ! O, O-B, qc, err, O-A rv
  81                            iv%radar(n)%rrno(k), &  ! retrieved rain water
  82                            iv%radar(n)%rrn(k)%inv,iv%radar(n)%rrn(k)%qc,iv%radar(n)%rrn(k)%error,&
  83                            re%radar(n)%rrn(k),  &  ! O, O-B, qc, err, O-A rrn
  84                            iv%radar(n)%rsno(k), &  ! retrieved snow
  85                            iv%radar(n)%rsn(k)%inv,iv%radar(n)%rsn(k)%qc,iv%radar(n)%rsn(k)%error,&
  86                            re%radar(n)%rsn(k),  &  ! O, O-B, qc, err, O-A rsn
  87                            iv%radar(n)%rgro(k), &  ! retrieved graupel
  88                            iv%radar(n)%rgr(k)%inv,iv%radar(n)%rgr(k)%qc,iv%radar(n)%rgr(k)%error,&
  89                            re%radar(n)%rgr(k),  &  ! O, O-B, qc, err, O-A rgr
  90                            iv%radar(n)%rqvo(k), &  ! retrieved water vapor
  91                            iv%radar(n)%rqv(k)%inv,iv%radar(n)%rqv(k)%qc,iv%radar(n)%rqv(k)%error,&
  92                            re%radar(n)%rqv(k)      ! O, O-B, qc, err, O-A rqv
  93                      else ! not use_radar_rqv
  94                         write(oma_radar_unit,'(2i8,a5,2f9.2,f17.7,4(2f17.7,i8,2f17.7))')&
  95                            num_obs , k, 'RADAR', &
  96                            iv % info (radar) % lat(1,n), &  ! Latitude
  97                            iv % info (radar) % lon(1,n), &  ! Longitude
  98                            iv % radar(n) % height(k),    &  ! Obs height in m
  99                            ob%radar(n)%rv(k),   &  ! radial velocity
 100                            iv%radar(n)%rv(k)%inv,iv%radar(n)%rv(k)%qc,iv%radar(n)%rv(k)%error,&
 101                            re%radar(n)%rv(k),   &  ! O, O-B, qc, err, O-A rv
 102                            iv%radar(n)%rrno(k), &  ! retrieved rain water
 103                            iv%radar(n)%rrn(k)%inv,iv%radar(n)%rrn(k)%qc,iv%radar(n)%rrn(k)%error,&
 104                            re%radar(n)%rrn(k),  &  ! O, O-B, qc, err, O-A rrn
 105                            iv%radar(n)%rsno(k), &  ! retrieved snow
 106                            iv%radar(n)%rsn(k)%inv,iv%radar(n)%rsn(k)%qc,iv%radar(n)%rsn(k)%error,&
 107                            re%radar(n)%rsn(k),  &  ! O, O-B, qc, err, O-A rsn
 108                            iv%radar(n)%rgro(k), &  ! retrieved graupel
 109                            iv%radar(n)%rgr(k)%inv,iv%radar(n)%rgr(k)%qc,iv%radar(n)%rgr(k)%error,&
 110                            re%radar(n)%rgr(k)      ! O, O-B, qc, err, O-A rgr
 111                      end if ! check use_radar_rqv
 112                   end if ! check use_radar_rhv
 113                end do ! leve loop
 114             end if ! check proc_domain
 115          end do ! nlocal loop
 116       end if ! check num_obs
 117    end if ! check nlocal
 118
 119 #ifdef DM_PARALLEL
 120    ! Wait to ensure all temporary files have been written
 121    call mpi_barrier(comm, ierr)
 122 #endif
 123
 124    close (oma_radar_unit)
 125    call da_free_unit(oma_radar_unit)
 126
 127    IF (rootproc) THEN
 128       call da_get_unit(omb_radar_unit)
 129       allocate (filename1(0:num_procs-1))
 130       do k = 0,num_procs-1
 131          write(unit=filename1(k),fmt ='(a,i2.2,a,i4.4)')'radar_omb_oma_',it,'.',k
 132       end do
 133       call da_get_unit(omb_radar_iter_unit)
 134       write(unit=outfname,fmt ='(a,i2.2)')'radar_omb_oma_',it
 135       open(omb_radar_iter_unit,file=trim(outfname),form='formatted', status='replace', iostat=ios)
 136       if (ios /= 0) call da_error(__FILE__,__LINE__, &
 137           (/"Cannot open file "//trim(outfname)/))
 138    ENDIF
 139
 140    num_obs = 0
 141    IF (iv % info(radar)%nlocal  >0 ) then
 142       do n = 1, iv% info(radar)%nlocal
 143          if (iv%info(radar)%proc_domain(1,n)) num_obs=num_obs+1
 144       end do
 145    ENDIF
 146    call da_proc_sum_int(num_obs)
 147    IF (num_obs > 0 .and. rootproc) then
 148       write(omb_radar_iter_unit,'(a20,i8)')'radar', num_obs
 149       num_obs = 0
 150       num = 0
 151       do k = 0,num_procs-1
 152          open(omb_radar_unit,file=trim(filename1(k)),status='old',iostat=ios)
 153          if (ios /= 0) call da_error(__FILE__,__LINE__, &
 154                                        (/"Cannot open file "//trim(filename1(k))/))
 155          read(omb_radar_unit, '(20x,i8)', iostat=ios)num_obs
 156          IF(ios /= 0)THEN
 157             write(unit=message(1),fmt='(A,A)') 'Nothing to read from ',filename1(k)
 158             call da_message(message(1:1))
 159             cycle
 160          ENDIF
 161          if (num_obs > 0) then
 162             do n = 1, num_obs
 163                read(omb_radar_unit,'(i8)') levels
 164                write(omb_radar_iter_unit,'(i8)') levels
 165                num = num + 1
 166                do m = 1, levels
 167                   if ( .not. use_radar_rhv ) then
 168                      read(omb_radar_unit,'(2i8,a5,2f9.2,f17.7,2(2f17.7,i8,2f17.7),f17.7)')&
 169                          kk,l, stn_id,     &  ! Station
 170                          lat, lon, height, &  ! Lat/lon, height
 171                          rv_obs, rv_inv, rv_qc, rv_error, rv_inc, &
 172                          rf_obs, rf_inv, rf_qc, rf_error, rf_inc, zk
 173                      write(omb_radar_iter_unit,'(2i8,a5,2f9.2,f17.7,2(2f17.7,i8,2f17.7),f17.7)')&
 174                          num,m,stn_id,     &  ! Station
 175                          lat, lon, height, &  ! Lat/lon, height
 176                          rv_obs, rv_inv, rv_qc, rv_error, rv_inc, &
 177                          rf_obs, rf_inv, rf_qc, rf_error, rf_inc, zk
 178                   else ! use_radar_rhv
 179                      if ( use_radar_rqv ) then
 180                         read(omb_radar_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')&
 181                             kk,l, stn_id,     &  ! Station
 182                             lat, lon, height, &  ! Lat/lon, height
 183                             rv_obs, rv_inv, rv_qc, rv_error, rv_inc, &
 184                             rn_obs, rn_inv, rn_qc, rn_error, rn_inc, &
 185                             sn_obs, sn_inv, sn_qc, sn_error, sn_inc, &
 186                             gr_obs, gr_inv, gr_qc, gr_error, gr_inc, &
 187                             qv_obs, qv_inv, qv_qc, qv_error, qv_inc
 188                         write(omb_radar_iter_unit,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))')&
 189                             num,m,stn_id,     &  ! Station
 190                             lat, lon, height, &  ! Lat/lon, height
 191                             rv_obs, rv_inv, rv_qc, rv_error, rv_inc, &
 192                             rn_obs, rn_inv, rn_qc, rn_error, rn_inc, &
 193                             sn_obs, sn_inv, sn_qc, sn_error, sn_inc, &
 194                             gr_obs, gr_inv, gr_qc, gr_error, gr_inc, &
 195                             qv_obs, qv_inv, qv_qc, qv_error, qv_inc
 196                      else ! not use_radar_rqv
 197                         read(omb_radar_unit,'(2i8,a5,2f9.2,f17.7,4(2f17.7,i8,2f17.7))')&
 198                             kk,l, stn_id,     &  ! Station
 199                             lat, lon, height, &  ! Lat/lon, height
 200                             rv_obs, rv_inv, rv_qc, rv_error, rv_inc, &
 201                             rn_obs, rn_inv, rn_qc, rn_error, rn_inc, &
 202                             sn_obs, sn_inv, sn_qc, sn_error, sn_inc, &
 203                             gr_obs, gr_inv, gr_qc, gr_error, gr_inc
 204                         write(omb_radar_iter_unit,'(2i8,a5,2f9.2,f17.7,4(2f17.7,i8,2f17.7))')&
 205                             num,m,stn_id,     &  ! Station
 206                             lat, lon, height, &  ! Lat/lon, height
 207                             rv_obs, rv_inv, rv_qc, rv_error, rv_inc, &
 208                             rn_obs, rn_inv, rn_qc, rn_error, rn_inc, &
 209                             sn_obs, sn_inv, sn_qc, sn_error, sn_inc, &
 210                             gr_obs, gr_inv, gr_qc, gr_error, gr_inc
 211                      end if ! check use_radar_rqv
 212                   end if ! check use_radar_rhv
 213                end do ! level loop
 214             end do ! num_obs loop
 215          end if ! check num_obs
 216       end do ! num_proc loop
 217    END IF
 218
 219    if (rootproc) then
 220       close(omb_radar_unit)
 221       close(omb_radar_iter_unit)
 222       call da_free_unit(omb_radar_unit)
 223       call da_free_unit(omb_radar_iter_unit)
 224       deallocate (filename1)
 225    end if
 226
 227    if (trace_use) call da_trace_exit("da_write_oa_radar_ascii")
 228
 229 end subroutine da_write_oa_radar_ascii