var/da/da_radiance/da_rttov_k.inc

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