wrfv2_fire/var/da/da_radiance/da_write_iv_rad_ascii.inc

   1 subroutine da_write_iv_rad_ascii (it, ob, iv )
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: write out innovation vector structure for radiance data.
   5    !---------------------------------------------------------------------------
   6
   7    implicit none
   8
   9    integer      ,     intent(in)  :: it       ! outer loop count
  10    type (y_type),     intent(in)  :: ob       ! Observation structure.
  11    type (iv_type),    intent(in)  :: iv       ! O-B structure.
  12
  13    integer                        :: n        ! Loop counter.
  14    integer                        :: i, k, l  ! Index dimension.
  15    integer                        :: nlevelss ! Number of obs levels.
  16
  17    integer            :: ios, innov_rad_unit
  18    character(len=filename_len)  :: filename
  19    character(len=7)   :: surftype
  20    integer            :: ndomain
  21
  22    if (trace_use) call da_trace_entry("da_write_iv_rad_ascii")
  23
  24    do i = 1, iv%num_inst
  25       if (iv%instid(i)%num_rad < 1) cycle
  26
  27       ! count number of obs within the loc%proc_domain
  28       ! ---------------------------------------------
  29       ndomain = 0
  30       do n =1,iv%instid(i)%num_rad
  31          if (iv%instid(i)%info%proc_domain(1,n)) then
  32             ndomain = ndomain + 1
  33          end if
  34       end do
  35       if (ndomain < 1) cycle
  36
  37       write(unit=filename, fmt='(i2.2,a,i4.4)') it,'_inv_'//trim(iv%instid(i)%rttovid_string)//'.', myproc
  38
  39       call da_get_unit(innov_rad_unit)
  40       open(unit=innov_rad_unit,file=trim(filename),form='formatted',iostat=ios)
  41       if (ios /= 0 ) then
  42          call da_error(__FILE__,__LINE__, &
  43             (/"Cannot open innovation radiance file"//filename/))
  44       Endif
  45       write(unit=innov_rad_unit,fmt='(a,a,i7,a,i5,a)') trim(iv%instid(i)%rttovid_string), &
  46                         ' number-of-pixels : ', ndomain, &
  47                         ' channel-number-of-each-pixel : ', iv%instid(i)%nchan, &
  48                         ' index-of-channels : '
  49       write(unit=innov_rad_unit,fmt='(10i5)') iv%instid(i)%ichan
  50
  51       write(unit=innov_rad_unit,fmt=*) ' pixel-info : i date scanpos landsea_mask  elv lat lon  satzen satazi'
  52       write(unit=innov_rad_unit,fmt=*) ' grid%xb-surf-info : i t2m mr2m(ppmv) u10 v10 ps ts smois tslb snowh isflg &
  53                     & soiltyp vegtyp vegfra elev clwp'
  54       ndomain = 0
  55       do n =1,iv%instid(i)%num_rad
  56          if (iv%instid(i)%info%proc_domain(1,n)) then
  57             ndomain=ndomain+1
  58             write(unit=innov_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,4f8.2)') 'INFO : ', ndomain, &
  59                                 iv%instid(i)%info%date_char(n), &
  60                                 iv%instid(i)%scanpos(n),   &
  61                                 iv%instid(i)%landsea_mask(n), &
  62                                 iv%instid(i)%info%elv(n),  &
  63                                 iv%instid(i)%info%lat(1,n),  &
  64                                 iv%instid(i)%info%lon(1,n), &
  65                                 iv%instid(i)%satzen(n),    &
  66                                 iv%instid(i)%satazi(n)
  67             select case (iv%instid(i)%isflg(n))
  68             case (0) ;
  69                surftype = ' SEA : '
  70             case (1) ;
  71                surftype = ' ICE : '
  72             case (2) ;
  73                surftype = 'LAND : '
  74             case (3) ;
  75                surftype = 'SNOW : '
  76             case (4) ;
  77                surftype = 'MSEA : '
  78             case (5) ;
  79                surftype = 'MICE : '
  80             case (6) ;
  81                surftype = 'MLND : '
  82             case (7) ;
  83                surftype = 'MSNO : '
  84             end select
  85             write(unit=innov_rad_unit,fmt='(a,i7,9f10.2,3i3,f8.3,f10.2,f8.3)') surftype, n, &
  86                              iv%instid(i)%t2m(n), &
  87                              iv%instid(i)%mr2m(n),   &
  88                              iv%instid(i)%u10(n), &
  89                              iv%instid(i)%v10(n),  &
  90                              iv%instid(i)%ps(n),  &
  91                              iv%instid(i)%ts(n),  &
  92                              iv%instid(i)%smois(n),  &
  93                              iv%instid(i)%tslb(n),  &
  94                              iv%instid(i)%snowh(n), &
  95                              iv%instid(i)%isflg(n), &
  96                              nint(iv%instid(i)%soiltyp(n)), &
  97                              nint(iv%instid(i)%vegtyp(n)), &
  98                              iv%instid(i)%vegfra(n), &
  99                              iv%instid(i)%elevation(n), &
 100                              iv%instid(i)%clwp(n)
 101
 102             write(unit=innov_rad_unit,fmt='(a)') 'OBS  : '
 103             write(unit=innov_rad_unit,fmt='(10f11.2)') ob%instid(i)%tb(:,n)
 104             write(unit=innov_rad_unit,fmt='(a)') 'BAK  : '
 105             write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%tb_xb(:,n)
 106             write(unit=innov_rad_unit,fmt='(a)') 'IVBC : '
 107             write(unit=innov_rad_unit,fmt='(10f11.2)')  iv%instid(i)%tb_inv(:,n)
 108             write(unit=innov_rad_unit,fmt='(a)') 'EMS  : '
 109             write(unit=innov_rad_unit,fmt='(10f11.2)')  iv%instid(i)%emiss(1:iv%instid(i)%nchan,n)
 110             if (rtm_option==rtm_option_crtm .and. write_jacobian) then
 111                write(unit=innov_rad_unit,fmt='(a)') 'EMS_JACOBIAN : '
 112                write(unit=innov_rad_unit,fmt='(10f10.3)') iv%instid(i)%emiss_jacobian(1:iv%instid(i)%nchan,n)
 113             end if
 114             write(unit=innov_rad_unit,fmt='(a)') 'ERR  : '
 115             write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%tb_error(:,n)
 116             write(unit=innov_rad_unit,fmt='(a)') 'QC   : '
 117             write(unit=innov_rad_unit,fmt='(10i11)') iv%instid(i)%tb_qc(:,n)
 118
 119             if (write_profile) then
 120                nlevelss  = iv%instid(i)%nlevels
 121                if ( rtm_option == rtm_option_rttov ) then
 122 #ifdef RTTOV
 123                   ! first, write RTTOV levels
 124                   write(unit=innov_rad_unit,fmt=*) 'RTM_level pres(mb) T(k) Q(ppmv)'
 125                   do k = 1, nlevelss
 126                      write(unit=innov_rad_unit,fmt='(i3,f10.2,f8.2,e11.4)') &
 127                         k, &                             ! RTTOV levels
 128                         coefs(i) % ref_prfl_p(k) , &
 129                         iv%instid(i)%t(k,n) , &
 130                         iv%instid(i)%mr(k,n)
 131                   end do  ! end loop RTTOV level
 132                   ! second, write WRF model levels
 133                   write(unit=innov_rad_unit,fmt=*) &
 134                      'WRF_level pres(mb) T(k) q(g/kg) clw(g/kg) rain(g/kg)'
 135                   do k=kts,kte
 136                      write(unit=innov_rad_unit,fmt='(i3,f10.2,f8.2,3e11.4)') &
 137                         k,  &                     ! WRF model levels
 138                         iv%instid(i)%pm(k,n) , &
 139                         iv%instid(i)%tm(k,n) , &
 140                         iv%instid(i)%qm(k,n)*1000 , &
 141                         iv%instid(i)%qcw(k,n)*1000.0, &
 142                         iv%instid(i)%qrn(k,n)*1000.0
 143                   end do ! end loop model level
 144 #endif
 145                end if ! end if rtm_option_rttov
 146
 147                if ( rtm_option == rtm_option_crtm ) then
 148 #ifdef CRTM
 149                   write(unit=innov_rad_unit,fmt=*) &
 150                      'level fullp(mb) halfp(mb) t(k) q(g/kg) water(mm) ice(mm) rain(mm) snow(mm) graupel(mm) hail(mm)'
 151                   if (crtm_cloud) then
 152                      do k=1,iv%instid(i)%nlevels-1
 153                         write(unit=innov_rad_unit,fmt='(i3,2f10.2,f8.2,13f8.3)') &
 154                            k,  &
 155                            iv%instid(i)%pf(k,n), &
 156                            iv%instid(i)%pm(k,n), &
 157                            iv%instid(i)%tm(k,n), &
 158                            iv%instid(i)%qm(k,n), &
 159                            iv%instid(i)%qcw(k,n), &
 160                            iv%instid(i)%qci(k,n), &
 161                            iv%instid(i)%qrn(k,n), &
 162                            iv%instid(i)%qsn(k,n), &
 163                            iv%instid(i)%qgr(k,n), &
 164                            iv%instid(i)%qhl(k,n), &
 165                            iv%instid(i)%rcw(k,n), &
 166                            iv%instid(i)%rci(k,n), &
 167                            iv%instid(i)%rrn(k,n), &
 168                            iv%instid(i)%rsn(k,n), &
 169                            iv%instid(i)%rgr(k,n), &
 170                            iv%instid(i)%rhl(k,n)
 171                      end do ! end loop profile
 172                   else  ! no cloud
 173                      do k=1,iv%instid(i)%nlevels-1
 174                         write(unit=innov_rad_unit,fmt='(i3,2f10.2,f8.2,7f8.3)') &
 175                            k,  &
 176                            iv%instid(i)%pf(k,n), &
 177                            iv%instid(i)%pm(k,n), &
 178                            iv%instid(i)%tm(k,n), &
 179                            iv%instid(i)%qm(k,n), &
 180                            0.0, &
 181                            0.0, &
 182                            0.0, &
 183                            0.0, &
 184                            0.0, &
 185                            0.0
 186                      end do ! end loop profile
 187                   end if  ! end if crtm_cloud
 188 #endif
 189                end if ! end if rtm_option_crtm
 190
 191             end if  ! end if write_profile
 192
 193             if ( rtm_option == rtm_option_crtm .and. write_jacobian) then
 194 #ifdef CRTM
 195                write(unit=innov_rad_unit,fmt=*) &
 196                   'channel level halfp(mb) t(k) q(g/kg) water(mm) ice(mm) rain(mm) snow(mm) graupel(mm) hail(mm)'
 197                if (crtm_cloud) then
 198                   do l=1,iv%instid(i)%nchan
 199                      do k=1,iv%instid(i)%nlevels-1
 200                         write(unit=innov_rad_unit,fmt='(i5,i3,f10.2,8f8.3,6f8.3)') &
 201                            l, k,  &
 202                            iv%instid(i)%pm(k,n), &
 203                            iv%instid(i)%t_jacobian(l,k,n), &
 204                            iv%instid(i)%q_jacobian(l,k,n), &
 205                            iv%instid(i)%water_jacobian(l,k,n), &
 206                            iv%instid(i)%ice_jacobian(l,k,n), &
 207                            iv%instid(i)%rain_jacobian(l,k,n), &
 208                            iv%instid(i)%snow_jacobian(l,k,n), &
 209                            iv%instid(i)%graupel_jacobian(l,k,n), &
 210                            iv%instid(i)%hail_jacobian(l,k,n), &
 211                            iv%instid(i)%water_r_jacobian(l,k,n), &
 212                            iv%instid(i)%ice_r_jacobian(l,k,n), &
 213                            iv%instid(i)%rain_r_jacobian(l,k,n), &
 214                            iv%instid(i)%snow_r_jacobian(l,k,n), &
 215                            iv%instid(i)%graupel_r_jacobian(l,k,n), &
 216                            iv%instid(i)%hail_r_jacobian(l,k,n)
 217                     end do ! end loop profile
 218                  end do ! end loop channels
 219                else  ! no cloud
 220                   do l=1,iv%instid(i)%nchan
 221                      do k=1,iv%instid(i)%nlevels-1
 222                         write(unit=innov_rad_unit,fmt='(i5,i3,f10.2,8f8.3,6f8.3)') &
 223                            l, k,  &
 224                            iv%instid(i)%pm(k,n), &
 225                            iv%instid(i)%t_jacobian(l,k,n), &
 226                            iv%instid(i)%q_jacobian(l,k,n), &
 227                            0., &
 228                            0., &
 229                            0., &
 230                            0., &
 231                            0., &
 232                            0., &
 233                            0., &
 234                            0., &
 235                            0., &
 236                            0., &
 237                            0., &
 238                            0.
 239                      end do ! end loop profile
 240                   end do ! end loop channels
 241                end if  ! end if crtm_cloud
 242 #endif
 243             end if !  end if write_jacobian
 244
 245          end if ! end if proc_domain
 246       end do ! end do pixels
 247       close(unit=innov_rad_unit)
 248       call da_free_unit(innov_rad_unit)
 249    end do ! end do instruments
 250
 251    if (trace_use) call da_trace_exit("da_write_iv_rad_ascii")
 252
 253 end subroutine da_write_iv_rad_ascii
 254