var/da/da_radiance/da_rttov_ad.inc

   1 #ifdef RTTOV
   2 subroutine da_rttov_ad( inst, nchanl, nprofiles, con_vars, &
   3                       aux_vars, con_vars_ad, tb )
   4
   5    !---------------------------------------------------------------------------
   6    ! PURPOSE: interface to the adjoint subroutine of RTTOV
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    integer             ,  intent (in) :: inst, nchanl, nprofiles
  12    type (con_vars_type),  intent (in) :: con_vars (nprofiles)
  13    type (con_vars_type),  intent (inout) :: con_vars_ad (nprofiles)
  14    type (aux_vars_type),  intent (in) :: aux_vars (nprofiles)
  15    real                ,  intent (in) :: tb(nchanl,nprofiles)
  16
  17    ! local variables
  18    integer             :: n, j, asw, nlevels, nchanprof
  19    Integer             :: alloc_status(10)
  20
  21    ! RTTOV input parameters
  22    type (rttov_chanprof), allocatable :: chanprof(:)
  23    type (rttov_profile),   allocatable :: profiles(:), profiles_ad(:)
  24    type (rttov_emissivity), allocatable :: emissivity(:), emissivity_ad(:)
  25    logical,               allocatable :: calcemis(:)
  26
  27    ! RTTOV out parameters
  28    integer                            :: errorstatus
  29
  30    ! RTTOV inout parameters
  31    type (rttov_radiance)               :: radiance, radiance_ad
  32    type (rttov_transmission)           :: transmission, transmission_ad
  33
  34    call da_trace_entry("da_rttov_ad")
  35
  36    nchanprof = nchanl*nprofiles
  37
  38    ! generate the chanprof array
  39    allocate (chanprof(nchanprof))
  40    do n = 1, nprofiles
  41       chanprof((n-1)*nchanl+1:n*nchanl)%prof = n
  42       chanprof((n-1)*nchanl+1:n*nchanl)%chan = (/ (j, j=1,nchanl) /)
  43    end do
  44
  45    alloc_status (:) = 0.0
  46
  47    nlevels = con_vars(1) % nlevels
  48
  49    ! allocate input profile arrays with the number of levels
  50    asw = 1  ! allocate
  51    allocate ( profiles(nprofiles), stat= alloc_status(1))
  52    call rttov_alloc_prof(        &
  53       & errorstatus,             &
  54       & nprofiles,               &
  55       & profiles,                &
  56       & nlevels,                 &
  57       & opts(inst),              &
  58       & asw,                     &
  59       & coefs = coefs(inst),     &  ! mandatory if either opts%addclouds or opts%addaerosl is true
  60       & init = .true.            )  ! additionally initialize profiles structure
  61    if ( errorstatus /= errorstatus_success .or. alloc_status(1) /= 0 ) then
  62      call da_error(__FILE__,__LINE__, &
  63         (/"memory allocation error for profile arrays"/))
  64    end if
  65
  66    asw = 1  ! allocate
  67    allocate ( profiles_ad(nprofiles), stat= alloc_status(2))
  68    call rttov_alloc_prof(        &
  69       & errorstatus,             &
  70       & nprofiles,               &
  71       & profiles_ad,             &
  72       & nlevels,                 &
  73       & opts(inst),              &
  74       & asw,                     &
  75       & coefs = coefs(inst),     &  ! mandatory if either opts%addclouds or opts%addaerosl is true
  76       & init = .true.            )  ! additionally initialize profiles structure
  77    if ( errorstatus /= errorstatus_success .or. alloc_status(2) /= 0 ) then
  78      call da_error(__FILE__,__LINE__, &
  79         (/"memory allocation error for profile AD arrays"/))
  80    end if
  81
  82    do n = 1, nprofiles
  83       profiles(n) % p(:)       = coefs(inst)%coef%ref_prfl_p(:)
  84       profiles(n) % t(:)       = con_vars(n)%t(:)
  85       profiles(n) % q(:)       = con_vars(n)%q(:)
  86       if ( opts(inst) % rt_mw % clw_data ) profiles(n) % clw(:)     = 0.0 !con_vars(n)%clw(:)
  87       if ( opts_rt_ir(inst) % ozone_data ) profiles(n) % o3(:)      = 0.0 !con_vars(n)%o3(:)
  88       if ( opts_rt_ir(inst) % co2_data )   profiles(n) % co2(:)     = 0.0 !con_vars(n)%co2(:)
  89       if ( opts_rt_ir(inst) % n2o_data )   profiles(n) % n2o(:)     = 0.0
  90       if ( opts_rt_ir(inst) % co_data )    profiles(n) % co(:)      = 0.0
  91       if ( opts_rt_ir(inst) % co_data )    profiles(n) % ch4(:)     = 0.0
  92       if ( opts_rt_ir(inst) % addaerosl )  profiles(n) % aerosols(:,:)   = 0.0
  93       if ( opts_rt_ir(inst) % addclouds )  then
  94          profiles(n) % cloud(:,:)      = 0.0
  95          profiles(n) % cfrac(:)        = 0.0
  96          profiles(n) % idg             = 1
  97          profiles(n) % ice_scheme      = 1
  98       end if
  99
 100       profiles(n) % skin % surftype   = aux_vars (n) % surftype
 101       if ( profiles(n)% skin % surftype == 1 ) then
 102          if ( opts_rt_ir(inst) % addsolar ) then
 103             ! if refelcted solar radiation is to be included in the SWIR channels, then
 104             ! specification of fresh or salt water needs to be provided
 105             profiles(n) % skin % watertype = 1
 106          end if
 107       end if
 108
 109       if ( coefs(inst)%coef%id_sensor == sensor_id_mw .or. &
 110            coefs(inst)%coef%id_sensor == sensor_id_po ) then
 111          if ( profiles(n) % skin % surftype == 2 ) then
 112             profiles(n) % skin % fastem (1) = 2.2
 113             profiles(n) % skin % fastem (2) = 3.7
 114             profiles(n) % skin % fastem (3) = 122.0
 115             profiles(n) % skin % fastem (4) = 0.0
 116             profiles(n) % skin % fastem (5) = 0.15
 117          else if ( profiles(n) % skin % surftype == 0 ) then
 118             profiles(n) % skin % fastem (1) = 3.0
 119             profiles(n) % skin % fastem (2) = 5.0
 120             profiles(n) % skin % fastem (3) = 15.0
 121             profiles(n) % skin % fastem (4) = 0.1
 122             profiles(n) % skin % fastem (5) = 0.3
 123          end if
 124       end if
 125
 126       profiles(n) % skin % t          = aux_vars (n) % surft
 127       !profiles(n) % skin % fastem (:) = 0.0  ! aux_vars (n) % fastem (:)
 128
 129       profiles(n) % s2m  % t    = aux_vars (n) % t2m
 130       profiles(n) % s2m  % q    = aux_vars (n) % q2m
 131       profiles(n) % s2m  % o    = 0.0 !aux_vars (n) % o3
 132       profiles(n) % s2m  % p    = con_vars (n) % ps
 133       profiles(n) % s2m  % u    = aux_vars (n) % u10
 134       profiles(n) % s2m  % v    = aux_vars (n) % v10
 135
 136       profiles(n) % zenangle  = aux_vars (n) % satzen
 137       profiles(n) % elevation = 0.001* aux_vars(n) % elevation   ! km
 138       profiles(n) % latitude  = aux_vars(n) % rlat
 139
 140       if ( opts_rt_ir(inst) % addsolar ) then
 141          profiles(n) % azangle     = aux_vars (n) % satazi
 142          profiles(n) % sunzenangle = aux_vars (n) % solzen     !50.0
 143          profiles(n) % sunazangle  = aux_vars (n) % solazi     !86.0
 144          profiles(n) % s2m % wfetc = 100000.0  ! m
 145       end if
 146
 147       profiles(n) % Be          = 0.35   ! optional, for zeeman effect for ssmis and amsua
 148       profiles(n) % cosbk       = 0.0    ! optional, for zeeman effect for ssmis and amsua
 149
 150       profiles(n) % ctp         = 500.0  ! hPa, optional, for simple cloud
 151       profiles(n) % cfraction   = 0.0    ! 0-1, optional, for simple cloud
 152
 153    end do
 154
 155    allocate (emissivity(nchanprof), stat=alloc_status(3))
 156    allocate (emissivity_ad(nchanprof), stat=alloc_status(4))
 157    allocate (calcemis(nchanprof), stat=alloc_status(7))
 158    if ( any( alloc_status /= 0 ) ) then
 159       call da_error(__FILE__,__LINE__, &
 160          (/"memory allocation error for emissivity or calcemis arrays"/))
 161    end if
 162
 163    ! allocate transmittance structure
 164    asw = 1 ! allocate
 165    call rttov_alloc_transmission( &
 166       & errorstatus,              &
 167       & transmission,             &
 168       & nlevels,                  &
 169       & nchanprof,                &
 170       & asw,                      &
 171       & init = .true. )
 172    if ( errorstatus /= errorstatus_success ) then
 173       call da_error(__FILE__,__LINE__, &
 174          (/"memory allocation error for transmittance arrays"/))
 175    end if
 176
 177    asw = 1 ! allocate
 178    call rttov_alloc_transmission( &
 179       & errorstatus,              &
 180       & transmission_ad,          &
 181       & nlevels,                  &
 182       & nchanprof,                &
 183       & asw,                      &
 184       & init = .true. )
 185    if ( errorstatus /= errorstatus_success ) then
 186       call da_error(__FILE__,__LINE__, &
 187          (/"memory allocation error for transmittance AD arrays"/))
 188    end if
 189
 190    ! allocate radiance results arrays with number of channels
 191    asw = 1 ! allocate
 192    call rttov_alloc_rad( &
 193       & errorstatus,     &
 194       & nchanprof,       &
 195       & radiance,        &
 196       & nlevels,         &
 197       & asw,             &
 198       & init=.true. )
 199    if ( errorstatus /= errorstatus_success ) then
 200       call da_error(__FILE__,__LINE__, &
 201          (/"memory allocation error for radiance arrays"/))
 202    end if
 203
 204    asw = 1 ! allocate
 205    call rttov_alloc_rad( &
 206       & errorstatus,     &
 207       & nchanprof,       &
 208       & radiance_ad,     &
 209       & nlevels,         &
 210       & asw,             &
 211       & init=.true. )
 212    if ( errorstatus /= errorstatus_success ) then
 213       call da_error(__FILE__,__LINE__, &
 214          (/"memory allocation error for radiance AD arrays"/))
 215    end if
 216
 217    do n = 1, nprofiles
 218        radiance_ad % bt_clear ((n-1)*nchanl+1:n*nchanl) = tb(1:nchanl,n)
 219        radiance_ad % clear ((n-1)*nchanl+1:n*nchanl) = 0.0
 220    end do
 221
 222    if ( coefs(inst)%coef%id_sensor == sensor_id_ir .or. &
 223         coefs(inst)%coef%id_sensor == sensor_id_hi ) then  ! infrared sensor
 224       calcemis(1:nchanprof)   = .true.
 225       emissivity(1:nchanprof)%emis_in = 0.0
 226       emissivity_ad(1:nchanprof)%emis_in = 0.0
 227    else if ( coefs(inst)%coef%id_sensor == sensor_id_mw .or. &
 228              coefs(inst)%coef%id_sensor == sensor_id_po ) then ! microwave sensor
 229       do n = 1, nprofiles
 230          if ( profiles(n) % skin % surftype == 1) then  ! sea
 231             calcemis((n-1)*nchanl+1:n*nchanl) = .true.
 232             emissivity((n-1)*nchanl+1:n*nchanl)%emis_in = 0.0
 233             emissivity_ad((n-1)*nchanl+1:n*nchanl)%emis_in = 0.0
 234          else                                           ! 0:land ; 2:sea-ice
 235             calcemis((n-1)*nchanl+1:n*nchanl) = .false.
 236             emissivity((n-1)*nchanl+1:n*nchanl)%emis_in = 0.9
 237             emissivity_ad((n-1)*nchanl+1:n*nchanl)%emis_in = 0.0
 238          end if
 239       end do
 240    end if
 241
 242    call  rttov_ad(          &
 243       & errorstatus,       & ! out
 244       & chanprof,          & ! in
 245       & opts(inst),        & ! in
 246       & profiles,          & ! in
 247       & profiles_ad,       & ! inout
 248       & coefs(inst),       & ! in
 249       & transmission,      & ! inout
 250       & transmission_ad,   & ! inout
 251       & radiance,          & ! inout
 252       & radiance_ad,       & ! inout
 253       & calcemis,          & ! in,    optional
 254       & emissivity,        & ! inout, optional
 255       & emissivity_ad)       ! inout, optional
 256
 257    if ( print_detail_rad .or. errorstatus /= errorstatus_success ) then
 258        write (message( 1),*) 'rttov_ad error code    = ', errorstatus
 259        write (message( 2),*) 'nchanl                 = ', nchanl
 260        write (message( 3),*) 'nprofiles              = ', nprofiles
 261        write (message( 4),*) 'calcemis               = ', calcemis
 262        write (message( 5),*) 'profiles%s2m           = ', profiles(1)%s2m
 263        write (message( 6),*) 'profiles%skin          = ', profiles(1)%skin
 264        write (message( 7),*) 'profiles%zenangle      = ', profiles(1)%zenangle
 265        write (message( 8),*) 'profiles%azangle       = ', profiles(1)%azangle
 266        write (message( 9),*) 'profiles%p             = ', profiles(1)%p
 267        write (message(10),*) 'profiles%t             = ', profiles(1)%t
 268        write (message(11),*) 'profiles%q             = ', profiles(1)%q
 269        write (message(12),*) 'emissivity_out         = ', emissivity(:)%emis_out
 270        write (message(13),*) 'radiance               = ', radiance%bt_clear
 271        write (message(14),*) 'profiles_ad%s2m        = ', profiles_ad(1)%s2m
 272        write (message(15),*) 'profiles_ad%skin       = ', profiles_ad(1)%skin
 273        write (message(16),*) 'profiles_ad%zenangle   = ', profiles_ad(1)%zenangle
 274        write (message(17),*) 'profiles_ad%azangle    = ', profiles_ad(1)%azangle
 275        write (message(18),*) 'profiles_ad%p          = ', profiles_ad(1)%p
 276        write (message(19),*) 'profiles_ad%t          = ', profiles_ad(1)%t
 277        write (message(20),*) 'profiles_ad%q          = ', profiles_ad(1)%q
 278        write (message(21),*) 'emissivity_out_ad      = ', emissivity_ad(:)%emis_out
 279        write (message(22),*) 'radiance_ad            = ', radiance_ad%bt_clear
 280
 281        if ( errorstatus /= errorstatus_success ) call da_error(__FILE__,__LINE__,message(1:22))
 282    end if
 283
 284    do n = 1, nprofiles
 285       con_vars_ad(n)%t(:)         = profiles_ad(n) % t(:)
 286       con_vars_ad(n)%q(:)         = profiles_ad(n) % q(:)
 287       con_vars_ad(n)%ps           = profiles_ad(n) % s2m % p
 288    end do
 289
 290    deallocate (emissivity)
 291    deallocate (emissivity_ad)
 292    deallocate (calcemis)
 293    deallocate (chanprof)
 294
 295    asw = 0 ! deallocation
 296    ! deallocate radiance arrays
 297    call rttov_alloc_rad (errorstatus,nchanprof,radiance,nlevels,asw)
 298    if ( errorstatus /= errorstatus_success ) then
 299       call da_error(__FILE__,__LINE__, &
 300         (/"radiance deallocation error"/))
 301    end if
 302    call rttov_alloc_rad (errorstatus,nchanprof,radiance_ad,nlevels,asw)
 303    if ( errorstatus /= errorstatus_success ) then
 304       call da_error(__FILE__,__LINE__, &
 305         (/"radiance AD deallocation error"/))
 306    end if
 307
 308    ! deallocate transmission arrays
 309    call rttov_alloc_transmission (errorstatus,transmission,nlevels,nchanprof,asw)
 310    if ( errorstatus /= errorstatus_success ) then
 311       call da_error(__FILE__,__LINE__, &
 312         (/"transmission deallocation error"/))
 313    end if
 314    call rttov_alloc_transmission (errorstatus,transmission_ad,nlevels,nchanprof,asw)
 315    if ( errorstatus /= errorstatus_success ) then
 316       call da_error(__FILE__,__LINE__, &
 317         (/"transmission AD deallocation error"/))
 318    end if
 319
 320    ! deallocate profile arrays
 321    call rttov_alloc_prof (errorstatus,nprofiles,profiles,nlevels,opts(inst),asw)
 322    deallocate(profiles,stat=alloc_status(8))
 323    if ( errorstatus /= errorstatus_success .or. alloc_status(8) /= 0 ) then
 324       call da_error(__FILE__,__LINE__, &
 325         (/"profile deallocation error"/))
 326    end if
 327    call rttov_alloc_prof (errorstatus,nprofiles,profiles_ad,nlevels,opts(inst),asw)
 328    deallocate(profiles_ad,stat=alloc_status(9))
 329    if ( errorstatus /= errorstatus_success .or. alloc_status(9) /= 0 ) then
 330       call da_error(__FILE__,__LINE__, &
 331         (/"profile AD deallocation error"/))
 332    end if
 333
 334    call da_trace_exit("da_rttov_ad")
 335
 336 end subroutine da_rttov_ad
 337 #endif
 338