var/da/da_radiance/da_rttov_tl.inc

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