var/da/da_radiance/da_get_innov_vector_crtmk.inc

   1 subroutine da_get_innov_vector_crtmk ( it, grid, ob, iv )
   2
   3    !---------------------------------------------------------------------------
   4    !  PURPOSE: Calculate innovation vector for radiance data.
   5    !
   6    !  METHOD:  d = y - H(x)
   7    !       1. interpolate grid%xb to obs location
   8    !       2. call foreward RTM to get simulated bright temperature
   9    !       3. obs BT - simulated BT
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    integer,           intent(in)    :: it       ! External iteration.
  15    type (domain),     intent(in)    :: grid     ! first guess state.
  16    type (y_type),     intent(inout) :: ob       ! Observation structure.
  17    type (iv_type),    intent(inout) :: iv       ! O-B structure.
  18
  19    integer :: n, icld  ! Loop counter.
  20    integer :: i, j, k  ! Index dimension.
  21    integer :: l        ! Index dimension.
  22    integer :: num_levs ! Number of obs levels.
  23    real    :: dx, dxm  ! Interpolation weights.
  24    real    :: dy, dym  ! Interpolation weights.
  25    integer :: alloc_status(40)
  26
  27    real, allocatable :: model_u10(:)
  28    real, allocatable :: model_v10(:)
  29    real, allocatable :: model_psfc(:)
  30    real              :: model_ptop
  31    real, allocatable :: model_ts(:)
  32    real, allocatable :: model_elv(:)
  33    real, allocatable :: model_smois(:)
  34    real, allocatable :: model_tslb(:)
  35    real, allocatable :: model_snowh(:)
  36    real, allocatable :: model_vegfra(:)
  37    real    :: model_isltyp, model_ivgtyp
  38    integer :: model_isflg
  39
  40    real    :: v_p(kms:kme)  ! Model value p at ob hor. location.
  41    real    :: model_qcw(kms:kme)
  42    real    :: model_rho(kms:kme)
  43    real, allocatable :: model_snow(:)  ! snow water equivalent, different from model_snowh,
  44                                        ! used in calculating reff_water
  45
  46    ! real    :: model_tm(kms:kme)
  47
  48    integer :: inst, nchanl, n1,n2
  49
  50    ! variables for computing clwp
  51    real    :: clw(kms:kme), dpf(kms:kme)
  52    real    :: clwp
  53
  54    ! CRTM local varaibles and types
  55    integer :: wmo_sensor_id,Error_Status, Allocate_Status
  56
  57    type (CRTM_RTSolution_type), allocatable :: RTSolution(:,:), RTSolution_K(:,:)
  58    type (CRTM_Atmosphere_type)   :: Atmosphere(1)
  59    type (CRTM_Surface_type)      :: Surface(1)
  60    type (CRTM_Atmosphere_type), allocatable :: Atmosphere_K(:,:)
  61    type (CRTM_Surface_type),    allocatable :: Surface_K(:,:)
  62    type (CRTM_GeometryInfo_type) :: GeometryInfo(1)
  63
  64    real :: t(1), a(1)
  65
  66    ! WHY? use argument it
  67    if (it==0) then; write(unit=stdout,fmt='(A)') "WHY? have argument it to"//__FILE__; end if
  68
  69    alloc_status (:) = 0
  70
  71    if (trace_use) call da_trace_entry("da_get_innov_vector_crtmk")
  72
  73    ! CRTM allocation
  74
  75    ! Atmosphere structure
  76    Atmosphere(1)%n_Layers=(kte-kts)+1   ! number of vertical levels
  77    Atmosphere(1)%n_Absorbers=2
  78    Atmosphere(1)%n_Clouds=0
  79    Atmosphere(1)%n_Aerosols=0
  80    if (crtm_cloud) Atmosphere(1)%n_Clouds=6
  81
  82    Error_Status = CRTM_Allocate_Atmosphere( Atmosphere(1)%n_Layers, &
  83                                             Atmosphere(1)%n_Absorbers, &
  84                                             Atmosphere(1)%n_Clouds, &
  85                                             Atmosphere(1)%n_Aerosols, &
  86                                             Atmosphere)
  87    if (Error_Status /= 0) then
  88        call da_error(__FILE__,__LINE__, &
  89          (/"Error in allocating CRTM Atmosphere Structure"/))
  90    end if
  91
  92    Atmosphere(1)%Absorber_ID(1)=H2O_ID
  93    Atmosphere(1)%Absorber_ID(2)=O3_ID
  94
  95    if (crtm_cloud) then
  96       Atmosphere(1)%Cloud(1)%Type=WATER_CLOUD
  97       Atmosphere(1)%Cloud(2)%Type=ICE_CLOUD
  98       Atmosphere(1)%Cloud(3)%Type=RAIN_CLOUD
  99       Atmosphere(1)%Cloud(4)%Type=SNOW_CLOUD
 100       Atmosphere(1)%Cloud(5)%Type=GRAUPEL_CLOUD
 101       Atmosphere(1)%Cloud(6)%Type=HAIL_CLOUD
 102    end if
 103
 104    !------------------------------------------------------
 105    ! [1.0] calculate the background bright temperature
 106    !-------------------------------------------------------
 107    do inst = 1, iv%num_inst                 ! loop for sensor
 108       ! if ( iv%instid(inst)%num_rad < 1 ) cycle
 109       if (iv%instid(inst)%info%n2 < iv%instid(inst)%info%n1) cycle
 110       num_levs  = kte-kts+1
 111       ! CRTM channel information structure
 112       ! Error_Status = CRTM_Set_ChannelInfo(Sensor_Descriptor(inst),ChannelInfo)
 113       ! if (Error_Status /= 0) then
 114       !   call da_error(__FILE__,__LINE__, (/"Error in calling CRTM_Set_ChannelInfo"/))
 115       ! end if
 116       nchanl    = ChannelInfo(inst)%n_channels
 117
 118       ! Allocate forward model solution RTSolution array to number of channels
 119       allocate (RTSolution(ChannelInfo(inst)%n_Channels,1),   &
 120                 RTSolution_K(ChannelInfo(inst)%n_Channels,1), &
 121                 Atmosphere_K(ChannelInfo(inst)%n_Channels,1), &
 122                 Surface_K(ChannelInfo(inst)%n_Channels,1),    &
 123                 STAT = Allocate_Status )
 124       if ( Allocate_Status /= 0 ) then
 125          call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution"/))
 126       end if
 127
 128       ! CRTM Surface Structure
 129       if (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsua') then
 130          wmo_sensor_id=WMO_AMSUA
 131       elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsub') then
 132          wmo_sensor_id=WMO_AMSUB
 133       elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='amsre') then
 134          wmo_sensor_id=WMO_AMSRE
 135       elseif (trim(crtm_sensor_name(rtminit_sensor(inst))) =='ssmi') then
 136          wmo_sensor_id=WMO_SSMI
 137       else
 138          wmo_sensor_id=INVALID_WMO_SENSOR_ID
 139       end if
 140
 141       Error_Status = CRTM_Allocate_Surface( nchanl,     &  ! Input
 142                                    Surface)  ! Output
 143       if (Error_Status /= 0) then
 144          call da_error(__FILE__,__LINE__, &
 145             (/"Error in allocating CRTM Surface Structure Structure"/))
 146       end if
 147       ! do n= 1, iv%instid(inst)%num_rad           ! loop for pixel
 148
 149       n1 = iv%instid(inst)%info%n1
 150       n2 = iv%instid(inst)%info%n2
 151
 152       allocate (model_u10(n1:n2))
 153       allocate (model_v10(n1:n2))
 154       allocate (model_psfc(n1:n2))
 155       allocate (model_ts(n1:n2))
 156       allocate (model_elv(n1:n2))
 157       allocate (model_smois(n1:n2))
 158       allocate (model_tslb(n1:n2))
 159       allocate (model_snowh(n1:n2))
 160       allocate (model_snow(n1:n2))
 161       allocate (model_vegfra(n1:n2))
 162
 163       model_u10(:)    = 0.0
 164       model_v10(:)    = 0.0
 165       model_psfc(:)   = 0.0
 166       model_ts(:)     = 0.0
 167       model_elv(:)    = 0.0
 168       model_smois(:)  = 0.0
 169       model_tslb(:)   = 0.0
 170       model_snowh(:)  = 0.0
 171       model_snow(:)   = 0.0
 172       model_vegfra(:) = 0.0
 173
 174       do n=n1,n2
 175          ! if ( n > iv%instid(inst)%num_rad ) exit
 176
 177          ! [1.1] Get horizontal interpolation weights:
 178
 179          i = iv%instid(inst)%info%i(1,n)
 180          j = iv%instid(inst)%info%j(1,n)
 181          dx = iv%instid(inst)%info%dx(1,n)
 182          dy = iv%instid(inst)%info%dy(1,n)
 183          dxm = iv%instid(inst)%info%dxm(1,n)
 184          dym = iv%instid(inst)%info%dym(1,n)
 185
 186          ! horizontal interpolate grid%xb pressure to ob position for every grid%xb layer
 187          ! get CRTM pressure Layers
 188          do k=kts,kte ! from bottem to top
 189             v_p(k) = dym*(dxm*grid%xb%p(i,j  ,k) + dx*grid%xb%p(i+1,j  ,k)) &
 190                 + dy *(dxm*grid%xb%p(i,j+1,k) + dx*grid%xb%p(i+1,j+1,k))
 191             v_p(k) = 0.01*v_p(k)  ! convert Pa to hPa
 192             Atmosphere(1)%Pressure(kte-k+1)=v_p(k) ! from top to bottom
 193          end do
 194
 195          ! determine surface type of obs location
 196          !-----------------------------------------
 197          call da_detsurtyp ( grid%xb%snow, grid%xb%xice, grid%xb%landmask,  &
 198             grid%xb%ivgtyp, grid%xb%isltyp, &
 199             ims, ime, jms, jme, &
 200             i, j, dx, dy, dxm, dym, &
 201             model_isflg,model_ivgtyp, model_isltyp, &
 202             Surface(1)%Water_Coverage, Surface(1)%Ice_Coverage, &
 203             Surface(1)%Land_Coverage, Surface(1)%Snow_Coverage )
 204
 205          call da_interp_lin_2d_partial (grid%xb%snow, iv%instid(inst)%info, 1, n, n, model_snow(n:n) )
 206
 207          ! [1.2] Interpolate horizontally to ob:
 208          do k=kts,kte ! from bottom to top
 209             call da_interp_lin_2d_partial (grid%xb%t(:,:,k), iv%instid(inst)%info,k,n,n,t)
 210             call da_interp_lin_2d_partial (grid%xb%q(:,:,k), iv%instid(inst)%info,k,n,n,a)
 211             Atmosphere(1)%Temperature(kte-k+1) = t(1)
 212             Atmosphere(1)%Absorber(kte-k+1,1)  = a(1)
 213             Atmosphere(1)%Absorber(kte-k+1,1) = 1000.0*Atmosphere(1)%Absorber(kte-k+1,1) ! in g/kg
 214
 215             ! NOTE: WRF high-level q values seems too big, replaced by constants
 216             if (v_p(k) < 75.0) Atmosphere(1)%Absorber(kte-k+1,1) = 0.001
 217
 218             call da_interp_lin_2d_partial (grid%xb%qcw(:,:,k), iv%instid(inst)%info,k,n,n, model_qcw(kte-k+1:kte-k+1))
 219
 220             if (crtm_cloud) then
 221                Atmosphere(1)%Cloud(1)%Water_Content(kte-k+1) = model_qcw(kte-k+1)
 222
 223                call da_interp_lin_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n, &
 224                   Atmosphere(1)%Cloud(2)%Water_Content(kte-k+1:kte-k+1) )
 225
 226                call da_interp_lin_2d_partial (grid%xb%qrn(:,:,k), iv%instid(inst)%info,k,n,n, &
 227                   Atmosphere(1)%Cloud(3)%Water_Content(kte-k+1:kte-k+1) )
 228
 229                call da_interp_lin_2d_partial (grid%xb%qsn(:,:,k), iv%instid(inst)%info, k,n,n, &
 230                   Atmosphere(1)%Cloud(4)%Water_Content(kte-k+1:kte-k+1) )
 231
 232                call da_interp_lin_2d_partial (grid%xb%qgr(:,:,k), iv%instid(inst)%info, k,n,n, &
 233                   Atmosphere(1)%Cloud(5)%Water_Content(kte-k+1:kte-k+1) )
 234
 235                Atmosphere(1)%Cloud(6)%Water_Content(kte-k+1) = 0.0
 236
 237                call da_interp_lin_2d_partial (grid%xb%rho(:,:,k), iv%instid(inst)%info, k,n,n, &
 238                   model_rho(k:k) )
 239
 240                call da_cld_eff_radius(Atmosphere(1)%Temperature(kte-k+1),model_rho(k),&
 241                                       Atmosphere(1)%Cloud(2)%Water_Content(kte-k+1),  &  !qci
 242                                       Atmosphere(1)%Cloud(3)%Water_Content(kte-k+1),  &  !qrn
 243                                       Atmosphere(1)%Cloud(4)%Water_Content(kte-k+1),  &  !qsn
 244                                       Atmosphere(1)%Cloud(5)%Water_Content(kte-k+1),  &  !qgr
 245                                       model_snow(n),                                  &
 246                                       Surface(1)%Ice_Coverage, Surface(1)%Land_Coverage, 1, &
 247                                       Atmosphere(1)%Cloud(1)%Effective_Radius(kte-k+1), &
 248                                       Atmosphere(1)%Cloud(2)%Effective_Radius(kte-k+1), &
 249                                       Atmosphere(1)%Cloud(3)%Effective_Radius(kte-k+1), &
 250                                       Atmosphere(1)%Cloud(4)%Effective_Radius(kte-k+1), &
 251                                       Atmosphere(1)%Cloud(5)%Effective_Radius(kte-k+1) )
 252
 253                ! reset the da_cld_eff_radius calcualted effective radius to constants if desired
 254
 255                Atmosphere(1)%Cloud(1)%Effective_Radius(kte-k+1)=0.01  ! in mm
 256                Atmosphere(1)%Cloud(2)%Effective_Radius(kte-k+1)=0.03  ! in mm
 257                !Atmosphere(1)%Cloud(3)%Effective_Radius(kte-k+1)=0.3
 258                !Atmosphere(1)%Cloud(4)%Effective_Radius(kte-k+1)=0.6
 259                !Atmosphere(1)%Cloud(5)%Effective_Radius(kte-k+1)=0.6
 260                Atmosphere(1)%Cloud(6)%Effective_Radius(kte-k+1)=0.6
 261
 262             end if
 263          end do
 264
 265          call da_interp_lin_2d_partial (grid%xb%u10,  iv%instid(inst)%info, 1, n, n, model_u10(n:n))
 266          call da_interp_lin_2d_partial (grid%xb%v10,  iv%instid(inst)%info, 1, n, n, model_v10(n:n))
 267          call da_interp_lin_2d_partial (grid%xb%psfc, iv%instid(inst)%info, 1, n, n, model_psfc(n:n))
 268
 269          model_psfc(n) = 0.01*model_psfc(n)           ! convert to hPa
 270          model_ptop = 0.01*grid%xb%ptop
 271
 272          ! get CRTM levels (0.005hPa at top) /model full level
 273          Atmosphere(1)%Level_Pressure(0)=model_ptop   ! to sigma level 51-->sigmaf=0
 274          Atmosphere(1)%Level_Pressure(Atmosphere(1)%n_Layers)=model_psfc(n) ! to sigma level 1->sigmaf=1
 275
 276          do k=kts+1,kte
 277             Atmosphere(1)%Level_Pressure(kte-k+1)= grid%xb%sigmaf(k)*(model_psfc(n)-model_ptop)+model_ptop
 278          end do
 279
 280          ! convert cloud content unit from kg/kg to kg/m^2
 281          if (crtm_cloud) then
 282             do k=kts,kte
 283                do icld=1,Atmosphere(1)%n_Clouds
 284                   Atmosphere(1)%Cloud(icld)%Water_Content(k)= Atmosphere(1)%Cloud(icld)%Water_Content(k)* &
 285                      (Atmosphere(1)%Level_Pressure(k)- Atmosphere(1)%Level_Pressure(k-1))*100.0/gravity
 286                end do
 287             end do
 288          end if
 289
 290          if ( model_isflg == 0 ) then   ! over sea using SST
 291             call da_interp_lin_2d_partial (grid%xb % tgrn, iv%instid(inst)%info, 1, n, n, model_ts(n:n))
 292          else
 293             call da_interp_lin_2d_partial (grid%xb % tsk,  iv%instid(inst)%info, 1, n, n, model_ts(n:n))
 294          end if
 295
 296          call da_interp_lin_2d_partial (grid%xb % terr, iv%instid(inst)%info, 1, n, n, model_elv(n:n))
 297
 298          ! variables for emissivity calculations
 299          !----------------------------------------
 300          call da_interp_lin_2d_partial (grid%xb % smois,  iv%instid(inst)%info, 1, n, n, model_smois(n:n) )
 301          call da_interp_lin_2d_partial (grid%xb % tslb,   iv%instid(inst)%info, 1, n, n, model_tslb(n:n) )
 302          call da_interp_lin_2d_partial (grid%xb % snowh,  iv%instid(inst)%info, 1, n, n, model_snowh(n:n) )
 303          call da_interp_lin_2d_partial (grid%xb % vegfra, iv%instid(inst)%info, 1, n, n, model_vegfra(n:n) )
 304
 305          ! model_snowh(n) = model_snowh(n)*100.0   ! convert from m to mm
 306          model_vegfra(n) = 0.01*model_vegfra(n)  ! convert range to 0~1
 307
 308          ! ADD for computing cloud liquid water path (mm) from guess
 309          clwp = 0.0
 310          do k = kts,kte ! from top to bottom
 311             dpf(k) = 100.0*(Atmosphere(1)%level_pressure(k) - Atmosphere(1)%level_pressure(k-1))
 312             clw  (k) = model_qcw(k)*dpf(k)/gravity ! kg/m2 or mm
 313             if (Atmosphere(1)%pressure(k)<100.0) clw(k) = 0.0
 314             clwp  = clwp + clw(k)
 315          end do
 316
 317          ! CRTM GeometryInfo Structure
 318          GeometryInfo(1)%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
 319          GeometryInfo(1)%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
 320          ! GeometryInfo(1)%Satellite_Height=830.0
 321          ! GeometryInfo(1)%Sensor_Scan_Angle=
 322          ! GeometryInfo(1)%Sensor_Zenith_Angle=
 323          ! GeometryInfo(1)%Sensor_Scan_Angle=
 324          ! GeometryInfo(1)%Source_Zenith_Angle=
 325
 326          ! CRTM Surface parameter data
 327
 328          if (Surface(1)%Land_Coverage > 0.0) then
 329             Surface(1)%Land_Type=GRASS_SOIL           ! land type (User guide appendix 3)
 330             Surface(1)%Land_Temperature=model_ts(n)      ! K
 331             Surface(1)%Soil_Moisture_Content= model_smois(n) !0.05    ! volumetric water content (g/cm**3)
 332             ! Surface(1)%Canopy_Water_Content=0.05      ! gravimetric water content
 333             Surface(1)%Vegetation_Fraction=model_vegfra(n)
 334             Surface(1)%Soil_Temperature=model_tslb(n)
 335          end if
 336          if (Surface(1)%Water_Coverage > 0.0) then
 337             ! Surface%Water_Type=SEA_WATER          ! (Currently NOT used)
 338             Surface(1)%Water_Temperature=model_ts(n)     ! K
 339             Surface(1)%Wind_Speed=sqrt(model_u10(n)**2+model_v10(n)**2)  ! m/sec
 340             ! surface(1)%Wind_Direction=0.0            ! NOT used
 341             Surface(1)%Salinity=33.0                   ! ppmv
 342          end if
 343
 344          if (Surface(1)%Snow_Coverage > 0.0) then
 345             Surface(1)%Snow_Type=NEW_SNOW             ! User guide appendix 3
 346             Surface(1)%Snow_Temperature=model_ts(n)      ! K
 347             Surface(1)%Snow_Depth=model_snowh(n)         ! mm
 348             ! Surface(1)%Snow_Density=0.2               ! g/cm**3
 349             ! Surface(1)%Snow_Grain_Size=2.0            ! mm
 350          end if
 351          if (Surface(1)%Ice_Coverage > 0.0) then
 352             ! Surface(1)%Ice_Type=FRESH_ICE             ! NO Table offered, single example is FRESH_ICE
 353             Surface(1)%Ice_Temperature=model_ts(n)       ! K
 354             Surface(1)%Ice_Thickness=10.0              ! mm
 355             ! Surface(1)%Ice_Density=0.9                ! g/cm**3
 356             ! Surface(1)%Ice_Roughness=0.0               ! NO Table offered, single example is ZERO
 357
 358          end if
 359          if (nchanl > 0) then
 360             Surface(1)%SensorData%n_channels = nchanl
 361             Surface(1)%SensorData%Sensor_ID  = wmo_sensor_id
 362             Surface(1)%SensorData%Tb(1:nchanl) = ob%instid(inst)%tb(1:nchanl,n)
 363          end if
 364
 365          ! CRTM surface/atmosphere K initialization
 366          do l = 1, ChannelInfo(inst)%n_Channels
 367             ! -- Copy the adjoint atmosphere structure
 368             Error_Status = CRTM_Assign_Atmosphere( Atmosphere, Atmosphere_K(l,:) )
 369
 370             if ( Error_Status /= 0 ) then
 371                call da_error(__FILE__,__LINE__,  &
 372                   (/"Error copying Atmosphere_K structure"/))
 373             end if
 374
 375             ! -- Copy the adjoint surface structure
 376             Error_Status = CRTM_Assign_Surface( Surface, Surface_K(l,:) )
 377
 378             if ( Error_Status /= 0 ) then
 379                call da_error(__FILE__,__LINE__, &
 380                   (/"Error copying Surface_K structure"/))
 381             end if
 382          end do
 383
 384          ! -- Zero the Adjoint outputs
 385          ! Important: adjoint variables must be initialized
 386          call CRTM_Zero_Atmosphere( Atmosphere_K )
 387          call CRTM_Zero_Surface( Surface_K )
 388
 389          ! [1.3] assign tb = R^-1 Re :
 390
 391          RTSolution_K(:,1)%brightness_temperature = 1.
 392          RTSolution_K(:,1)%radiance = 0.
 393
 394          ! [1.3] Call RTM K-Matrix model
 395
 396          call da_crtm_k(1, nchanl, 1, Atmosphere,   &
 397                             Surface,      &
 398                             RTSolution_K,&
 399                             GeometryInfo, &
 400                             ChannelInfo(inst),  &
 401                             Atmosphere_K,&
 402                             Surface_K,   &
 403                             RTSolution) !,   &
 404                             !Options = Options)
 405
 406          !----------------------------------------------------------------
 407          ! [2.0] calculate components of innovation vector:
 408          !----------------------------------------------------------------
 409          do k = 1, nchanl
 410             if ( iv%instid(inst)%tb_inv(k,n) > missing_r ) then
 411                iv%instid(inst)%tb_xb(k,n)  = RTSolution(k,1)%Brightness_Temperature
 412                iv%instid(inst)%tb_inv(k,n) = &
 413                ob%instid(inst)%tb(k,n) - RTSolution(k,1)%Brightness_Temperature
 414             else
 415                iv%instid(inst)%tb_xb(k,n)    = RTSolution(k,1)%Brightness_Temperature
 416                iv%instid(inst)%tb_inv(k,n)   = missing_r
 417             end if
 418             iv%instid(inst)%emiss(k,n) = RTSolution(k,1)%surface_emissivity
 419             ! surface Jacobian
 420             iv%instid(inst)%ts_jacobian(k,n) = Surface_k(k,1)%water_temperature
 421             iv%instid(inst)%windspeed_jacobian(k,n) = Surface_k(k,1)%wind_speed
 422             iv%instid(inst)%emiss_jacobian(k,n) = RTSolution_k(k,1)%surface_emissivity
 423          end do
 424
 425          !----------------------------------------------------------------
 426          ! [3.0] store base state and Jacobian to innovation structure
 427          !----------------------------------------------------------------
 428          ! full level pressures
 429          iv%instid(inst)%pf(0,n)  = Atmosphere(1)%level_pressure(0)
 430          do k=1,Atmosphere(1)%n_layers
 431             iv%instid(inst)%pm(k,n)  = Atmosphere(1)%pressure(k)
 432             iv%instid(inst)%pf(k,n)  = Atmosphere(1)%level_pressure(k)
 433             iv%instid(inst)%tm(k,n)  = Atmosphere(1)%temperature(k)
 434             iv%instid(inst)%qm(k,n)  = Atmosphere(1)%absorber(k,1)
 435             ! T, Q Jacobian
 436             do l=1,nchanl
 437                iv%instid(inst)%t_jacobian(l,k,n) = Atmosphere_k(l,1)%temperature(k)
 438                iv%instid(inst)%q_jacobian(l,k,n) = Atmosphere_k(l,1)%absorber(k,1)
 439             end do
 440             if (crtm_cloud) then
 441                iv%instid(inst)%qcw(k,n) = Atmosphere(1)%cloud(1)%water_content(k)
 442                iv%instid(inst)%qci(k,n) = Atmosphere(1)%cloud(2)%water_content(k)
 443                iv%instid(inst)%qrn(k,n) = Atmosphere(1)%cloud(3)%water_content(k)
 444                iv%instid(inst)%qsn(k,n) = Atmosphere(1)%cloud(4)%water_content(k)
 445                iv%instid(inst)%qgr(k,n) = Atmosphere(1)%cloud(5)%water_content(k)
 446                iv%instid(inst)%qhl(k,n) = Atmosphere(1)%cloud(6)%water_content(k)
 447                iv%instid(inst)%rcw(k,n) = Atmosphere(1)%cloud(1)%effective_radius(k)
 448                iv%instid(inst)%rci(k,n) = Atmosphere(1)%cloud(2)%effective_radius(k)
 449                iv%instid(inst)%rrn(k,n) = Atmosphere(1)%cloud(3)%effective_radius(k)
 450                iv%instid(inst)%rsn(k,n) = Atmosphere(1)%cloud(4)%effective_radius(k)
 451                iv%instid(inst)%rgr(k,n) = Atmosphere(1)%cloud(5)%effective_radius(k)
 452                iv%instid(inst)%rhl(k,n) = Atmosphere(1)%cloud(6)%effective_radius(k)
 453                ! Cloud Jacobian
 454                do l=1,nchanl
 455                   iv%instid(inst)%water_jacobian(l,k,n) =   &
 456                      Atmosphere_k(l,1)%cloud(1)%water_content(k)
 457                   iv%instid(inst)%ice_jacobian(l,k,n) =     &
 458                      Atmosphere_k(l,1)%cloud(2)%water_content(k)
 459                   iv%instid(inst)%rain_jacobian(l,k,n) =    &
 460                      Atmosphere_k(l,1)%cloud(3)%water_content(k)
 461                   iv%instid(inst)%snow_jacobian(l,k,n) =    &
 462                      Atmosphere_k(l,1)%cloud(4)%water_content(k)
 463                   iv%instid(inst)%graupel_jacobian(l,k,n) = &
 464                      Atmosphere_k(l,1)%cloud(5)%water_content(k)
 465                   iv%instid(inst)%hail_jacobian(l,k,n) =    &
 466                      Atmosphere_k(l,1)%cloud(6)%water_content(k)
 467
 468                   iv%instid(inst)%water_r_jacobian(l,k,n) =   &
 469                      Atmosphere_k(l,1)%cloud(1)%effective_radius(k)
 470                   iv%instid(inst)%ice_r_jacobian(l,k,n) =     &
 471                      Atmosphere_k(l,1)%cloud(2)%effective_radius(k)
 472                   iv%instid(inst)%rain_r_jacobian(l,k,n) =    &
 473                      Atmosphere_k(l,1)%cloud(3)%effective_radius(k)
 474                   iv%instid(inst)%snow_r_jacobian(l,k,n) =    &
 475                      Atmosphere_k(l,1)%cloud(4)%effective_radius(k)
 476                   iv%instid(inst)%graupel_r_jacobian(l,k,n) = &
 477                      Atmosphere_k(l,1)%cloud(5)%effective_radius(k)
 478                   iv%instid(inst)%hail_r_jacobian(l,k,n) =    &
 479                      Atmosphere_k(l,1)%cloud(6)%effective_radius(k)
 480                end do
 481             end if
 482          end do
 483
 484          !----------------------------------------------
 485          ! [5.0] store surface information to innovation structure
 486          !----------------------------------------------
 487          iv%instid(inst)%u10(n)       = model_u10(n)
 488          iv%instid(inst)%v10(n)       = model_v10(n)
 489          iv%instid(inst)%t2m(n)       = 0.01*missing_r !model_t2m
 490          iv%instid(inst)%mr2m(n)      = 0.01*missing_r !model_mr2m
 491          iv%instid(inst)%ps(n)        = model_psfc(n)
 492          iv%instid(inst)%ts(n)        = model_ts(n)
 493          iv%instid(inst)%smois(n)     = model_smois(n)
 494          iv%instid(inst)%tslb(n)      = model_tslb(n)
 495          iv%instid(inst)%snowh(n)     = model_snowh(n)
 496          iv%instid(inst)%isflg(n)     = model_isflg
 497          iv%instid(inst)%elevation(n) = model_elv(n)
 498          iv%instid(inst)%soiltyp(n)   = model_isltyp
 499          iv%instid(inst)%vegtyp(n)    = model_ivgtyp
 500          iv%instid(inst)%vegfra(n)    = model_vegfra(n)
 501          iv%instid(inst)%clwp(n)      = clwp
 502          iv%instid(inst)%water_coverage(n) = Surface(1)%water_coverage
 503          iv%instid(inst)%land_coverage(n)  = Surface(1)%land_coverage
 504          iv%instid(inst)%ice_coverage(n)   = Surface(1)%ice_coverage
 505          iv%instid(inst)%snow_coverage(n)  = Surface(1)%snow_coverage
 506       end do       !  end loop for pixels
 507
 508       deallocate (model_u10)
 509       deallocate (model_v10)
 510       deallocate (model_psfc)
 511       deallocate (model_ts)
 512       deallocate (model_tslb)
 513       deallocate (model_snowh)
 514       deallocate (model_snow)
 515       deallocate (model_elv)
 516       deallocate (model_vegfra)
 517       deallocate (model_smois)
 518
 519       deallocate( RTSolution, RTSolution_K, STAT = Allocate_Status )
 520
 521       if (Allocate_Status /= 0) then
 522          call da_error(__FILE__,__LINE__, &
 523             (/"Error in deallocating RTSolution"/))
 524       end if
 525
 526       Error_Status = CRTM_Destroy_Surface(Surface)
 527       if (Error_Status /= 0) then
 528          call da_error(__FILE__,__LINE__, &
 529             (/"Error in deallocating CRTM Surface Structure"/))
 530       end if
 531
 532       Error_Status = CRTM_Destroy_Surface(Surface_K)
 533       if ( Error_Status /= 0 ) THEN
 534          call da_error(__FILE__,__LINE__, &
 535             (/"Error in deallocatting CRTM Surface_K Structure"/))
 536       endif
 537
 538       Error_Status = CRTM_Destroy_Atmosphere( Atmosphere_K )
 539       if ( Error_Status /= 0 ) THEN
 540          call da_error(__FILE__,__LINE__, &
 541             (/"Error in deallocatting CRTM Atmosphere_K Structure"/))
 542       endif
 543
 544       deallocate( Atmosphere_K, Surface_K, STAT = Allocate_Status )
 545       if ( Allocate_Status /= 0 ) THEN
 546          call da_error(__FILE__,__LINE__, &
 547             (/"Error in deallocatting CRTM Surface Structure"/))
 548       endif
 549
 550    end do        ! end loop for sensor
 551
 552    Error_Status = CRTM_Destroy_Atmosphere (Atmosphere)
 553    if (Error_Status /= 0) then
 554        call da_error(__FILE__,__LINE__, &
 555          (/"Error in deallocating CRTM Atmosphere Structure"/))
 556    end if
 557
 558    if (trace_use) call da_trace_exit("da_get_innov_vector_crtmk")
 559
 560 end subroutine da_get_innov_vector_crtmk
 561