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    logical            :: amsr2, ahi
  22    real               :: cip ! to output cloud-ice path
  23    integer            :: cloudflag ! to output cloudflag
  24    integer, dimension(1) :: maxl
  25    integer, allocatable :: level_max(:)
  26
  27    real, allocatable :: dtransmt(:,:), transmt_jac(:,:), transmt(:,:), lod(:,:), lod_jac(:,:)
  28
  29    if (trace_use) call da_trace_entry("da_write_iv_rad_ascii")
  30
  31    write(unit=message(1),fmt='(A)') 'Writing radiance OMB ascii file'
  32    call da_message(message(1:1))
  33
  34    do i = 1, iv%num_inst
  35       if (iv%instid(i)%num_rad < 1) cycle
  36
  37       ! count number of obs within the loc%proc_domain
  38       ! ---------------------------------------------
  39       ndomain = 0
  40       do n =1,iv%instid(i)%num_rad
  41          if (iv%instid(i)%info%proc_domain(1,n)) then
  42             ndomain = ndomain + 1
  43          end if
  44       end do
  45       if (ndomain < 1) cycle
  46
  47       if (rtm_option==rtm_option_crtm .and. write_jacobian ) then
  48          allocate ( dtransmt(iv%instid(i)%nchan,iv%instid(i)%nlevels) )
  49          allocate ( transmt_jac(iv%instid(i)%nchan,iv%instid(i)%nlevels) )
  50          allocate ( transmt(iv%instid(i)%nchan,iv%instid(i)%nlevels) )
  51          allocate ( lod(iv%instid(i)%nchan,iv%instid(i)%nlevels) )
  52          allocate ( lod_jac(iv%instid(i)%nchan,iv%instid(i)%nlevels) )
  53       end if
  54
  55       ! Get variables for maximum level of weighting function
  56       !  Only filled if calc_weightfunc .and. use_crtm_kmatrix both are true
  57       if ( rtm_option == rtm_option_crtm .and. calc_weightfunc .and. use_crtm_kmatrix ) then
  58          allocate(level_max(iv%instid(i)%nchan))
  59       endif
  60
  61       amsr2 = index(iv%instid(i)%rttovid_string,'amsr2') > 0
  62       ahi   = index(iv%instid(i)%rttovid_string,'ahi') > 0
  63
  64       write(unit=filename, fmt='(i2.2,a,i4.4)') it,'_inv_'//trim(iv%instid(i)%rttovid_string)//'.', myproc
  65
  66       call da_get_unit(innov_rad_unit)
  67       open(unit=innov_rad_unit,file=trim(filename),form='formatted',iostat=ios)
  68       if (ios /= 0 ) then
  69          call da_error(__FILE__,__LINE__, &
  70             (/"Cannot open innovation radiance file"//filename/))
  71       Endif
  72       write(unit=innov_rad_unit,fmt='(a,a,i7,a,i5,a)') trim(iv%instid(i)%rttovid_string), &
  73                         ' number-of-pixels : ', ndomain, &
  74                         ' channel-number-of-each-pixel : ', iv%instid(i)%nchan, &
  75                         ' index-of-channels : '
  76       write(unit=innov_rad_unit,fmt='(10i5)') iv%instid(i)%ichan
  77       if ( amsr2 ) then
  78          write(unit=innov_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask  elv lat lon  satzen satazi solzen solazi clw'
  79       else
  80          write(unit=innov_rad_unit,fmt='(a)') ' pixel-info : i date scanpos landsea_mask  elv lat lon  satzen satazi solzen solazi'
  81       end if
  82       write(unit=innov_rad_unit,fmt='(a)') ' grid%xb-surf-info : i t2m mr2m(ppmv) u10 v10 ps ts smois tslb snowh isflg &
  83                     & soiltyp vegtyp vegfra elev clwp cloud_frac cip cloudflag'
  84       ndomain = 0
  85       do n =1,iv%instid(i)%num_rad
  86          if (iv%instid(i)%info%proc_domain(1,n)) then
  87             ndomain=ndomain+1
  88             if ( amsr2 ) then ! write out clw
  89                write(unit=innov_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,6f8.2,f8.3)') 'INFO : ', ndomain, &
  90                                 iv%instid(i)%info%date_char(n), &
  91                                 iv%instid(i)%scanpos(n),        &
  92                                 iv%instid(i)%landsea_mask(n),   &
  93                                 iv%instid(i)%info%elv(n),       &
  94                                 iv%instid(i)%info%lat(1,n),     &
  95                                 iv%instid(i)%info%lon(1,n),     &
  96                                 iv%instid(i)%satzen(n),         &
  97                                 iv%instid(i)%satazi(n),         &
  98                                 iv%instid(i)%solzen(n),         &
  99                                 iv%instid(i)%solazi(n),         &
 100                                 iv%instid(i)%clw(n)
 101             else ! no clw info
 102                write(unit=innov_rad_unit,fmt='(a,i7,2x,a,i6,i3,f6.0,6f8.2)') 'INFO : ', ndomain, &
 103                                 iv%instid(i)%info%date_char(n), &
 104                                 iv%instid(i)%scanpos(n),        &
 105                                 iv%instid(i)%landsea_mask(n),   &
 106                                 iv%instid(i)%info%elv(n),       &
 107                                 iv%instid(i)%info%lat(1,n),     &
 108                                 iv%instid(i)%info%lon(1,n),     &
 109                                 iv%instid(i)%satzen(n),         &
 110                                 iv%instid(i)%satazi(n),         &
 111                                 iv%instid(i)%solzen(n),         &
 112                                 iv%instid(i)%solazi(n)
 113             end if
 114             select case (iv%instid(i)%isflg(n))
 115             case (0) ;
 116                surftype = ' SEA : '
 117             case (1) ;
 118                surftype = ' ICE : '
 119             case (2) ;
 120                surftype = 'LAND : '
 121             case (3) ;
 122                surftype = 'SNOW : '
 123             case (4) ;
 124                surftype = 'MSEA : '
 125             case (5) ;
 126                surftype = 'MICE : '
 127             case (6) ;
 128                surftype = 'MLND : '
 129             case (7) ;
 130                surftype = 'MSNO : '
 131             end select
 132
 133             ! Output cloud-ice path, cloudflag, pressure of weighting function peak
 134             if (rtm_option==rtm_option_crtm .and. crtm_cloud ) then
 135                cip = iv%instid(i)%cip(n)
 136             else
 137                cip = -9999.0
 138             endif
 139             if ( ahi ) then
 140                cloudflag = iv%instid(i)%cloudflag(n)
 141             else
 142                cloudflag = -999
 143             endif
 144
 145             ! %nlevels is unstaggered, so subtract 1 to get number of mass points
 146             ! Order is top down in CRTM
 147             if ( rtm_option == rtm_option_crtm .and. calc_weightfunc .and. use_crtm_kmatrix ) then
 148                level_max(:) = -1
 149                do l=1,iv%instid(i)%nchan
 150                   maxl(:) = maxloc(iv%instid(i)%der_trans(l,:,n)) ! returns index of maximum value; if a tie, returns first occurrence of maximum
 151                   level_max(l) = maxl(1) ! model level of weighting function peak for this pixel and channel, going from top-->down
 152                enddo
 153               !level_max(:) = ( iv%instid(i)%nlevels - 1 ) - level_max(:) + 1 ! make order bottom-up
 154             endif
 155
 156             write(unit=innov_rad_unit,fmt='(a,i7,9f10.2,3i3,f8.3,f10.2,f8.3,f15.5,f15.5,i7)') surftype, n, &
 157                              iv%instid(i)%t2m(n), &
 158                              iv%instid(i)%mr2m(n),   &
 159                              iv%instid(i)%u10(n), &
 160                              iv%instid(i)%v10(n),  &
 161                              iv%instid(i)%ps(n),  &
 162                              iv%instid(i)%ts(n),  &
 163                              iv%instid(i)%smois(n),  &
 164                              iv%instid(i)%tslb(n),  &
 165                              iv%instid(i)%snowh(n), &
 166                              iv%instid(i)%isflg(n), &
 167                              nint(iv%instid(i)%soiltyp(n)), &
 168                              nint(iv%instid(i)%vegtyp(n)), &
 169                              iv%instid(i)%vegfra(n), &
 170                              iv%instid(i)%elevation(n), &
 171                              iv%instid(i)%clwp(n), &
 172                              iv%instid(i)%cloud_frac(n), &
 173                              cip, &
 174                              cloudflag
 175
 176             write(unit=innov_rad_unit,fmt='(a)') 'OBS  : '
 177             write(unit=innov_rad_unit,fmt='(10f11.2)') ob%instid(i)%tb(:,n)
 178             write(unit=innov_rad_unit,fmt='(a)') 'BAK  : '
 179             write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%tb_xb(:,n)
 180             if (rtm_option==rtm_option_crtm .and. crtm_cloud .and. (amsr2 .or. ahi) ) then
 181                write(unit=innov_rad_unit,fmt='(a)') 'BAK_clr : '
 182                write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%tb_xb_clr(:,n)
 183             endif
 184             if ( rtm_option == rtm_option_crtm .and. calc_weightfunc .and. use_crtm_kmatrix ) then
 185                write(unit=innov_rad_unit,fmt='(a)') 'WEIGHTFUNC_PEAK : '
 186                write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%pm( (/level_max/), n )
 187             endif
 188             write(unit=innov_rad_unit,fmt='(a)') 'IVBC : '
 189             write(unit=innov_rad_unit,fmt='(10f11.2)')  iv%instid(i)%tb_inv(:,n)
 190             write(unit=innov_rad_unit,fmt='(a)') 'EMS  : '
 191             write(unit=innov_rad_unit,fmt='(10f11.2)')  iv%instid(i)%emiss(1:iv%instid(i)%nchan,n)
 192             if (rtm_option==rtm_option_crtm .and. write_jacobian) then
 193                write(unit=innov_rad_unit,fmt='(a)') 'EMS_JACOBIAN : '
 194                write(unit=innov_rad_unit,fmt='(10f10.3)') iv%instid(i)%emiss_jacobian(1:iv%instid(i)%nchan,n)
 195             end if
 196             write(unit=innov_rad_unit,fmt='(a)') 'ERR  : '
 197             write(unit=innov_rad_unit,fmt='(10f11.2)') iv%instid(i)%tb_error(:,n)
 198             write(unit=innov_rad_unit,fmt='(a)') 'QC   : '
 199             write(unit=innov_rad_unit,fmt='(10i11)') iv%instid(i)%tb_qc(:,n)
 200
 201             if (write_profile) then
 202                nlevelss  = iv%instid(i)%nlevels
 203                if ( rtm_option == rtm_option_rttov ) then
 204 #ifdef RTTOV
 205                   ! first, write RTTOV levels
 206                   write(unit=innov_rad_unit,fmt='(a)') 'RTM_level pres(mb) T(k) Q(ppmv)'
 207                   do k = 1, nlevelss
 208                      write(unit=innov_rad_unit,fmt='(i3,f10.2,f8.2,e11.4)') &
 209                         k, &                             ! RTTOV levels
 210                         coefs(i) % coef % ref_prfl_p(k) , &
 211                         iv%instid(i)%t(k,n) , &
 212                         iv%instid(i)%mr(k,n)
 213                   end do  ! end loop RTTOV level
 214                   ! second, write WRF model levels
 215                   write(unit=innov_rad_unit,fmt='(a)') &
 216                      'WRF_level pres(mb) T(k) q(g/kg) clw(g/kg) rain(g/kg)'
 217                   do k=kts,kte
 218                      write(unit=innov_rad_unit,fmt='(i3,f10.2,f8.2,3e11.4)') &
 219                         k,  &                     ! WRF model levels
 220                         iv%instid(i)%pm(k,n) , &
 221                         iv%instid(i)%tm(k,n) , &
 222                         iv%instid(i)%qm(k,n)*1000 , &
 223                         iv%instid(i)%qcw(k,n)*1000.0, &
 224                         iv%instid(i)%qrn(k,n)*1000.0
 225                   end do ! end loop model level
 226 #endif
 227                end if ! end if rtm_option_rttov
 228
 229                if ( rtm_option == rtm_option_crtm ) then
 230 #ifdef CRTM
 231                   write(unit=innov_rad_unit,fmt='(a)') &
 232                      'level fullp(mb) halfp(mb) t(k) q(g/kg) water(mm) ice(mm) rain(mm) snow(mm) graupel(mm) hail(mm)'
 233                   if (crtm_cloud) then
 234                      do k=1,iv%instid(i)%nlevels-1
 235                         write(unit=innov_rad_unit,fmt='(i3,2f10.2,f8.2,13f8.3)') &
 236                            k,  &
 237                            iv%instid(i)%pf(k,n), &
 238                            iv%instid(i)%pm(k,n), &
 239                            iv%instid(i)%tm(k,n), &
 240                            iv%instid(i)%qm(k,n), &
 241                            iv%instid(i)%qcw(k,n), &
 242                            iv%instid(i)%qci(k,n), &
 243                            iv%instid(i)%qrn(k,n), &
 244                            iv%instid(i)%qsn(k,n), &
 245                            iv%instid(i)%qgr(k,n), &
 246                            iv%instid(i)%qhl(k,n), &
 247                            iv%instid(i)%rcw(k,n), &
 248                            iv%instid(i)%rci(k,n), &
 249                            iv%instid(i)%rrn(k,n), &
 250                            iv%instid(i)%rsn(k,n), &
 251                            iv%instid(i)%rgr(k,n), &
 252                            iv%instid(i)%rhl(k,n)
 253                      end do ! end loop profile
 254                   else  ! no cloud
 255                      do k=1,iv%instid(i)%nlevels-1
 256                         write(unit=innov_rad_unit,fmt='(i3,2f10.2,f8.2,7f8.3)') &
 257                            k,  &
 258                            iv%instid(i)%pf(k,n), &
 259                            iv%instid(i)%pm(k,n), &
 260                            iv%instid(i)%tm(k,n), &
 261                            iv%instid(i)%qm(k,n), &
 262                            0.0, &
 263                            0.0, &
 264                            0.0, &
 265                            0.0, &
 266                            0.0, &
 267                            0.0
 268                      end do ! end loop profile
 269                   end if  ! end if crtm_cloud
 270 #endif
 271                end if ! end if rtm_option_crtm
 272
 273             end if  ! end if write_profile
 274
 275             if ( rtm_option == rtm_option_crtm .and. write_jacobian) then
 276 #ifdef CRTM
 277
 278                if ( calc_weightfunc ) then
 279                   dtransmt(:,:)    = iv%instid(i)%der_trans(:,:,n)
 280                   transmt(:,:)     = iv%instid(i)%trans(:,:,n)
 281                   transmt_jac(:,:) = iv%instid(i)%trans_jacobian(:,:,n)
 282                   lod(:,:)         = iv%instid(i)%lod(:,:,n)
 283                   lod_jac(:,:)     = iv%instid(i)%lod_jacobian(:,:,n)
 284                else
 285                   dtransmt(:,:)    = 0.0
 286                   transmt(:,:)     = 0.0
 287                   transmt_jac(:,:) = 0.0
 288                   lod(:,:)         = 0.0
 289                   lod_jac(:,:)     = 0.0
 290                end if
 291
 292                write(unit=innov_rad_unit,fmt='(a)') &
 293                   'channel level halfp(mb) t(k) q(g/kg) der_trans trans_jac trans lod_jac lod water(mm) ice(mm) rain(mm) snow(mm) graupel(mm) hail(mm)'
 294                if (crtm_cloud) then
 295                   do l=1,iv%instid(i)%nchan
 296                      do k=1,iv%instid(i)%nlevels-1
 297                         write(unit=innov_rad_unit,fmt='(i5,i3,f10.2,13f14.7,6f14.7)') &
 298                            l, k,  &
 299                            iv%instid(i)%pm(k,n), &
 300                            iv%instid(i)%t_jacobian(l,k,n), &
 301                            iv%instid(i)%q_jacobian(l,k,n), &
 302                            dtransmt(l,k),&
 303                            transmt_jac(l,k),&
 304                            transmt(l,k),&
 305                            lod_jac(l,k),&
 306                            lod(l,k),&
 307                            iv%instid(i)%water_jacobian(l,k,n), &
 308                            iv%instid(i)%ice_jacobian(l,k,n), &
 309                            iv%instid(i)%rain_jacobian(l,k,n), &
 310                            iv%instid(i)%snow_jacobian(l,k,n), &
 311                            iv%instid(i)%graupel_jacobian(l,k,n), &
 312                            iv%instid(i)%hail_jacobian(l,k,n), &
 313                            iv%instid(i)%water_r_jacobian(l,k,n), &
 314                            iv%instid(i)%ice_r_jacobian(l,k,n), &
 315                            iv%instid(i)%rain_r_jacobian(l,k,n), &
 316                            iv%instid(i)%snow_r_jacobian(l,k,n), &
 317                            iv%instid(i)%graupel_r_jacobian(l,k,n), &
 318                            iv%instid(i)%hail_r_jacobian(l,k,n)
 319                     end do ! end loop profile
 320                  end do ! end loop channels
 321                else  ! no cloud
 322                   do l=1,iv%instid(i)%nchan
 323                      do k=1,iv%instid(i)%nlevels-1
 324                         write(unit=innov_rad_unit,fmt='(i5,i3,f10.2,13f14.7,6f14.7)') &
 325                            l, k,  &
 326                            iv%instid(i)%pm(k,n), &
 327                            iv%instid(i)%t_jacobian(l,k,n), &
 328                            iv%instid(i)%q_jacobian(l,k,n), &
 329                            dtransmt(l,k),&
 330                            transmt_jac(l,k),&
 331                            transmt(l,k),&
 332                            lod_jac(l,k),&
 333                            lod(l,k),&
 334                            0., &
 335                            0., &
 336                            0., &
 337                            0., &
 338                            0., &
 339                            0., &
 340                            0., &
 341                            0., &
 342                            0., &
 343                            0., &
 344                            0., &
 345                            0.
 346                      end do ! end loop profile
 347                   end do ! end loop channels
 348                end if  ! end if crtm_cloud
 349 #endif
 350             end if !  end if write_jacobian
 351
 352          end if ! end if proc_domain
 353       end do ! end do pixels
 354       if (rtm_option==rtm_option_crtm .and. write_jacobian ) then
 355          deallocate ( dtransmt )
 356          deallocate ( transmt_jac )
 357          deallocate ( transmt )
 358          deallocate ( lod )
 359          deallocate ( lod_jac )
 360       end if
 361       if ( rtm_option == rtm_option_crtm .and. calc_weightfunc .and. use_crtm_kmatrix ) then
 362          deallocate(level_max)
 363       endif
 364       close(unit=innov_rad_unit)
 365       call da_free_unit(innov_rad_unit)
 366    end do ! end do instruments
 367
 368    if (trace_use) call da_trace_exit("da_write_iv_rad_ascii")
 369
 370 end subroutine da_write_iv_rad_ascii
 371