var/da/da_radiance/da_rttov_direct.inc

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