var/da/da_radiance/da_get_innov_vector_crtm.inc

   1 subroutine da_get_innov_vector_crtm ( 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    !  HISTORY: 12/15/2008 effective radius unit is micron     Zhiquan Liu
  11    !---------------------------------------------------------------------------
  12
  13    implicit none
  14
  15    integer,           intent(in)    :: it       ! External iteration.
  16    type (domain),     intent(in)    :: grid     ! first guess state.
  17    type (y_type),     intent(inout) :: ob       ! Observation structure.
  18    type (iv_type),    intent(inout) :: iv       ! O-B structure.
  19
  20    integer, parameter             :: AIRS_Max_Channels = 281
  21
  22    integer :: n, icld  ! Loop counter.
  23    integer :: i, j, k  ! Index dimension.
  24    integer :: l        ! Index dimension.
  25    integer :: num_levs ! Number of obs levels.
  26    real    :: dx, dxm  ! Interpolation weights.
  27    real    :: dy, dym  ! Interpolation weights.
  28    integer :: alloc_status(40)
  29
  30    real, allocatable :: model_u10(:)
  31    real, allocatable :: model_v10(:)
  32    real, allocatable :: model_psfc(:)
  33    real              :: model_ptop
  34    real, allocatable :: model_ts(:)
  35    real, allocatable :: model_elv(:)
  36    real, allocatable :: model_smois(:)
  37    real, allocatable :: model_tslb(:)
  38    real, allocatable :: model_snowh(:)
  39    real, allocatable :: model_vegfra(:)
  40    real    :: model_isltyp, model_ivgtyp
  41    integer :: model_isflg
  42
  43    real    :: model_qcw(kms:kme)
  44    real    :: model_qci(kms:kme)
  45    real    :: model_rho(kms:kme)
  46    real, allocatable :: model_snow(:)  ! snow water equivalent, different from model_snowh,
  47                                        ! used in calculating reff_water
  48    real, allocatable :: em_mspps(:)    ! emissivity caluclated using MSPPS algorithm
  49    real              :: ts_mspps       ! surface temperature calcualted using MSPPS algorithm
  50    character(len=19) :: date_char
  51
  52    integer :: month, crtm_climat
  53
  54    ! real    :: model_tm(kms:kme)
  55
  56    integer :: inst, nchanl, n1,n2
  57    integer :: ipred, npred, gammapred
  58
  59    ! variables for computing clwp (clw/clwp) and cloud ice (ci/cip)
  60    real    :: clw(kms:kme), dpf(kms:kme), ci(kms:kme)
  61    real    :: clwp, cip
  62
  63    real    :: total_od
  64
  65    ! CRTM local varaibles and types
  66    integer :: Allocate_Status
  67    integer :: n_layers, n_absorbers, n_clouds, n_aerosols
  68 !! for ecmwf cloud
  69    real    :: t_tropp(ims:ime,jms:jme),temp1
  70    integer :: k_tropp(ims:ime,jms:jme)
  71    type (CRTM_RTSolution_type), allocatable :: RTSolution(:,:)
  72    type (CRTM_Atmosphere_type)   :: Atmosphere(1)
  73    type (CRTM_Surface_type)      :: Surface(1)
  74    type (CRTM_Geometry_type)     :: GeometryInfo(1)
  75    type (CRTM_Options_type)      :: Options(1)
  76    type (CRTM_RTSolution_type), allocatable :: RTSolution_K(:,:)
  77    type (CRTM_Atmosphere_type), allocatable :: Atmosphere_K(:,:)
  78    type (CRTM_Surface_type),    allocatable :: Surface_K(:,:)
  79
  80    real :: t(1), a(1), p(1)
  81
  82    integer, allocatable :: wrf_to_crtm_mw(:)
  83    integer :: n_vegtype
  84
  85 !! for crtm cloud
  86    real :: qcw(1),qrn(1), qci(1),qsn(1),qgr(1)
  87    logical :: calc_tb_clr
  88    type (CRTM_RTSolution_type), allocatable :: RTSolution_clr(:,:)
  89
  90    ! Initializations for Cloud Detection
  91    integer           :: ilev, jlev, nclouds
  92    real, allocatable :: hessian(:,:)
  93    real*8, allocatable :: eignvec(:,:), eignval(:)
  94    real              :: rad_clr, rad_ovc_ilev, rad_ovc_jlev
  95
  96    integer           :: Band_Size(5), Bands(AIRS_Max_Channels,5)
  97 !For Zhuge and Zou cloud detection
  98    real, allocatable :: geoht_full(:,:,:)
  99    real              :: geoht_pixel(kts:min(kte,kme-1))
 100    real              :: tt_pixel(kts:min(kte,kme-1))
 101    real              :: pp_pixel(kts:min(kte,kme-1))
 102       Band_Size(1:5) = (/86, 0, 0, 16, 0 /)
 103       Bands(:,:)     = 0
 104       Bands(1:Band_Size(1),1) = &
 105 &    (/                                                 &              !&      1,   6,   7,  10,  11,  15,  16,  17,  20,  21, &
 106 &                                                       &              !&     22,  24,  27,  28,  30,  36,  39,  40,  42,  51, &
 107 &                                                       &              !&     52,  54,  55,  56,  59,  62,  63,  68,  69,  71, &
 108 &                                                       &              !&     72,  73,  74,  75,  76,  77,  78,  79,  80,  82, &
 109 &                     92,  93,  98,  99, 101, 104, 105, &              !&     83,  84,  86,  92,  93,  98,  99, 101, 104, 105, &
 110 &     108, 110, 111, 113, 116, 117, 123, 124, 128, 129, &
 111 &     138, 139, 144, 145, 150, 151, 156, 157, 159, 162, &
 112 &     165, 168, 169, 170, 172, 173, 174, 175, 177, 179, &
 113 &     180, 182, 185, 186, 190, 192,      198, 201, 204, &              !&     180, 182, 185, 186, 190, 192, 193, 198, 201, 204, &
 114 &     207, 210,      215, 216,      221,      226, 227, &              !&     207, 210, 213, 215, 216, 218, 221, 224, 226, 227, &
 115 &     232,                     252, 253, 256, 257, 261, &              !&     232, 239, 248, 250, 251, 252, 253, 256, 257, 261, &
 116 &     262, 267, 272, 295, 299,      305,           310, &              !&     262, 267, 272, 295, 299, 300, 305, 308, 309, 310, &
 117 &          321, 325, 333, 338, 355, 362, 375, 453, 475, &              !&     318, 321, 325, 333, 338, 355, 362, 375, 453, 475, &
 118 &     484, 497, 528, 587, 672, 787, 791, 843, 870, 914, &
 119 &     950 /)
 120
 121 !      Bands(1:Band_Size(2),2) = &
 122 !&    (/ 1003, 1012, 1019, 1024, 1030, 1038, 1048, 1069, 1079, 1082,  &
 123 !&       1083, 1088, 1090, 1092, 1095, 1104, 1111, 1115, 1116, 1119,  &
 124 !&       1120, 1123, 1130, 1138, 1142, 1178, 1199, 1206, 1221, 1237,  &
 125 !&       1252, 1260, 1263, 1266, 1278, 1285 /)
 126
 127 !      Bands(1:Band_Size(3),3) = &
 128 !&    (/       1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, &  !&    1290, 1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, &
 129 !&       1466,       1477,             1500, 1519,       1538, 1545, &  !&    1466, 1471, 1477, 1479, 1488, 1500, 1519, 1520, 1538, 1545, &
 130 !&       1565, 1574, 1583, 1593,       1627, 1636,       1652, 1669, &  !&    1565, 1574, 1583, 1593, 1614, 1627, 1636, 1644, 1652, 1669, &
 131 !&                   1694, 1708,       1723, 1740, 1748,       1756, &  !&    1674, 1681, 1694, 1708, 1717, 1723, 1740, 1748, 1751, 1756, &
 132 !&             1766, 1771, 1777,       1783, 1794, 1800,       1806, &  !&    1763, 1766, 1771, 1777, 1780, 1783, 1794, 1800, 1803, 1806, &
 133 !&             1826, 1843  /)                                           !&    1812, 1826, 1843  /)
 134
 135       Bands(1:Band_Size(4),4) = &
 136 &    (/ 1852, 1865, 1866,       1868, 1869, 1872, 1873,       1876, &  !&    1852, 1865, 1866, 1867, 1868, 1869, 1872, 1873, 1875, 1876,
 137 &             1881, 1882, 1883,                   1911, 1917, 1918, &  !&    1877, 1881, 1882, 1883, 1884, 1897, 1901, 1911, 1917, 1918, &
 138 &                   1924, 1928        /)                               !&    1921, 1923, 1924, 1928, 1937  /)
 139
 140 !      Bands(1:Band_Size(5),5) = &
 141 !&    (/ 1938, 1939, 1941, 1946, 1947, 1948, 1958, 1971, 1973, 1988, &
 142 !&       1995, 2084, 2085, 2097, 2098, 2099, 2100, 2101, 2103, 2104, &
 143 !&       2106, 2107, 2108, 2109, 2110, 2111, 2112, 2113, 2114, 2115, &
 144 !&       2116, 2117, 2118, 2119, 2120, 2121, 2122, 2123, 2128, 2134, &
 145 !&       2141, 2145, 2149, 2153, 2164, 2189, 2197, 2209, 2226, 2234, &
 146 !&       2280, 2318, 2321, 2325, 2328, 2333, 2339, 2348, 2353, 2355, &
 147 !&       2363, 2370, 2371, 2377  /)
 148
 149    ! WHY? use argument it
 150    if (it==0) then; write(unit=stdout,fmt='(A)') "WHY? have argument it to"//__FILE__; end if
 151
 152    alloc_status (:) = 0
 153
 154    if (trace_use) call da_trace_entry("da_get_innov_vector_crtm")
 155
 156    ! CRTM allocation
 157
 158    ! CRTM Atmosphere
 159    n_layers    = (kte-kts)+1   ! number of vertical levels
 160    n_absorbers = 2
 161    n_aerosols  = 0
 162    n_clouds    = 0
 163    if ( crtm_cloud ) n_clouds = 6
 164
 165    call CRTM_Atmosphere_Create ( Atmosphere(1), &
 166                                  n_layers,      &
 167                                  n_absorbers,   &
 168                                  n_clouds,      &
 169                                  n_aerosols )
 170    IF ( .NOT. CRTM_Atmosphere_Associated( Atmosphere(1) ) ) THEN
 171        call da_error(__FILE__,__LINE__, &
 172          (/"Error in allocating CRTM Atmosphere Structure"/))
 173    END IF
 174
 175    Atmosphere(1)%Absorber_ID(1)=H2O_ID
 176    Atmosphere(1)%Absorber_ID(2)=O3_ID
 177    Atmosphere(1)%Absorber_Units(1) = MASS_MIXING_RATIO_UNITS
 178    Atmosphere(1)%Absorber_Units(2) = VOLUME_MIXING_RATIO_UNITS
 179
 180    if (crtm_cloud) then
 181       Atmosphere(1)%Cloud(1)%Type=WATER_CLOUD
 182       Atmosphere(1)%Cloud(2)%Type=ICE_CLOUD
 183       Atmosphere(1)%Cloud(3)%Type=RAIN_CLOUD
 184       Atmosphere(1)%Cloud(4)%Type=SNOW_CLOUD
 185       Atmosphere(1)%Cloud(5)%Type=GRAUPEL_CLOUD
 186       Atmosphere(1)%Cloud(6)%Type=HAIL_CLOUD
 187    end if
 188
 189    select case (trim(grid%mminlu))
 190       case ('USGS')
 191          n_vegtype = USGS_n_type
 192       case ('MODIFIED_IGBP_MODIS_NOAH')
 193          n_vegtype = IGBP_n_type
 194       case default
 195          call da_error(__FILE__,__LINE__,(/"Unknown land use type"/))
 196    end select
 197    allocate(wrf_to_crtm_mw(n_vegtype))
 198    if ( n_vegtype == USGS_n_type ) then
 199          wrf_to_crtm_mw = usgs_to_crtm_mw
 200    else if ( n_vegtype == IGBP_n_type ) then
 201          wrf_to_crtm_mw = igbp_to_crtm_mw
 202    end if
 203    if ( use_clddet_zz ) then
 204       allocate ( geoht_full(ims:ime,jms:jme,kms:kme-1) )
 205       do k = kms, kme-1
 206          do j = jms, jme
 207             do i = ims, ime
 208                geoht_full(i,j,k) = 0.5 * ( grid%ph_2(i,j,k)   + grid%phb(i,j,k) + &
 209                                            grid%ph_2(i,j,k+1) + grid%phb(i,j,k+1) ) / gravity
 210             end do
 211          end do
 212       end do
 213    end if
 214    !------------------------------------------------------
 215    ! [1.0] calculate the background bright temperature
 216    !-------------------------------------------------------
 217    do inst = 1, iv%num_inst                 ! loop for sensor
 218       ! if ( iv%instid(inst)%num_rad < 1 ) cycle
 219       if (iv%instid(inst)%info%n2 < iv%instid(inst)%info%n1) cycle
 220       num_levs  = kte-kts+1
 221       nchanl    = ChannelInfo(inst)%n_channels
 222
 223       if ( ChannelInfo(inst)%Sensor_Type == 2 ) then  ! INFRARED_SENSOR=2
 224          if ( index(grid%mminlu, trim(CRTM_IRlandCoeff_Classification())) <= 0 ) then
 225             call da_error(__FILE__,__LINE__, &
 226               (/"modify crtm_irland_coef to use the coef file that is consistent with your model land use type"/))
 227          end if
 228       end if
 229
 230       ! Allocate forward model solution RTSolution array to number of channels
 231       allocate (RTSolution(ChannelInfo(inst)%n_Channels,1), STAT = Allocate_Status )
 232       if ( Allocate_Status /= 0 ) then
 233          call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution"/))
 234       end if
 235
 236       call CRTM_RTSolution_Create( RTSolution , &                  !Output
 237                                    Atmosphere(1)%n_layers )        !Input
 238       if ( .NOT.(ANY(crtm_rtsolution_associated(RTSolution))) ) then
 239          call da_error(__FILE__,__LINE__, &
 240             (/"Error in allocating CRTM RTSolution Stucture"/))
 241       end if
 242
 243       calc_tb_clr = .false.
 244       if ( crtm_cloud .and. &
 245            ( trim( crtm_sensor_name(rtminit_sensor(inst))) == 'amsr2' .or. &
 246              trim( crtm_sensor_name(rtminit_sensor(inst))) == 'ahi') ) then
 247          !Tb_clear_sky is only needed for symmetric obs error model
 248          !symmetric obs error model only implemented for amsr2 for now
 249          calc_tb_clr = .true.
 250       end if
 251
 252       if ( calc_tb_clr ) then
 253          allocate (RTSolution_clr(ChannelInfo(inst)%n_Channels,1), STAT = Allocate_Status )
 254          if ( Allocate_Status /= 0 ) then
 255             call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution_clr"/))
 256          end if
 257          call CRTM_RTSolution_Create( RTSolution_clr , &              !Output
 258                                       Atmosphere(1)%n_layers )        !Input
 259          if ( .NOT.(ANY(crtm_rtsolution_associated(RTSolution_clr))) ) then
 260             call da_error(__FILE__,__LINE__, &
 261                (/"Error in allocating CRTM RTSolution_clr Stucture"/))
 262          end if
 263       end if !calc_tb_clr
 264
 265       if (use_crtm_kmatrix) then
 266          allocate (RTSolution_K(ChannelInfo(inst)%n_Channels,1), &
 267                    Atmosphere_K(ChannelInfo(inst)%n_Channels,1), &
 268                    Surface_K(ChannelInfo(inst)%n_Channels,1),    &
 269                    STAT = Allocate_Status )
 270          if ( Allocate_Status /= 0 ) then
 271             call da_error(__FILE__,__LINE__, (/"Error in allocating RTSolution_K"/))
 272          end if
 273          call CRTM_RTSolution_Create( RTSolution_K , &                !Output
 274                                       Atmosphere(1)%n_layers )        !Input
 275          if ( .NOT.(ANY(crtm_rtsolution_associated(RTSolution_K))) ) then
 276             call da_error(__FILE__,__LINE__, &
 277                (/"Error in allocating CRTM RTSolution_K Stucture"/))
 278          end if
 279       end if
 280
 281       call CRTM_Surface_Create ( Surface(1), & ! Output
 282                                  nchanl )      ! Input
 283       IF ( .NOT. CRTM_Surface_Associated( Surface(1) ) ) THEN
 284          call da_error(__FILE__,__LINE__, &
 285             (/"Error in allocating CRTM Surface Structure"/))
 286       end if
 287
 288       ! CRTM Options structure
 289       call CRTM_Options_Create( Options(1), &  ! Output
 290                                 nchanl )       ! Input
 291       IF ( .NOT. CRTM_Options_Associated( Options(1) ) ) THEN
 292         call da_error(__FILE__,__LINE__, &
 293           (/"Error in allocatting CRTM Options Structure"/))
 294       endif
 295       if ( use_antcorr(inst) ) Options(1)%Use_Antenna_Correction = .true.
 296
 297       ! do n= 1, iv%instid(inst)%num_rad           ! loop for pixel
 298
 299       n1 = iv%instid(inst)%info%n1
 300       n2 = iv%instid(inst)%info%n2
 301
 302       allocate (model_u10(n1:n2))
 303       allocate (model_v10(n1:n2))
 304       allocate (model_psfc(n1:n2))
 305       allocate (model_ts(n1:n2))
 306       allocate (model_elv(n1:n2))
 307       allocate (model_smois(n1:n2))
 308       allocate (model_tslb(n1:n2))
 309       allocate (model_snowh(n1:n2))
 310       allocate (model_snow(n1:n2))
 311       allocate (model_vegfra(n1:n2))
 312
 313       model_u10(:)    = 0.0
 314       model_v10(:)    = 0.0
 315       model_psfc(:)   = 0.0
 316       model_ts(:)     = 0.0
 317       model_elv(:)    = 0.0
 318       model_smois(:)  = 0.0
 319       model_tslb(:)   = 0.0
 320       model_snowh(:)  = 0.0
 321       model_snow(:)   = 0.0
 322       model_vegfra(:) = 0.0
 323
 324       ! Gamma correction from VarBC
 325       !----------------------------
 326 !#ifdef CRTM_MODIF
 327       if ( use_varbc .or. freeze_varbc ) then
 328          gammapred = iv%instid(inst)%varbc_info%gammapred
 329          do k = 1, nchanl
 330             npred = iv%instid(inst)%varbc(k)%npred
 331             if (npred <= 0) cycle                                           ! VarBC channels only
 332             if (iv%instid(inst)%varbc_info%npredmax < gammapred)  cycle     ! Gamma channels only
 333             if (iv%instid(inst)%varbc(k)%pred_use(gammapred) < 0) cycle     ! Gamma channels only
 334             do ipred = 1, npred
 335                if (iv%instid(inst)%varbc(k)%ipred(ipred) /= gammapred) cycle
 336                RTSolution(k,1)%Gamma = iv%instid(inst)%varbc(k)%param(ipred)
 337             end do
 338          end do
 339       end if
 340       RTSolution_K(:,1)%Gamma = 0.0
 341 !#endif
 342
 343       if ( use_mspps_emis(inst) ) allocate(em_mspps(nchanl))
 344
 345       pixel_loop: do n=n1,n2
 346          ! if ( n > iv%instid(inst)%num_rad ) exit
 347
 348          ! [1.1] Get horizontal interpolation weights:
 349
 350          i = iv%instid(inst)%info%i(1,n)
 351          j = iv%instid(inst)%info%j(1,n)
 352          dx = iv%instid(inst)%info%dx(1,n)
 353          dy = iv%instid(inst)%info%dy(1,n)
 354          dxm = iv%instid(inst)%info%dxm(1,n)
 355          dym = iv%instid(inst)%info%dym(1,n)
 356
 357          ! determine CRTM climatology based on latitude and month
 358          date_char=iv%instid(inst)%info%date_char(n)
 359          read(unit=date_char(6:7), fmt='(i2.2)') month
 360          call da_det_crtm_climat(iv%instid(inst)%info%lat(1,n), month, crtm_climat)
 361          iv%instid(inst)%crtm_climat(n) = crtm_climat
 362          Atmosphere(1)%Climatology = iv%instid(inst)%crtm_climat(n)
 363
 364          ! determine surface type of obs location
 365          !-----------------------------------------
 366          call da_detsurtyp ( grid%xb%snow, grid%xb%xice, grid%xb%landmask,  &
 367             grid%xb%ivgtyp, grid%xb%isltyp, &
 368             ims, ime, jms, jme, &
 369             i, j, dx, dy, dxm, dym, &
 370             model_isflg,model_ivgtyp, model_isltyp, &
 371             Surface(1)%Water_Coverage, Surface(1)%Ice_Coverage, &
 372             Surface(1)%Land_Coverage, Surface(1)%Snow_Coverage )
 373
 374          call da_interp_2d_partial (grid%xb%snow, iv%instid(inst)%info, 1, n, n, model_snow(n:n))
 375
 376          ! [1.2] Interpolate horizontally to ob:
 377          do k=kts,kte ! from bottom to top
 378             call da_interp_2d_partial (grid%xb%p(:,:,k), iv%instid(inst)%info,k,n,n,p)
 379             call da_interp_2d_partial (grid%xb%t(:,:,k), iv%instid(inst)%info,k,n,n,t)
 380             call da_interp_2d_partial (grid%xb%q(:,:,k), iv%instid(inst)%info,k,n,n,a)
 381
 382             Atmosphere(1)%Pressure(kte-k+1)    = 0.01   * p(1)  ! convert Pa to hPa
 383             Atmosphere(1)%Temperature(kte-k+1) = t(1)
 384             Atmosphere(1)%Absorber(kte-k+1,1)  = 1000.0 * a(1)/(1.0-a(1))  ! in g/kg mixing ratio
 385
 386             ! NOTE: WRF high-level q values seems too big, replaced by constants
 387             if (p(1)*0.01 < 75.0) Atmosphere(1)%Absorber(kte-k+1,1) = 0.001
 388
 389             call da_interp_2d_partial (grid%xb%qcw(:,:,k), iv%instid(inst)%info,k,n,n, model_qcw(kte-k+1:kte-k+1))
 390
 391             if (crtm_cloud) then
 392
 393                call da_interp_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n,qci)
 394                call da_interp_2d_partial (grid%xb%qci(:,:,k), iv%instid(inst)%info,k,n,n,model_qci(kte-k+1:kte-k+1))
 395
 396                call da_interp_2d_partial (grid%xb%qrn(:,:,k), iv%instid(inst)%info,k,n,n,qrn)
 397
 398                call da_interp_2d_partial (grid%xb%qsn(:,:,k), iv%instid(inst)%info, k,n,n,qsn)
 399
 400                call da_interp_2d_partial (grid%xb%qgr(:,:,k), iv%instid(inst)%info, k,n,n,qgr)
 401
 402                Atmosphere(1)%Cloud(1)%Water_Content(kte-k+1)=model_qcw(kte-k+1)
 403                Atmosphere(1)%Cloud(2)%Water_Content(kte-k+1)=qci(1)
 404                Atmosphere(1)%Cloud(3)%Water_Content(kte-k+1)=qrn(1)
 405                Atmosphere(1)%Cloud(4)%Water_Content(kte-k+1)=qsn(1)
 406                Atmosphere(1)%Cloud(5)%Water_Content(kte-k+1)=qgr(1)
 407                Atmosphere(1)%Cloud(6)%Water_Content(kte-k+1)=0.0
 408
 409                call da_interp_2d_partial (grid%xb%rho(:,:,k), iv%instid(inst)%info, k,n,n, &
 410                   model_rho(k:k) )
 411
 412                call da_cld_eff_radius(Atmosphere(1)%Temperature(kte-k+1),model_rho(k),&
 413                                       Atmosphere(1)%Cloud(2)%Water_Content(kte-k+1),  &  !qci
 414                                       Atmosphere(1)%Cloud(3)%Water_Content(kte-k+1),  &  !qrn
 415                                       Atmosphere(1)%Cloud(4)%Water_Content(kte-k+1),  &  !qsn
 416                                       Atmosphere(1)%Cloud(5)%Water_Content(kte-k+1),  &  !qgr
 417                                       model_snow(n),                                  &
 418                                       Surface(1)%Ice_Coverage, Surface(1)%Land_Coverage, 1, &
 419                                       Atmosphere(1)%Cloud(1)%Effective_Radius(kte-k+1), &
 420                                       Atmosphere(1)%Cloud(2)%Effective_Radius(kte-k+1), &
 421                                       Atmosphere(1)%Cloud(3)%Effective_Radius(kte-k+1), &
 422                                       Atmosphere(1)%Cloud(4)%Effective_Radius(kte-k+1), &
 423                                       Atmosphere(1)%Cloud(5)%Effective_Radius(kte-k+1) )
 424
 425                ! reset the da_cld_eff_radius calcualted effective radius to constants if desired
 426
 427                Atmosphere(1)%Cloud(1)%Effective_Radius(kte-k+1)=10.0  ! in micron
 428                Atmosphere(1)%Cloud(2)%Effective_Radius(kte-k+1)=30.0  ! in micron
 429                !Atmosphere(1)%Cloud(3)%Effective_Radius(kte-k+1)=300
 430                !Atmosphere(1)%Cloud(4)%Effective_Radius(kte-k+1)=600
 431                !Atmosphere(1)%Cloud(5)%Effective_Radius(kte-k+1)=600
 432                Atmosphere(1)%Cloud(6)%Effective_Radius(kte-k+1)=600
 433
 434             end if
 435          end do
 436          if (use_clddet_zz) then
 437             ! Find tropopause temperature for Zhuge and Zou Cloud Detection
 438             do k = kts, min(kte,kme-1)
 439                tt_pixel(k) = Atmosphere(1)%Temperature(kte-k+1)
 440                pp_pixel(k) = Atmosphere(1)%Pressure(kte-k+1)
 441                call da_interp_2d_partial ( geoht_full(:,:,k), iv%instid(inst)%info, k, n, n, geoht_pixel(k) )
 442             end do
 443             call da_trop_wmo ( tt_pixel, geoht_pixel, pp_pixel, (min(kte,kme-1)-kts+1), tropt = iv%instid(inst)%tropt(n) )
 444          end if
 445
 446
 447          call da_interp_2d_partial (grid%xb%u10,  iv%instid(inst)%info, 1, n, n, model_u10(n:n))
 448          call da_interp_2d_partial (grid%xb%v10,  iv%instid(inst)%info, 1, n, n, model_v10(n:n))
 449          call da_interp_2d_partial (grid%xb%psfc, iv%instid(inst)%info, 1, n, n, model_psfc(n:n))
 450
 451          model_psfc(n) = 0.01*model_psfc(n)           ! convert to hPa
 452          model_ptop    = 0.01*grid%xb%ptop
 453
 454          ! get CRTM levels (0.005hPa at top) /model full level
 455          Atmosphere(1)%Level_Pressure(0)                      = model_ptop    ! to sigma level 51-->sigmaf=0
 456          Atmosphere(1)%Level_Pressure(Atmosphere(1)%n_Layers) = model_psfc(n) ! to sigma level 1->sigmaf=1
 457
 458          ! calculating full-level P from half-level P using
 459          ! linear interpolation: (pf-pf1)/(sigmaf-sigmaf1)=(ph2-ph1)/(sigmah2-sigmah1)
 460          ! note: sigma is from bottem to top; pressure is from top to bottem
 461          do k=kts,kte-1 ! from top to bottem
 462             Atmosphere(1)%Level_Pressure(k)= (Atmosphere(1)%Pressure(k+1) - Atmosphere(1)%Pressure(k))/ &
 463                                              (grid%xb%sigmah(kte-k) - grid%xb%sigmah(kte-k+1)) * &
 464                                              (grid%xb%sigmaf(kte-k+1) - grid%xb%sigmah(kte-k+1)) + &
 465                                               Atmosphere(1)%Pressure(k)
 466          end do
 467
 468          ! check for pressure monotonicity
 469          do k = kte, kts+1, -1  ! from bottom to top
 470             if ( Atmosphere(1)%Level_Pressure(k) <= Atmosphere(1)%Level_Pressure(k-1) ) then
 471                write(message(1),'(a,2i5.0,a)') ' Skipping the pixel at loc ', i, j,  &
 472                   ' where the pressure is not monotonic'
 473                call da_warning(__FILE__,__LINE__,message(1:1))
 474                iv%instid(inst)%tb_inv(:,n) = missing_r
 475                iv%instid(inst)%info%proc_domain(:,n) = .false.
 476                cycle pixel_loop
 477             end if
 478          end do
 479
 480          ! convert cloud content unit from kg/kg to kg/m^2
 481          if (crtm_cloud) then
 482             do k=kts,kte
 483                do icld=1,Atmosphere(1)%n_Clouds
 484                   Atmosphere(1)%Cloud(icld)%Water_Content(k)= Atmosphere(1)%Cloud(icld)%Water_Content(k)* &
 485                      (Atmosphere(1)%Level_Pressure(k)- Atmosphere(1)%Level_Pressure(k-1))*100.0/gravity
 486                end do
 487             end do
 488          end if
 489
 490          if ( model_isflg == 0 ) then   ! over sea using SST
 491             call da_interp_2d_partial (grid%xb % tgrn, iv%instid(inst)%info, 1, n, n, model_ts(n:n))
 492          else
 493             call da_interp_2d_partial (grid%xb % tsk,  iv%instid(inst)%info, 1, n, n, model_ts(n:n))
 494          end if
 495
 496          call da_interp_2d_partial (grid%xb % terr, iv%instid(inst)%info, 1, n, n, model_elv(n:n))
 497
 498          ! variables for emissivity calculations
 499          !----------------------------------------
 500          call da_interp_2d_partial (grid%xb % smois,  iv%instid(inst)%info, 1, n, n, model_smois(n:n) )
 501          call da_interp_2d_partial (grid%xb % tslb,   iv%instid(inst)%info, 1, n, n, model_tslb(n:n) )
 502          call da_interp_2d_partial (grid%xb % snowh,  iv%instid(inst)%info, 1, n, n, model_snowh(n:n) )
 503          call da_interp_2d_partial (grid%xb % vegfra, iv%instid(inst)%info, 1, n, n, model_vegfra(n:n) )
 504
 505          ! model_snowh(n) = model_snowh(n)*100.0   ! convert from m to mm
 506          model_vegfra(n) = 0.01*model_vegfra(n)  ! convert range to 0~1
 507
 508          ! ADD for computing cloud liquid water path (mm) from guess
 509          clwp = 0.0
 510          do k = kts,kte ! from top to bottom
 511             dpf(k) = 100.0*(Atmosphere(1)%level_pressure(k) - Atmosphere(1)%level_pressure(k-1))
 512             clw  (k) = model_qcw(k)*dpf(k)/gravity ! kg/m2 or mm
 513             if (Atmosphere(1)%pressure(k)<100.0) clw(k) = 0.0
 514             clwp  = clwp + clw(k)
 515          end do
 516
 517          ! computing ice path from guess
 518          if (crtm_cloud) then
 519             cip = 0.0
 520             do k = kts,kte ! from top to bottom
 521                dpf(k) = 100.0*(Atmosphere(1)%level_pressure(k) - Atmosphere(1)%level_pressure(k-1))
 522                ci(k) = model_qci(k)*dpf(k)/gravity ! kg/m2 or mm
 523                cip  = cip + ci(k)
 524             end do
 525          endif
 526
 527          ! CRTM GeometryInfo Structure
 528          GeometryInfo(1)%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
 529          GeometryInfo(1)%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
 530          GeometryInfo(1)%iFOV=iv%instid(inst)%scanpos(n)
 531          ! GeometryInfo(1)%Satellite_Height=830.0
 532          ! GeometryInfo(1)%Sensor_Scan_Angle=
 533          ! GeometryInfo(1)%Sensor_Zenith_Angle=
 534          ! GeometryInfo(1)%Sensor_Scan_Angle=
 535          ! GeometryInfo(1)%Source_Zenith_Angle=
 536
 537          ! CRTM Surface parameter data
 538
 539          if ( use_mspps_ts(inst) ) then
 540             ! only for AMSU-A over land
 541             if ( trim(crtm_sensor_name(rtminit_sensor(inst))) == 'amsua' &
 542                  .and. model_isflg == 2 ) then
 543                call da_mspps_ts(ob%instid(inst)%tb(1:nchanl,n), nchanl,  &
 544                                 iv%instid(inst)%satzen(n), ts_mspps)
 545                ! ts_mspps is initilaized as negative values in the
 546                ! da_mspps_ts subroutine.  Apply only valid values here.
 547                if ( ts_mspps > 0.0 ) then
 548                   model_ts(n) = ts_mspps
 549                end if
 550             end if
 551          end if
 552
 553          if (Surface(1)%Land_Coverage > 0.0) then
 554             Surface(1)%Land_Type       = min(11,nint(model_ivgtyp))
 555             Surface(1)%Vegetation_Type = max(1,wrf_to_crtm_mw(nint(model_ivgtyp)))
 556             Surface(1)%Soil_Type       = max(1,wrf_to_crtm_soil(nint(model_isltyp)))
 557             ! glacial land ice soil type and glacial vegetation type are not valid land types
 558             if ( surface(1)%Soil_Type == 9 .or. surface(1)%Vegetation_Type == 13 ) then
 559                surface(1)%ice_coverage = min(surface(1)%ice_coverage + surface(1)%land_coverage, 1.0)
 560                surface(1)%land_coverage = 0
 561             end if
 562             ! reset water-dominated land type
 563             if ( Surface(1)%Land_Type == grid%iswater ) Surface(1)%Land_Type = 1
 564             Surface(1)%Land_Temperature=model_ts(n)      ! K
 565             Surface(1)%Soil_Moisture_Content= model_smois(n) !0.05    ! volumetric water content (g/cm**3)
 566             ! Surface(1)%Canopy_Water_Content=0.05      ! gravimetric water content
 567             Surface(1)%Vegetation_Fraction=model_vegfra(n)
 568             Surface(1)%Soil_Temperature=model_tslb(n)
 569          end if
 570          if (Surface(1)%Water_Coverage > 0.0) then
 571             ! Surface%Water_Type=SEA_WATER          ! (Currently NOT used)
 572             Surface(1)%Water_Temperature=model_ts(n)     ! K
 573             Surface(1)%Wind_Speed=sqrt(model_u10(n)**2+model_v10(n)**2)  ! m/sec
 574             ! surface(1)%Wind_Direction=0.0            ! NOT used
 575             Surface(1)%Salinity=33.0                   ! ppmv
 576          end if
 577
 578          if (Surface(1)%Snow_Coverage > 0.0) then
 579             Surface(1)%Snow_Type=2 !NEW_SNOW
 580             Surface(1)%Snow_Temperature=model_ts(n)      ! K
 581             Surface(1)%Snow_Depth=model_snowh(n)         ! mm
 582             ! Surface(1)%Snow_Density=0.2               ! g/cm**3
 583             ! Surface(1)%Snow_Grain_Size=2.0            ! mm
 584          end if
 585          if (Surface(1)%Ice_Coverage > 0.0) then
 586             ! Surface(1)%Ice_Type=FRESH_ICE             ! NO Table offered, single example is FRESH_ICE
 587             Surface(1)%Ice_Temperature=model_ts(n)       ! K
 588             Surface(1)%Ice_Thickness=10.0              ! mm
 589             ! Surface(1)%Ice_Density=0.9                ! g/cm**3
 590             ! Surface(1)%Ice_Roughness=0.0               ! NO Table offered, single example is ZERO
 591
 592          end if
 593          if (nchanl > 0) then
 594             Surface(1)%SensorData%n_channels        = nchanl
 595             Surface(1)%SensorData%Sensor_Channel    = ChannelInfo(inst)%Sensor_Channel
 596             Surface(1)%SensorData%Sensor_Id         = ChannelInfo(inst)%Sensor_Id
 597             Surface(1)%SensorData%WMO_Satellite_Id  = ChannelInfo(inst)%WMO_Satellite_Id
 598             Surface(1)%SensorData%WMO_Sensor_Id     = ChannelInfo(inst)%WMO_Sensor_Id
 599             Surface(1)%SensorData%Tb(1:nchanl)      = ob%instid(inst)%tb(1:nchanl,n)
 600          end if
 601
 602          ! user emissivity options
 603          Options(1)%use_emissivity = .false.
 604          if ( use_mspps_emis(inst) ) then
 605             ! Only for AMSU-A over land
 606             if ( trim(crtm_sensor_name(rtminit_sensor(inst))) == 'amsua' &
 607                  .and. model_isflg == 2 ) then
 608                call da_mspps_emis(ob%instid(inst)%tb(1:nchanl,n), nchanl, em_mspps)
 609                Options(1)%use_emissivity = .true.
 610                Options(1)%emissivity(1:nchanl) = em_mspps(1:nchanl)
 611             end if
 612          end if
 613
 614          if (use_crtm_kmatrix) then
 615
 616             ! CRTM surface/atmosphere K initialization
 617             do l = 1, ChannelInfo(inst)%n_Channels
 618                ! -- Copy the adjoint atmosphere structure
 619                Atmosphere_K(l,1) = Atmosphere(1)
 620
 621                ! -- Copy the adjoint surface structure
 622                Surface_K(l,1) = Surface(1)
 623             end do
 624
 625             ! -- Zero the Adjoint outputs
 626             ! Important: adjoint variables must be initialized
 627             call CRTM_Atmosphere_Zero( Atmosphere_K )
 628             call CRTM_Surface_Zero( Surface_K )
 629
 630             ! Assign tb = R^-1 Re :
 631             RTSolution_K(:,1)%brightness_temperature = 1.
 632             RTSolution_K(:,1)%radiance = 0.
 633
 634             ! [1.3] Call RTM K-Matrix model
 635             call da_crtm_k(1, nchanl, 1, Atmosphere,   &
 636                                Surface,      &
 637                                RTSolution_K,&
 638                                GeometryInfo, &
 639                                ChannelInfo(inst),  &
 640                                Atmosphere_K,&
 641                                Surface_K,   &
 642                                RTSolution,  &
 643                                Options)
 644          else
 645
 646             ! [1.3] Call RTM forward model
 647             call da_crtm_direct(1, nchanl, 1, Atmosphere,   &
 648                Surface,      &
 649                GeometryInfo, &
 650                ChannelInfo(inst:inst),  &
 651                RTSolution,              &
 652                Options)
 653
 654          end if
 655
 656 ! Compute Overcast Radiances for MMR in Auligné (2014).or. PF in Xu et al., (2016).or. MW03 ecmwf method in McNally AP and Watts PD (2003)
 657          !----------------------------------------------------------------
 658          if (use_clddet==1 .or. use_clddet==2.or. use_clddet==3 ) then
 659             do i = 1, nchanl
 660                !crtm_2.2.3 contains Upwelling_Overcast_Radiance
 661                iv%instid(inst)%rad_ovc(i,kts:kte,n)=RTSolution(i,1)%Upwelling_Overcast_Radiance(:)
 662             end do
 663          end if
 664
 665    !---------------------------------------------------------------------------
 666    ! Precondition satellite control variable (Cloud Cover(s)):
 667    !---------------------------------------------------------------------------
 668          if (.false.) then
 669 !         if (use_satcv(2)) then
 670             nclouds = iv%instid(inst)%cv_index(n)%nclouds
 671             allocate (hessian(nclouds,nclouds))
 672             allocate (eignvec(nclouds,nclouds))
 673             allocate (eignval(nclouds))
 674             hessian(:,:)= 0.0
 675             do ilev=1, nclouds
 676                do jlev=ilev, nclouds
 677                   do k = 1, nchanl
 678                      if (ALL(iv%instid(inst)%ichan(k) /=  Bands(:,1))) cycle   ! Only Channels in Band 1
 679                         rad_clr      = iv%instid(inst)%rad_xb(k,n)
 680                         rad_ovc_ilev = iv%instid(inst)%rad_ovc(k,kte-nclouds+ilev,n)
 681                         rad_ovc_jlev = iv%instid(inst)%rad_ovc(k,kte-nclouds+jlev,n)
 682                         hessian(ilev,jlev) = hessian(ilev,jlev) + (rad_ovc_ilev-rad_clr) * (rad_ovc_jlev-rad_clr)
 683                   end do
 684                  hessian(jlev,ilev) = hessian(ilev,jlev)
 685                end do
 686             end do
 687
 688             call da_eof_decomposition(nclouds, hessian, eignvec, eignval)
 689
 690 !           if (ANY(eignval <= 0)) write(unit=stdout,fmt='(3A,I4,A,100F18.5)')      &
 691 !              'SATCV: non-positive Hessian for ', trim(iv%instid(inst)%rttovid_string), ' ,pixel ',n,'--> Eigenvalues =',eignval
 692
 693             do ilev = 1, nclouds
 694                do jlev = ilev, nclouds
 695                   iv%instid(inst)%cv_index(n)%vtox(ilev,jlev) = &
 696                      sum( eignvec(ilev,:) * sqrt(1.0/eignval(:)) * eignvec(jlev,:), mask = eignval >0 )
 697
 698                   iv%instid(inst)%cv_index(n)%vtox(jlev,ilev) = iv%instid(inst)%cv_index(n)%vtox(ilev,jlev)
 699                end do
 700             end do
 701
 702             deallocate(hessian,eignvec,eignval)
 703          end if
 704
 705          !----------------------------------------------------------------
 706          ! [2.0] calculate components of innovation vector:
 707          !----------------------------------------------------------------
 708          do k = 1, nchanl
 709             iv%instid(inst)%tb_xb(k,n)  = RTSolution(k,1)%Brightness_Temperature
 710             iv%instid(inst)%rad_xb(k,n) = RTSolution(k,1)%Radiance
 711             iv%instid(inst)%emiss(k,n)  = RTSolution(k,1)%surface_emissivity
 712
 713             if (use_pseudo_rad .or. use_simulated_rad) then ! input is innovation
 714               ob%instid(inst)%tb(k,n) = missing_r
 715               if ( iv%instid(inst)%tb_inv(k,n) > missing_r ) &
 716                 ob%instid(inst)%tb(k,n) = RTSolution(k,1)%Brightness_Temperature + iv%instid(inst)%tb_inv(k,n)
 717             else
 718                if ( iv%instid(inst)%tb_inv(k,n) > missing_r ) &
 719                   iv%instid(inst)%tb_inv(k,n) = ob%instid(inst)%tb(k,n) - iv%instid(inst)%tb_xb(k,n)
 720             end if
 721
 722             if (use_crtm_kmatrix) then
 723                ! surface Jacobian
 724                iv%instid(inst)%ts_jacobian(k,n) = Surface_k(k,1)%water_temperature
 725                iv%instid(inst)%windspeed_jacobian(k,n) = Surface_k(k,1)%wind_speed
 726                iv%instid(inst)%emiss_jacobian(k,n) = RTSolution_k(k,1)%surface_emissivity
 727 !#ifdef CRTM_MODIF
 728                iv%instid(inst)%gamma_jacobian(k,n) = RTSolution_k(k,1)%Gamma
 729 !#endif
 730             end if
 731          end do
 732
 733          !----------------------------------------------------------------
 734          ! [3.0] store base state (and Jacobian) to innovation structure
 735          !----------------------------------------------------------------
 736          ! full level pressures
 737          iv%instid(inst)%pf(0,n)  = Atmosphere(1)%level_pressure(0)
 738          if (use_crtm_kmatrix) then
 739             ! PS Jacobian
 740             do l=1,nchanl
 741                iv%instid(inst)%ps_jacobian(l,n)  = Atmosphere_k(l,1)%level_pressure(Atmosphere(1)%n_layers)
 742             end do
 743          end if
 744          do k=1,Atmosphere(1)%n_layers
 745             iv%instid(inst)%pm(k,n)  = Atmosphere(1)%pressure(k)
 746             iv%instid(inst)%pf(k,n)  = Atmosphere(1)%level_pressure(k)
 747             iv%instid(inst)%tm(k,n)  = Atmosphere(1)%temperature(k)
 748             iv%instid(inst)%qm(k,n)  = Atmosphere(1)%absorber(k,1)
 749
 750             if (use_crtm_kmatrix) then
 751                ! T, Q Jacobian
 752                do l=1,nchanl
 753                   iv%instid(inst)%t_jacobian(l,k,n) = Atmosphere_k(l,1)%temperature(k)
 754                   iv%instid(inst)%q_jacobian(l,k,n) = Atmosphere_k(l,1)%absorber(k,1)
 755                end do
 756             end if
 757          end do
 758
 759          if (crtm_cloud) then
 760             do k=1,Atmosphere(1)%n_layers
 761                iv%instid(inst)%qcw(k,n) = Atmosphere(1)%cloud(1)%water_content(k)
 762                iv%instid(inst)%qci(k,n) = Atmosphere(1)%cloud(2)%water_content(k)
 763                iv%instid(inst)%qrn(k,n) = Atmosphere(1)%cloud(3)%water_content(k)
 764                iv%instid(inst)%qsn(k,n) = Atmosphere(1)%cloud(4)%water_content(k)
 765                iv%instid(inst)%qgr(k,n) = Atmosphere(1)%cloud(5)%water_content(k)
 766                iv%instid(inst)%qhl(k,n) = Atmosphere(1)%cloud(6)%water_content(k)
 767                iv%instid(inst)%rcw(k,n) = Atmosphere(1)%cloud(1)%effective_radius(k)
 768                iv%instid(inst)%rci(k,n) = Atmosphere(1)%cloud(2)%effective_radius(k)
 769                iv%instid(inst)%rrn(k,n) = Atmosphere(1)%cloud(3)%effective_radius(k)
 770                iv%instid(inst)%rsn(k,n) = Atmosphere(1)%cloud(4)%effective_radius(k)
 771                iv%instid(inst)%rgr(k,n) = Atmosphere(1)%cloud(5)%effective_radius(k)
 772                iv%instid(inst)%rhl(k,n) = Atmosphere(1)%cloud(6)%effective_radius(k)
 773
 774                if (use_crtm_kmatrix) then
 775                   ! Cloud Jacobian
 776                   do l=1,nchanl
 777                      iv%instid(inst)%water_jacobian(l,k,n)    = Atmosphere_k(l,1)%cloud(1)%water_content(k)
 778                      iv%instid(inst)%ice_jacobian(l,k,n)      = Atmosphere_k(l,1)%cloud(2)%water_content(k)
 779                      iv%instid(inst)%rain_jacobian(l,k,n)     = Atmosphere_k(l,1)%cloud(3)%water_content(k)
 780                      iv%instid(inst)%snow_jacobian(l,k,n)     = Atmosphere_k(l,1)%cloud(4)%water_content(k)
 781                      iv%instid(inst)%graupel_jacobian(l,k,n)  = Atmosphere_k(l,1)%cloud(5)%water_content(k)
 782                      iv%instid(inst)%hail_jacobian(l,k,n)     = Atmosphere_k(l,1)%cloud(6)%water_content(k)
 783
 784                      iv%instid(inst)%water_r_jacobian(l,k,n)  = Atmosphere_k(l,1)%cloud(1)%effective_radius(k)
 785                      iv%instid(inst)%ice_r_jacobian(l,k,n)    = Atmosphere_k(l,1)%cloud(2)%effective_radius(k)
 786                      iv%instid(inst)%rain_r_jacobian(l,k,n)   = Atmosphere_k(l,1)%cloud(3)%effective_radius(k)
 787                      iv%instid(inst)%snow_r_jacobian(l,k,n)   = Atmosphere_k(l,1)%cloud(4)%effective_radius(k)
 788                      iv%instid(inst)%graupel_r_jacobian(l,k,n)= Atmosphere_k(l,1)%cloud(5)%effective_radius(k)
 789                      iv%instid(inst)%hail_r_jacobian(l,k,n)   = Atmosphere_k(l,1)%cloud(6)%effective_radius(k)
 790                   end do
 791                end if
 792             end do
 793          end if
 794
 795          if ( calc_weightfunc .and. use_crtm_kmatrix ) then
 796             do l = 1, nchanl
 797                total_od = 0.0
 798                do k = 1, Atmosphere(1)%n_layers
 799                   iv%instid(inst)%lod(l,k,n) = RTSolution(l,1)%layer_optical_depth(k)
 800                   iv%instid(inst)%lod_jacobian(l,k,n) = RTSolution_K(l,1)%layer_optical_depth(k)
 801                   total_od = total_od + RTSolution(l,1)%layer_optical_depth(k)
 802                   iv%instid(inst)%trans(l,k,n) = exp(-total_od)
 803                   iv%instid(inst)%trans_jacobian(l,k,n) = -iv%instid(inst)%lod_jacobian(l,k,n)*exp(iv%instid(inst)%lod(l,k,n))
 804                end do
 805                iv%instid(inst)%der_trans(l,1,n) = 0.0
 806                do k= 2, Atmosphere(1)%n_layers
 807                   iv%instid(inst)%der_trans(l,k,n) =  abs( (iv%instid(inst)%trans(l,k-1,n)-iv%instid(inst)%trans(l,k,n))/ &
 808                                                            (log(iv%instid(inst)%pm(k-1,n))-log(iv%instid(inst)%pm(k,n))) )
 809                end do
 810             end do
 811          end if
 812
 813          if (use_clddet==3) then  !begin ecmwf cloud detection
 814             !kmin_t,kmax_p corresponding tropopause and KBPL
 815             ! find the k index of tropopause
 816             t_tropp = 999.0  ! initialize
 817             k_tropp = kte    ! initialize the k index of tropopause
 818             do k = kte, kts, -1
 819                do j = jts, jte
 820                   do i = its, ite
 821                      if ( grid%xb%t(i,j,k) < t_tropp(i,j) .and.  &
 822                           grid%xb%p(i,j,k) > 5000.0 ) then  ! tropopause should not be higher than 50 hPa
 823                          t_tropp(i,j) = grid%xb%t(i,j,k)
 824                          k_tropp(i,j) = k
 825                          iv%instid(inst)%kmin_t(n)=Atmosphere(1)%n_layers-k_tropp(i,j)+1
 826                      end if
 827                   end do
 828                end do
 829             end do
 830             iv%instid(inst)%kmax_p(n)=Atmosphere(1)%n_layers-grid%kpbl(i,j)+1
 831             !write(unit=stdout,fmt='(A)') &
 832             !   'preparing ECMWF cloud detection input using cloud sensitivity to rank channels'
 833             do l=1,nchanl
 834                temp1=0
 835                do k=1,Atmosphere(1)%n_layers  !rad_ovc(l,1,n) top
 836                   iv%instid(inst)%sensitivity_ratio(l,k,n) = abs(iv%instid(inst)%rad_xb(l,n)-iv%instid(inst)%rad_ovc(l,k,n))
 837                   iv%instid(inst)%sensitivity_ratio(l,k,n) = iv%instid(inst)%sensitivity_ratio(l,k,n)/iv%instid(inst)%rad_xb(l,n)
 838                end do
 839             end do
 840             do l=1,nchanl
 841                temp1=0
 842                do k=Atmosphere(1)%n_layers,1,-1 ! sensitivity_ratio(l,1,n) top
 843                   if( iv%instid(inst)%sensitivity_ratio(l,k,n) > 0.01 ) then
 844                      iv%instid(inst)%p_chan_level(l,n) = real(k) ! the peak level
 845                      exit
 846                   end if
 847                end do
 848             end do
 849          end if     !end ecmwf cloud detection
 850
 851          !----------------------------------------------
 852          ! [4.0] store surface information to innovation structure
 853          !----------------------------------------------
 854          iv%instid(inst)%u10(n)       = model_u10(n)
 855          iv%instid(inst)%v10(n)       = model_v10(n)
 856          iv%instid(inst)%t2m(n)       = 0.01*missing_r !model_t2m
 857          iv%instid(inst)%mr2m(n)      = 0.01*missing_r !model_mr2m
 858          iv%instid(inst)%ps(n)        = model_psfc(n)
 859          iv%instid(inst)%ts(n)        = model_ts(n)
 860          iv%instid(inst)%smois(n)     = model_smois(n)
 861          iv%instid(inst)%tslb(n)      = model_tslb(n)
 862          iv%instid(inst)%snowh(n)     = model_snowh(n)
 863          iv%instid(inst)%isflg(n)     = model_isflg
 864          iv%instid(inst)%elevation(n) = model_elv(n)
 865          iv%instid(inst)%soiltyp(n)   = model_isltyp
 866          iv%instid(inst)%vegtyp(n)    = model_ivgtyp
 867          iv%instid(inst)%vegfra(n)    = model_vegfra(n)
 868          iv%instid(inst)%clwp(n)      = clwp
 869          if (crtm_cloud) iv%instid(inst)%cip(n) = cip
 870          iv%instid(inst)%water_coverage(n) = Surface(1)%water_coverage
 871          iv%instid(inst)%land_coverage(n)  = Surface(1)%land_coverage
 872          iv%instid(inst)%ice_coverage(n)   = Surface(1)%ice_coverage
 873          iv%instid(inst)%snow_coverage(n)  = Surface(1)%snow_coverage
 874
 875          if ( calc_tb_clr ) then
 876             atmosphere(1)%n_clouds = 0
 877             call da_crtm_direct(1, nchanl, 1, Atmosphere,   &
 878                Surface,                                     &
 879                GeometryInfo,                                &
 880                ChannelInfo(inst:inst),                      &
 881                RTSolution_clr,                              &
 882                Options)
 883             do k = 1, nchanl
 884                iv%instid(inst)%tb_xb_clr(k,n) = RTSolution_clr(k,1)%Brightness_Temperature
 885             end do
 886             ! resetting n_clouds
 887             atmosphere(1)%n_clouds = n_clouds
 888          end if !calc_tb_clr
 889
 890       end do pixel_loop      !  end loop for pixels
 891
 892       if ( use_mspps_emis(inst) ) deallocate(em_mspps)
 893
 894       deallocate (model_u10)
 895       deallocate (model_v10)
 896       deallocate (model_psfc)
 897       deallocate (model_ts)
 898       deallocate (model_tslb)
 899       deallocate (model_snowh)
 900       deallocate (model_snow)
 901       deallocate (model_elv)
 902       deallocate (model_vegfra)
 903       deallocate (model_smois)
 904
 905       call CRTM_RTSolution_Destroy(RTSolution)
 906       deallocate( RTSolution, STAT = Allocate_Status )
 907       if (Allocate_Status /= 0) &
 908          call da_error(__FILE__,__LINE__, &
 909             (/"Error in deallocating RTSolution"/))
 910
 911       if ( calc_tb_clr ) then
 912          call CRTM_RTSolution_Destroy(RTSolution_clr)
 913          deallocate( RTSolution_clr, STAT = Allocate_Status )
 914          if (Allocate_Status /= 0) &
 915             call da_error(__FILE__,__LINE__, &
 916                (/"Error in deallocating RTSolution_clr"/))
 917       end if !calc_tb_clr
 918
 919       if (use_crtm_kmatrix) then
 920          call CRTM_RTSolution_Destroy(RTSolution_K)
 921          deallocate( RTSolution_K, STAT = Allocate_Status )
 922          if (Allocate_Status /= 0) &
 923             call da_error(__FILE__,__LINE__, &
 924                (/"Error in deallocating RTSolution_K"/))
 925       end if
 926
 927       call CRTM_Options_Destroy(Options)
 928       IF ( CRTM_Options_Associated( Options(1) ) ) &
 929          call da_error(__FILE__,__LINE__, &
 930             (/"Error in deallocating CRTM Options Structure"/))
 931
 932       call CRTM_Surface_Destroy(Surface)
 933       IF ( CRTM_Surface_Associated( Surface(1) ) ) &
 934          call da_error(__FILE__,__LINE__, &
 935             (/"Error in deallocating CRTM Surface Structure"/))
 936
 937       if (use_crtm_kmatrix) then
 938          call CRTM_Surface_Destroy(Surface_K)
 939          IF ( ANY( CRTM_Surface_Associated( Surface_K(1,:)) ) ) &
 940             call da_error(__FILE__,__LINE__, &
 941                (/"Error in deallocatting CRTM Surface_K Structure"/))
 942       endif
 943
 944       if (use_crtm_kmatrix) then
 945          call CRTM_Atmosphere_Destroy( Atmosphere_K )
 946          IF ( ANY( CRTM_Atmosphere_Associated( Atmosphere_K(1,:)) ) ) &
 947             call da_error(__FILE__,__LINE__, &
 948                (/"Error in deallocatting CRTM Atmosphere_K Structure"/))
 949       endif
 950
 951       if (use_crtm_kmatrix) then
 952          deallocate( Atmosphere_K, Surface_K, STAT = Allocate_Status )
 953          if ( Allocate_Status /= 0 ) &
 954             call da_error(__FILE__,__LINE__, &
 955                (/"Error in deallocatting CRTM Surface_K Structure"/))
 956       endif
 957
 958    end do        ! end loop for sensor
 959
 960    deallocate (wrf_to_crtm_mw)
 961    if ( use_clddet_zz ) deallocate ( geoht_full )
 962    call CRTM_Atmosphere_Destroy (Atmosphere)
 963    IF ( CRTM_Atmosphere_Associated( Atmosphere(1) ) ) &
 964        call da_error(__FILE__,__LINE__, &
 965          (/"Error in deallocating CRTM Atmosphere Structure"/))
 966
 967    if (trace_use) call da_trace_exit("da_get_innov_vector_crtm")
 968
 969 end subroutine da_get_innov_vector_crtm
 970