var/da/da_radiance/da_transform_xtoy_crtm_adj.inc

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