var/da/da_radiance/da_transform_xtoy_crtm.inc

   1 subroutine da_transform_xtoy_crtm (cv_size, cv, grid, iv, y )
   2
   3    !---------------------------------------------------------------------------
   4    !  PURPOSE: transform from analysis increment to
   5    !                          pertubation radiance.
   6    !
   7    !  METHOD:  delta_y = H delta_x
   8    !           1. input reference state of CRTM_TL
   9    !           2. interpolate analysis increment to obs location
  10    !           3. Call CRTM_TL
  11    !
  12    !  HISTORY: 11/16/2006 - Creation                        Zhiquan Liu
  13    !           11/25/2008 - Zero Jacobian for top level     Tom Auligne
  14    !
  15    !---------------------------------------------------------------------------
  16
  17    implicit none
  18
  19    integer, intent(in)            :: cv_size         ! Size of cv array.
  20    real, intent(in)               :: cv(1:cv_size)   ! control variables.
  21    type (domain),  intent(in)     :: grid
  22    type (y_type),  intent(inout)  :: y        ! H' delta_x
  23    type (iv_type), intent(in)     :: iv       ! O-B structure.
  24
  25    integer, parameter             :: AIRS_Max_Channels = 281
  26
  27 #ifdef CRTM
  28    integer                        :: k, l  ! Index dimension.
  29
  30    integer           :: inst, num_rad, nchanl, n, icld
  31    integer           :: ipred, npred, gammapred, id
  32    real, allocatable :: temperature(:,:)
  33    real, allocatable :: absorber(:,:), xb_q(:,:)
  34    real, allocatable :: psfc(:)
  35 !! for crtm_cloud
  36    real, allocatable :: qcw(:,:),qci(:,:),qrn(:,:),qsn(:,:),qgr(:,:)
  37
  38    ! CRTM local variables and types
  39    integer :: Allocate_Status
  40    integer :: n_layers, n_absorbers, n_clouds, n_aerosols
  41    type( CRTM_RTSolution_type ), ALLOCATABLE :: RTSolution(:,:),RTSolution_TL(:,:)
  42    type( CRTM_Atmosphere_type ), allocatable :: Atmosphere(:), Atmosphere_TL(:)
  43    type( CRTM_Surface_type ),    allocatable :: Surface(:), Surface_TL(:)
  44    type( CRTM_Geometry_type ),   allocatable :: GeometryInfo(:)
  45    type( CRTM_Options_type ), allocatable    :: Options(:)
  46
  47    integer                        :: ts_index
  48    integer                        :: nclouds, ncv
  49    real, allocatable              :: cc_tl(:)
  50    real*8                         :: rad_clr, rad_cld   ! RT clear/cloudy radiances
  51    real, allocatable              :: rad_ovc(:)         ! RT overcast radiances
  52    real*8                         :: rad_tl, tb_tl
  53
  54    ! Initializations for AIRS (MMR) Cloud Detection
  55    integer                        :: Band_Size(5), Bands(AIRS_Max_Channels,5)
  56
  57       Band_Size(1:5) = (/86, 0, 0, 16, 0 /)
  58       Bands(:,:)     = 0
  59       Bands(1:Band_Size(1),1) = &
  60 &    (/                                                 &              !&      1,   6,   7,  10,  11,  15,  16,  17,  20,  21, &
  61 &                                                       &              !&     22,  24,  27,  28,  30,  36,  39,  40,  42,  51, &
  62 &                                                       &              !&     52,  54,  55,  56,  59,  62,  63,  68,  69,  71, &
  63 &                                                       &              !&     72,  73,  74,  75,  76,  77,  78,  79,  80,  82, &
  64 &                     92,  93,  98,  99, 101, 104, 105, &              !&     83,  84,  86,  92,  93,  98,  99, 101, 104, 105, &
  65 &     108, 110, 111, 113, 116, 117, 123, 124, 128, 129, &
  66 &     138, 139, 144, 145, 150, 151, 156, 157, 159, 162, &
  67 &     165, 168, 169, 170, 172, 173, 174, 175, 177, 179, &
  68 &     180, 182, 185, 186, 190, 192,      198, 201, 204, &              !&     180, 182, 185, 186, 190, 192, 193, 198, 201, 204, &
  69 &     207, 210,      215, 216,      221,      226, 227, &              !&     207, 210, 213, 215, 216, 218, 221, 224, 226, 227, &
  70 &     232,                     252, 253, 256, 257, 261, &              !&     232, 239, 248, 250, 251, 252, 253, 256, 257, 261, &
  71 &     262, 267, 272, 295, 299,      305,           310, &              !&     262, 267, 272, 295, 299, 300, 305, 308, 309, 310, &
  72 &          321, 325, 333, 338, 355, 362, 375, 453, 475, &              !&     318, 321, 325, 333, 338, 355, 362, 375, 453, 475, &
  73 &     484, 497, 528, 587, 672, 787, 791, 843, 870, 914, &
  74 &     950 /)
  75
  76 !      Bands(1:Band_Size(2),2) = &
  77 !&    (/ 1003, 1012, 1019, 1024, 1030, 1038, 1048, 1069, 1079, 1082,  &
  78 !&       1083, 1088, 1090, 1092, 1095, 1104, 1111, 1115, 1116, 1119,  &
  79 !&       1120, 1123, 1130, 1138, 1142, 1178, 1199, 1206, 1221, 1237,  &
  80 !&       1252, 1260, 1263, 1266, 1278, 1285 /)
  81
  82 !      Bands(1:Band_Size(3),3) = &
  83 !&    (/       1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, &  !&    1290, 1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, &
  84 !&       1466,       1477,             1500, 1519,       1538, 1545, &  !&    1466, 1471, 1477, 1479, 1488, 1500, 1519, 1520, 1538, 1545, &
  85 !&       1565, 1574, 1583, 1593,       1627, 1636,       1652, 1669, &  !&    1565, 1574, 1583, 1593, 1614, 1627, 1636, 1644, 1652, 1669, &
  86 !&                   1694, 1708,       1723, 1740, 1748,       1756, &  !&    1674, 1681, 1694, 1708, 1717, 1723, 1740, 1748, 1751, 1756, &
  87 !&             1766, 1771, 1777,       1783, 1794, 1800,       1806, &  !&    1763, 1766, 1771, 1777, 1780, 1783, 1794, 1800, 1803, 1806, &
  88 !&             1826, 1843  /)                                           !&    1812, 1826, 1843  /)
  89
  90       Bands(1:Band_Size(4),4) = &
  91 &    (/ 1852, 1865, 1866,       1868, 1869, 1872, 1873,       1876, &  !&    1852, 1865, 1866, 1867, 1868, 1869, 1872, 1873, 1875, 1876,
  92 &             1881, 1882, 1883,                   1911, 1917, 1918, &  !&    1877, 1881, 1882, 1883, 1884, 1897, 1901, 1911, 1917, 1918, &
  93 &                   1924, 1928        /)                               !&    1921, 1923, 1924, 1928, 1937  /)
  94
  95 !      Bands(1:Band_Size(5),5) = &
  96 !&    (/ 1938, 1939, 1941, 1946, 1947, 1948, 1958, 1971, 1973, 1988, &
  97 !&       1995, 2084, 2085, 2097, 2098, 2099, 2100, 2101, 2103, 2104, &
  98 !&       2106, 2107, 2108, 2109, 2110, 2111, 2112, 2113, 2114, 2115, &
  99 !&       2116, 2117, 2118, 2119, 2120, 2121, 2122, 2123, 2128, 2134, &
 100 !&       2141, 2145, 2149, 2153, 2164, 2189, 2197, 2209, 2226, 2234, &
 101 !&       2280, 2318, 2321, 2325, 2328, 2333, 2339, 2348, 2353, 2355, &
 102 !&       2363, 2370, 2371, 2377  /)
 103
 104 !---------------------------------------------------------
 105
 106    if ( iv%num_inst < 1 ) return
 107
 108    if (trace_use) call da_trace_entry("da_transform_xtoy_crtm")
 109
 110    sensor_loop : do inst = 1, iv%num_inst                 ! loop for sensor
 111
 112       num_rad = iv%instid(inst)%info%n2 - iv%instid(inst)%info%n1 + 1
 113       if ( num_rad < 1 ) cycle
 114
 115       allocate (Atmosphere   (iv%instid(inst)%info%n1:iv%instid(inst)%info%n2), &
 116                 Atmosphere_TL(iv%instid(inst)%info%n1:iv%instid(inst)%info%n2), &
 117                 STAT = Allocate_Status )
 118       if ( Allocate_Status /= 0 ) then
 119          call da_error(__FILE__,__LINE__, &
 120           (/"Error in allocatting Atmosphere"/))
 121       end if
 122       allocate (Surface      (iv%instid(inst)%info%n1:iv%instid(inst)%info%n2), &
 123                 Surface_TL   (iv%instid(inst)%info%n1:iv%instid(inst)%info%n2), &
 124                 STAT = Allocate_Status )
 125       if ( Allocate_Status /= 0 ) then
 126          call da_error(__FILE__,__LINE__, &
 127           (/"Error in allocatting Surface"/))
 128       end if
 129       allocate (GeometryInfo (iv%instid(inst)%info%n1:iv%instid(inst)%info%n2), &
 130                 STAT = Allocate_Status )
 131       if ( Allocate_Status /= 0 ) then
 132          call da_error(__FILE__,__LINE__, &
 133           (/"Error in allocatting GeometryInfo"/))
 134       end if
 135       allocate (Options      (iv%instid(inst)%info%n1:iv%instid(inst)%info%n2), &
 136                 STAT = Allocate_Status )
 137       if ( Allocate_Status /= 0 ) then
 138          call da_error(__FILE__,__LINE__, &
 139           (/"Error in allocatting Options"/))
 140       end if
 141
 142 !----------------------------------------------------------------------------
 143 ! CRTM allocation
 144 !
 145 ! Atmosphere structure
 146
 147       do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2 ! loop for pixel
 148
 149          n_layers    = (kte-kts)+1   ! number of vertical levels
 150          n_absorbers = 2
 151          n_aerosols  = 0
 152          n_clouds    = 0
 153          if ( crtm_cloud ) n_clouds = 6
 154
 155          call CRTM_Atmosphere_Create( Atmosphere(n), &
 156                                       n_layers,      &
 157                                       n_absorbers,   &
 158                                       n_clouds,      &
 159                                       n_aerosols )
 160          IF ( .NOT. CRTM_Atmosphere_Associated( Atmosphere(n) ) ) THEN
 161             call da_error(__FILE__,__LINE__, &
 162                          (/"Error in allocatting CRTM Atmosphere Structure"/))
 163          end if
 164
 165          Atmosphere(n)%Absorber_ID(1)=H2O_ID
 166          Atmosphere(n)%Absorber_ID(2)=O3_ID
 167          Atmosphere(n)%Absorber_Units(1) = MASS_MIXING_RATIO_UNITS
 168          Atmosphere(n)%Absorber_Units(2) = VOLUME_MIXING_RATIO_UNITS
 169          Atmosphere(n)%Climatology=iv%instid(inst)%crtm_climat(n)
 170
 171          if (crtm_cloud) then
 172             Atmosphere(n)%Cloud(1)%Type=WATER_CLOUD
 173             Atmosphere(n)%Cloud(2)%Type=ICE_CLOUD
 174             Atmosphere(n)%Cloud(3)%Type=RAIN_CLOUD
 175             Atmosphere(n)%Cloud(4)%Type=SNOW_CLOUD
 176             Atmosphere(n)%Cloud(5)%Type=GRAUPEL_CLOUD
 177             Atmosphere(n)%Cloud(6)%Type=HAIL_CLOUD
 178          end if
 179
 180       end do
 181
 182 !-------------------------------------------------------------------------------
 183
 184       nchanl    = ChannelInfo(inst)%n_channels
 185
 186   ! Allocate forward model solution RTSolution array to number of channels
 187       allocate( RTSolution( ChannelInfo(inst)%n_Channels, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2 ), &
 188                 RTSolution_TL( ChannelInfo(inst)%n_Channels, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2 ), &
 189                 STAT = Allocate_Status )
 190       if ( Allocate_Status /= 0 ) then
 191          call da_error(__FILE__,__LINE__, &
 192                       (/"Error in allocatting RTSolution"/))
 193       END IF
 194
 195       do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2 ! loop for pixel
 196
 197          call CRTM_Surface_Create( Surface(n),  & ! Output
 198                                    nchanl )       ! Input
 199          IF ( .NOT. CRTM_Surface_Associated( Surface(n) ) ) THEN
 200             call da_error(__FILE__,__LINE__, &
 201                          (/"Error in allocatting CRTM Surface Structure"/))
 202          end if
 203
 204   ! CRTM Options structure
 205          call CRTM_Options_Create( Options(n),  &  ! Output
 206                                    nchanl )        ! Input
 207          IF ( .NOT. CRTM_Options_Associated( Options(n) ) ) THEN
 208             call da_error(__FILE__,__LINE__, &
 209                          (/"Error in allocatting CRTM Options Structure"/))
 210          endif
 211          if ( use_antcorr(inst) ) Options(n)%Use_Antenna_Correction = .true.
 212
 213       end do
 214
 215       allocate (temperature(Atmosphere(iv%instid(inst)%info%n1)%n_layers, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 216       allocate (absorber(Atmosphere(iv%instid(inst)%info%n1)%n_layers, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 217       allocate (xb_q(Atmosphere(iv%instid(inst)%info%n1)%n_layers, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 218       allocate (psfc(iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 219
 220       temperature(:,:) = 0.0
 221       absorber(:,:)    = 0.0
 222       xb_q(:,:)        = 0.0
 223       psfc(:)          = 0.0
 224
 225       if (crtm_cloud) then
 226
 227          allocate (qcw(Atmosphere(iv%instid(inst)%info%n1)%n_layers, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 228          allocate (qci(Atmosphere(iv%instid(inst)%info%n1)%n_layers, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 229          allocate (qrn(Atmosphere(iv%instid(inst)%info%n1)%n_layers, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 230          allocate (qsn(Atmosphere(iv%instid(inst)%info%n1)%n_layers, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 231          allocate (qgr(Atmosphere(iv%instid(inst)%info%n1)%n_layers, iv%instid(inst)%info%n1:iv%instid(inst)%info%n2))
 232          qcw(:,:) = 0.0
 233          qci(:,:) = 0.0
 234          qrn(:,:) = 0.0
 235          qsn(:,:) = 0.0
 236          qgr(:,:) = 0.0
 237
 238       end if
 239
 240       do k=kts,kte ! from bottom to top
 241          call da_interp_2d_partial (grid%xa%t(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%info%n1, iv%instid(inst)%info%n2, &
 242             temperature(kte-k+1,:))
 243
 244          ! dq=xa%q and absorber both specific humidity (kg h2o / kg moist air)
 245          call da_interp_2d_partial (grid%xa%q(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%info%n1, iv%instid(inst)%info%n2, &
 246             absorber(kte-k+1,:))
 247
 248          if ( crtm_cloud .and. cloud_cv_options > 0 ) then
 249
 250             call da_interp_2d_partial (grid%xa%qcw(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%info%n1, iv%instid(inst)%info%n2, &
 251                qcw(kte-k+1,:))
 252             call da_interp_2d_partial (grid%xa%qrn(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%info%n1, iv%instid(inst)%info%n2, &
 253                qrn(kte-k+1,:))
 254             if ( cloud_cv_options > 1 ) then
 255                call da_interp_2d_partial (grid%xa%qci(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%info%n1, iv%instid(inst)%info%n2, &
 256                   qci(kte-k+1,:))
 257                call da_interp_2d_partial (grid%xa%qsn(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%info%n1, iv%instid(inst)%info%n2, &
 258                   qsn(kte-k+1,:))
 259                call da_interp_2d_partial (grid%xa%qgr(:,:,k), iv%instid(inst)%info, k, iv%instid(inst)%info%n1, iv%instid(inst)%info%n2, &
 260                   qgr(kte-k+1,:))
 261             end if
 262
 263          end if
 264
 265       end do
 266
 267       call da_interp_2d_partial (grid%xa%psfc, iv%instid(inst)%info, 1, iv%instid(inst)%info%n1, iv%instid(inst)%info%n2, psfc(:))
 268
 269       ! Gamma correction from VarBC
 270       !----------------------------
 271 !#ifdef CRTM_MODIF
 272       if (use_varbc) then
 273         gammapred = iv%instid(inst)%varbc_info%gammapred
 274         do k = 1, nchanl
 275            npred = iv%instid(inst)%varbc(k)%npred
 276            if (npred <= 0) cycle                                               ! VarBC channels only
 277            if (iv%instid(inst)%varbc_info%npredmax < gammapred)  cycle         ! Gamma channels only
 278            if (iv%instid(inst)%varbc(k)%pred_use(gammapred) < 0) cycle         ! Gamma channels only
 279            do ipred = 1, npred
 280               if (iv%instid(inst)%varbc(k)%ipred(ipred) /= gammapred) cycle
 281               id = iv%instid(inst)%varbc(k)%index(ipred)
 282               RTSolution(k,:)%Gamma    = iv%instid(inst)%varbc(k)%param(ipred)
 283               RTSolution_TL(k,:)%Gamma = SUM(cv(id) * iv%instid(inst)%varbc(k)%vtox(ipred,1:npred))
 284            end do
 285         end do
 286       end if
 287 !#endif
 288
 289       do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2 ! loop for pixel
 290
 291       ! [1.1] Get horizontal interpolation weights:
 292
 293       ! [1.3] Extract base state Atmosphere variables
 294          Atmosphere(n)%level_pressure(0) = iv%instid(inst)%pf(0,n)
 295          do k=1,Atmosphere(n)%n_layers
 296             Atmosphere(n)%pressure(k)       = iv%instid(inst)%pm(k,n)
 297             Atmosphere(n)%level_pressure(k) = iv%instid(inst)%pf(k,n)
 298             Atmosphere(n)%temperature(k)    = iv%instid(inst)%tm(k,n)
 299             Atmosphere(n)%absorber(k,1)     = iv%instid(inst)%qm(k,n)
 300             ! convert mixing ratio (g h2o / kg dry air) to specific humidity (kg h2o / kg moist air)
 301             xb_q(k,n) = 0.001*iv%instid(inst)%qm(k,n)/(1.0+0.001*iv%instid(inst)%qm(k,n)) ! specific humidity
 302          end do
 303          if (crtm_cloud) then
 304             do k=1,Atmosphere(n)%n_layers
 305               Atmosphere(n)%cloud(1)%water_content(k)=iv%instid(inst)%qcw(k,n)
 306               Atmosphere(n)%cloud(2)%water_content(k)=iv%instid(inst)%qci(k,n)
 307               Atmosphere(n)%cloud(3)%water_content(k)=iv%instid(inst)%qrn(k,n)
 308               Atmosphere(n)%cloud(4)%water_content(k)=iv%instid(inst)%qsn(k,n)
 309               Atmosphere(n)%cloud(5)%water_content(k)=iv%instid(inst)%qgr(k,n)
 310               Atmosphere(n)%cloud(6)%water_content(k)=iv%instid(inst)%qhl(k,n)
 311               Atmosphere(n)%cloud(1)%effective_radius(k)=iv%instid(inst)%rcw(k,n)
 312               Atmosphere(n)%cloud(2)%effective_radius(k)=iv%instid(inst)%rci(k,n)
 313               Atmosphere(n)%cloud(3)%effective_radius(k)=iv%instid(inst)%rrn(k,n)
 314               Atmosphere(n)%cloud(4)%effective_radius(k)=iv%instid(inst)%rsn(k,n)
 315               Atmosphere(n)%cloud(5)%effective_radius(k)=iv%instid(inst)%rgr(k,n)
 316               Atmosphere(n)%cloud(6)%effective_radius(k)=iv%instid(inst)%rhl(k,n)
 317            end do
 318          end if
 319
 320       ! [1.4] User-supplied emissivity
 321            !Options%emissivity_switch = 1
 322            !Options%emissivity(1:Options%n_channels) = &
 323            !    iv%instid(inst)%emiss(1:Options%n_channels,n)
 324
 325       ! [1.4] CRTM Surface parameter data
 326          Surface(n)%Land_Coverage=iv%instid(inst)%land_coverage(n)
 327          Surface(n)%Water_Coverage=iv%instid(inst)%water_coverage(n)
 328          Surface(n)%Snow_Coverage=iv%instid(inst)%snow_coverage(n)
 329          Surface(n)%Ice_Coverage=iv%instid(inst)%ice_coverage(n)
 330
 331          if (Surface(n)%Land_Coverage > 0.0) then
 332             Surface(n)%Land_Type=11 !GRASS_SOIL
 333             Surface(n)%Land_Temperature=iv%instid(inst)%ts(n)      ! K
 334             Surface(n)%Soil_Moisture_Content=iv%instid(inst)%smois(n)  !0.05    ! volumetric water content (g/cm**3)
 335             !Surface(n)%Canopy_Water_Content=0.05      ! gravimetric water content
 336             Surface(n)%Vegetation_Fraction=iv%instid(inst)%vegfra(n)
 337             Surface(n)%Soil_Temperature=iv%instid(inst)%tslb(n)
 338          end if
 339          if (Surface(n)%Water_Coverage > 0.0) then
 340             !Surface(n)%Water_Type=SEA_WATER          ! (Currently NOT used)
 341             Surface(n)%Water_Temperature=iv%instid(inst)%ts(n)     ! K
 342             Surface(n)%Wind_Speed=sqrt((iv%instid(inst)%u10(n))**2+ &
 343                                       (iv%instid(inst)%v10(n))**2)  ! m/sec
 344             !surface(n)%Wind_Direction=0.0            ! NOT used
 345             Surface(n)%Salinity=33.                   ! ppmv
 346          end if
 347          if (Surface(n)%Snow_Coverage > 0.0) then
 348             Surface(n)%Snow_Type=2 !NEW_SNOW
 349             Surface(n)%Snow_Temperature=iv%instid(inst)%ts(n)      ! K
 350             Surface(n)%Snow_Depth=iv%instid(inst)%snowh(n)         ! mm
 351             !Surface(n)%Snow_Density=0.2               ! g/cm**3
 352             !Surface(n)%Snow_Grain_Size=2.0            ! mm
 353          end if
 354          if (Surface(n)%Ice_Coverage > 0.0) then
 355             !Surface(n)%Ice_Type=FRESH_ICE             ! NO Table offered, single egrid%xample is FRESH_ICE
 356             Surface(n)%Ice_Temperature=iv%instid(inst)%ts(n)       ! K
 357             Surface(n)%Ice_Thickness=10.0              ! mm
 358             !Surface(n)%Ice_Density=0.9                ! g/cm**3
 359             !Surface(n)%Ice_Roughness=0.0               ! NO Table offered, single egrid%xample is ZERO
 360          end if
 361          Surface(n)%SensorData%n_channels        = nchanl
 362          Surface(n)%SensorData%Sensor_Channel    = ChannelInfo(inst)%Sensor_Channel
 363          Surface(n)%SensorData%Sensor_Id         = ChannelInfo(inst)%Sensor_Id
 364          Surface(n)%SensorData%WMO_Satellite_Id  = ChannelInfo(inst)%WMO_Satellite_Id
 365          Surface(n)%SensorData%WMO_Sensor_Id     = ChannelInfo(inst)%WMO_Sensor_Id
 366          Surface(n)%SensorData%Tb(1:nchanl) = iv%instid(inst)%tb_inv(1:nchanl,n) + &
 367                                               iv%instid(inst)%tb_xb(1:nchanl,n)
 368
 369       ! -- Copy the TL atmosphere structure
 370          Atmosphere_TL(n) = Atmosphere(n)
 371
 372       ! -- Copy the TL surface structure
 373          Surface_TL(n) = Surface(n)
 374
 375     ! -- Zero the TL inputs
 376     ! Important: adjoint variables must be initialized
 377         call CRTM_Atmosphere_Zero( Atmosphere_TL(n) )
 378         call CRTM_Surface_Zero( Surface_TL(n) )
 379
 380         ! Humidity in g/kg. Zero Jacobian for top level(s) (above 75hPa)
 381         do k = kts, kte
 382           if ( iv%instid(inst)%pm(k,n) < 75.0 ) absorber(k,n) = 0.0
 383         end do
 384         ! compute mixing ratio increment, dr, (g h2o / kg dry air) from xb_q and dq (both kg h2o / kg moist air)
 385         Atmosphere_TL(n)%Absorber(kts+1:kte,1)  = 1000.0 / (1.0-xb_q(kts+1:kte,n))**2  &
 386                                                   * absorber(kts+1:kte,n)
 387
 388         Atmosphere_TL(n)%Temperature(kts+1:kte) = temperature(kts+1:kte,n)       ! Zero Jacobian for top level
 389         Atmosphere_TL(n)%Level_Pressure(Atmosphere_TL(n)%n_Layers) = 0.01 * psfc(n)
 390
 391         if (crtm_cloud) then
 392            Atmosphere_TL(n)%Cloud(1)%Water_Content(kts:kte)             = qcw(kts:kte,n)
 393            Atmosphere_TL(n)%Cloud(2)%Water_Content(kts:kte)             = qci(kts:kte,n)
 394            Atmosphere_TL(n)%Cloud(3)%Water_Content(kts:kte)             = qrn(kts:kte,n)
 395            Atmosphere_TL(n)%Cloud(4)%Water_Content(kts:kte)             = qsn(kts:kte,n)
 396            Atmosphere_TL(n)%Cloud(5)%Water_Content(kts:kte)             = qgr(kts:kte,n)
 397            Atmosphere_TL(n)%Cloud(6)%Water_Content(kts:kte)             = 0.
 398         ! convert cloud content unit from kg/kg to kg/m^2
 399            do k=kts,kte
 400               do icld=1,Atmosphere(n)%n_Clouds
 401
 402                  Atmosphere_TL(n)%Cloud(icld)%Water_Content(k)=  Atmosphere_TL(n)%Cloud(icld)%Water_Content(k)* &
 403                                 (Atmosphere(n)%Level_Pressure(k)- Atmosphere(n)%Level_Pressure(k-1))*100./gravity
 404
 405               enddo
 406            enddo
 407         end if
 408
 409       ! Skin Temperature
 410       !-----------------
 411         if (use_satcv(1)) then
 412            ts_index = iv%instid(inst)%cv_index(n)%ts
 413            if (Surface(n)%Land_Coverage  > 0.0) Surface_TL(n)%Land_Temperature  = cv(ts_index)
 414            if (Surface(n)%Water_Coverage > 0.0) Surface_TL(n)%Water_Temperature = cv(ts_index)
 415            if (Surface(n)%Snow_Coverage  > 0.0) Surface_TL(n)%Snow_Temperature  = cv(ts_index)
 416            if (Surface(n)%Ice_Coverage   > 0.0) Surface_TL(n)%Ice_Temperature   = cv(ts_index)
 417         end if
 418
 419       ! Cloud cover(s)
 420       !---------------
 421         if (use_satcv(2)) then
 422            nclouds = iv%instid(inst)%cv_index(n)%nclouds
 423            ncv     = iv%instid(inst)%cv_index(n)%ncv
 424            allocate (rad_ovc (nclouds))
 425            allocate (cc_tl   (nclouds))
 426
 427            cc_tl(:) = cv(iv%instid(inst)%cv_index(n)%cc(:))
 428 !           !---------------------------------------------------------------
 429 !           ! Change of variable (preconditioning)
 430 !           !---------------------------------------------------------------
 431 !          do icld = 1, nclouds
 432 !             cc_tl(icld) = SUM( cv(iv%instid(inst)%cv_index(n)%cc) * &
 433 !                                   iv%instid(inst)%cv_index(n)%vtox(icld,1:ncv) )
 434 !          end do
 435
 436            do k = 1, nchanl
 437               if (ALL(iv%instid(inst)%ichan(k) /=  Bands(:,1))) cycle   ! Only Channels in Band 1
 438               rad_clr    = iv%instid(inst)%rad_xb(k,n)
 439               rad_ovc(:) = iv%instid(inst)%rad_ovc(k,kte-nclouds+1:kte,n)
 440               rad_cld    = SUM(cc_tl*rad_ovc) + (1.0 - SUM(cc_tl))*rad_clr
 441               rad_tl     = SUM(cc_tl*(rad_ovc-rad_clr))
 442               call CRTM_Planck_Temperature_TL(inst,k,rad_clr,rad_tl,tb_tl)
 443               y%instid(inst)%tb(k,n) = y%instid(inst)%tb(k,n) + tb_tl
 444            end do
 445
 446            deallocate(cc_tl, rad_ovc)
 447         else
 448            y%instid(inst)%tb(:,n) = 0.0
 449         end if
 450
 451       ! [1.5] CRTM GeometryInfo Structure
 452         GeometryInfo(n)%Sensor_Zenith_Angle=iv%instid(inst)%satzen(n)
 453         GeometryInfo(n)%Source_Zenith_Angle=iv%instid(inst)%solzen(n)
 454         GeometryInfo(n)%iFOV=iv%instid(inst)%scanpos(n)
 455    !    GeometryInfo(n)%Satellite_Height=830.0
 456    !    GeometryInfo(n)%Sensor_Scan_Angle=
 457    !    GeometryInfo(n)%Sensor_Zenith_Angle=
 458    !    GeometryInfo(n)%Sensor_Scan_Angle=
 459    !    GeometryInfo(n)%Source_Zenith_Angle=
 460
 461       end do
 462
 463
 464       ! [1.6] Call CRTM_TL model
 465
 466       !$OMP PARALLEL DO &
 467       !$OMP PRIVATE ( n )
 468       do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2 ! loop for pixel
 469          if (.not. use_crtm_kmatrix) then
 470             call da_crtm_tl (1, nchanl, 1, Atmosphere(n),   &
 471                                Surface(n),      &
 472                                Atmosphere_TL(n),&
 473                                Surface_TL(n),   &
 474                                GeometryInfo(n), &
 475                                ChannelInfo(inst:inst),  &
 476                                RTSolution(:,n),   &
 477                                RTSolution_TL(:,n),&
 478                                Options(n))
 479          else
 480             RTSolution_TL(:,n)%brightness_temperature = 0.
 481             do l = 1, ChannelInfo(inst)%n_Channels
 482                do k = kts , kte
 483                   RTSolution_TL(l,n)%brightness_temperature = RTSolution_TL(l,n)%brightness_temperature + &
 484                                        iv%instid(inst)%t_jacobian(l,k,n) * Atmosphere_TL(n)%Temperature(k) +  &
 485                                        iv%instid(inst)%q_jacobian(l,k,n) * Atmosphere_TL(n)%absorber(k,1)
 486                   if (crtm_cloud) then
 487                      RTSolution_TL(l,n)%brightness_temperature = RTSolution_TL(l,n)%brightness_temperature + &
 488                         iv%instid(inst)%water_jacobian(l,k,n)   * Atmosphere_TL(n)%Cloud(1)%Water_Content(k) + &
 489                         iv%instid(inst)%ice_jacobian(l,k,n)     * Atmosphere_TL(n)%Cloud(2)%Water_Content(k) + &
 490                         iv%instid(inst)%rain_jacobian(l,k,n)    * Atmosphere_TL(n)%Cloud(3)%Water_Content(k) + &
 491                         iv%instid(inst)%snow_jacobian(l,k,n)    * Atmosphere_TL(n)%Cloud(4)%Water_Content(k) + &
 492                         iv%instid(inst)%graupel_jacobian(l,k,n) * Atmosphere_TL(n)%Cloud(5)%Water_Content(k) + &
 493                         iv%instid(inst)%hail_jacobian(l,k,n)    * Atmosphere_TL(n)%Cloud(6)%Water_Content(k)
 494                   end if
 495                end do
 496                RTSolution_TL(l,n)%brightness_temperature = RTSolution_TL(l,n)%brightness_temperature + &
 497                                        iv%instid(inst)%ps_jacobian(l,n) * Atmosphere_TL(n)%Level_Pressure(Atmosphere_TL(n)%n_Layers)
 498             end do
 499          endif
 500       end do
 501      !$OMP END PARALLEL DO
 502
 503       !-------------------------------------------------------------------
 504       ! [1.7] assign Hdx :
 505       !-------------------------------------------------------------------
 506       do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2 ! loop for pixel
 507
 508          y%instid(inst)%tb(:,n) = y%instid(inst)%tb(:,n) + RTSolution_TL(:,n)%brightness_temperature
 509
 510          call CRTM_Atmosphere_Destroy( Atmosphere_TL(n) )
 511          IF ( CRTM_Atmosphere_Associated(Atmosphere_TL(n)) ) THEN
 512             call da_error(__FILE__,__LINE__, &
 513                (/"Error in deallocatting CRTM Atmosphere_TL Structure"/))
 514          end if
 515
 516          call CRTM_Surface_Destroy(Surface_TL(n))
 517          IF ( CRTM_Surface_Associated(Surface_TL(n)) ) THEN
 518             call da_error(__FILE__,__LINE__, &
 519                (/"Error in deallocatting CRTM Surface_TL Structure"/))
 520          end if
 521
 522       end do  ! end loop for pixels
 523
 524       deallocate (temperature)
 525       deallocate (absorber)
 526       deallocate (xb_q)
 527       deallocate (psfc)
 528
 529       if (crtm_cloud) then
 530          deallocate (qcw)
 531          deallocate (qci)
 532          deallocate (qrn)
 533          deallocate (qsn)
 534          deallocate (qgr)
 535       end if
 536
 537       do n = iv%instid(inst)%info%n1, iv%instid(inst)%info%n2 ! loop for pixel
 538       !-------------------------------------------------------------------
 539       ! [2.0] Deallocating CRTM structures
 540       ! -------------------------------------------------------------------
 541          call CRTM_Options_Destroy(Options(n))
 542          IF ( CRTM_Options_Associated(Options(n)) ) THEN
 543             call da_error(__FILE__,__LINE__, &
 544                          (/"Error in deallocatting CRTM Options Structure"/))
 545          end if
 546
 547          call CRTM_Surface_Destroy(Surface(n))
 548          IF ( CRTM_Surface_Associated(Surface(n)) ) THEN
 549             call da_error(__FILE__,__LINE__, &
 550                          (/"Error in deallocatting CRTM Surface Structure"/))
 551          end if
 552       !-------------------------------------------------------------------
 553       ! [3.0] Deallocating CRTM Atmosphere structures
 554       !-------------------------------------------------------------------
 555          call CRTM_Atmosphere_Destroy( Atmosphere(n) )
 556          IF ( CRTM_Atmosphere_Associated(Atmosphere(n)) ) THEN
 557             call da_error(__FILE__,__LINE__, &
 558                          (/"Error in deallocatting CRTM Atmosphere Structure"/))
 559          end if
 560
 561       end do
 562
 563       deallocate( RTSolution, RTSolution_TL, STAT = Allocate_Status )
 564       if ( Allocate_Status /= 0 ) then
 565          call da_error(__FILE__,__LINE__, &
 566                       (/"Error in deallocatting RTSolution"/))
 567       END IF
 568
 569       deallocate( Atmosphere, Atmosphere_TL, STAT = Allocate_Status )
 570       if ( Allocate_Status /= 0 ) then
 571          call da_error(__FILE__,__LINE__, &
 572                       (/"Error in deallocatting Atmosphere"/))
 573       END IF
 574
 575       deallocate( Surface, Surface_TL, STAT = Allocate_Status )
 576       if ( Allocate_Status /= 0 ) then
 577          call da_error(__FILE__,__LINE__, &
 578                       (/"Error in deallocatting Surface"/))
 579       END IF
 580
 581       deallocate( Options, STAT = Allocate_Status )
 582       if ( Allocate_Status /= 0 ) then
 583          call da_error(__FILE__,__LINE__, &
 584                       (/"Error in deallocatting Options"/))
 585       END IF
 586
 587       deallocate( GeometryInfo, STAT = Allocate_Status )
 588       if ( Allocate_Status /= 0 ) then
 589          call da_error(__FILE__,__LINE__, &
 590                       (/"Error in deallocatting GeometryInfo"/))
 591       END IF
 592
 593    end do sensor_loop     ! end loop for sensor
 594
 595    if (trace_use) call da_trace_exit("da_transform_xtoy_crtm")
 596 #else
 597     call da_error(__FILE__,__LINE__, &
 598        (/"Must compile with $CRTM option for radiances"/))
 599 #endif
 600
 601 end subroutine da_transform_xtoy_crtm
 602