ufo_radiance_tlad_mod.F90

Go to the documentation of this file.

! (C) Copyright 2017-2018 UCAR
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.

!> Fortran module to handle tl/ad for radiance observations

module ufo_radiance_tlad_mod

 use iso_c_binding
 use config_mod
 use kinds

 use ufo_geovals_mod, only: ufo_geovals, ufo_geoval, ufo_geovals_get_var
 use ufo_basis_tlad_mod, only: ufo_basis_tlad
 use ufo_vars_mod
 use ufo_radiance_utils_mod
 use crtm_module
 use obsspace_mod

 implicit none
 private

 !> Fortran derived type for radiance trajectory
 type, extends(ufo_basis_tlad), public :: ufo_radiance_tlad
 private
  type(rad_conf) :: rc
  integer :: n_profiles
  integer :: n_layers
  integer :: n_channels
  type(crtm_atmosphere_type), allocatable :: atm_k(:,:)
  type(crtm_surface_type), allocatable :: sfc_k(:,:)
 contains
  procedure :: setup  => ufo_radiance_tlad_setup
  procedure :: delete  => ufo_radiance_tlad_delete
  procedure :: settraj => ufo_radiance_tlad_settraj
  procedure :: simobs_tl  => ufo_radiance_simobs_tl
  procedure :: simobs_ad  => ufo_radiance_simobs_ad
 end type ufo_radiance_tlad

contains

! ------------------------------------------------------------------------------

subroutine ufo_radiance_tlad_setup(self, c_conf)

implicit none
class(ufo_radiance_tlad), intent(inout) :: self
type(c_ptr),              intent(in)    :: c_conf

 call rad_conf_setup(self%rc,c_conf)

end subroutine ufo_radiance_tlad_setup

! ------------------------------------------------------------------------------

subroutine ufo_radiance_tlad_delete(self)

implicit none
class(ufo_radiance_tlad), intent(inout) :: self

 self%ltraj = .false.
 call rad_conf_delete(self%rc)

 if (allocated(self%atm_k)) then
   call crtm_atmosphere_destroy(self%atm_K)
   deallocate(self%atm_k)
 endif

 if (allocated(self%sfc_k)) then
   call crtm_surface_destroy(self%sfc_K)
   deallocate(self%sfc_k)
 endif

end subroutine ufo_radiance_tlad_delete

! ------------------------------------------------------------------------------

subroutine ufo_radiance_tlad_settraj(self, geovals, obss)

implicit none

class(ufo_radiance_tlad), intent(inout) :: self
type(ufo_geovals),        intent(in)    :: geovals
type(c_ptr), value,       intent(in)    :: obss

! Local Variables
character(*), parameter :: PROGRAM_NAME = 'ufo_radiance_mod.F90'
character(255) :: message, version
integer        :: err_stat, alloc_stat
integer        :: n, k1, ierr
type(ufo_geoval), pointer :: temp

! Define the "non-demoninational" arguments
type(CRTM_ChannelInfo_type)             :: chinfo(self%rc%n_Sensors)
type(CRTM_Geometry_type),   allocatable :: geo(:)

! Define the FORWARD variables
type(CRTM_Atmosphere_type), allocatable :: atm(:)
type(CRTM_Surface_type),    allocatable :: sfc(:)
type(CRTM_RTSolution_type), allocatable :: rts(:,:)

! Define the K-MATRIX variables
type(CRTM_RTSolution_type), allocatable :: rts_K(:,:)


 ! Get number of profile and layers from geovals
 ! ---------------------------------------------
 self%n_Profiles = geovals%nobs
 call ufo_geovals_get_var(geovals, var_tv, temp, status=ierr)
 self%n_Layers = temp%nval
 nullify(temp)

 ! Program header
 ! --------------
 call crtm_version( version )
 call program_message( program_name, &
                       'Check/example program for the CRTM Forward and K-Matrix (setTraj) functions using '//&
                       trim(self%rc%ENDIAN_type)//' coefficient datafiles', &
                       'CRTM Version: '//trim(version) )


 ! Initialise all the sensors at once
 ! ----------------------------------
 !** NOTE: CRTM_Init points to the various binary files needed for CRTM.  See the
 !**       CRTM_Lifecycle.f90 for more details.

 write( *,'(/5x,"Initializing the CRTM (setTraj) ...")' )
 err_stat = crtm_init( self%rc%SENSOR_ID, &
            chinfo, &
            file_path=trim(self%rc%COEFFICIENT_PATH), &
            quiet=.true.)
 if ( err_stat /= success ) THEN
   message = 'Error initializing CRTM (setTraj)'
   call display_message( program_name, message, failure )
   stop
 end if

 ! Loop over all sensors. Not necessary if we're calling CRTM for each sensor
 ! ----------------------------------------------------------------------------
 sensor_loop:do n = 1, self%rc%n_Sensors


   ! Determine the number of channels for the current sensor
   ! -------------------------------------------------------
   self%N_Channels = crtm_channelinfo_n_channels(chinfo(n))


   ! Allocate the ARRAYS
   ! -------------------
   allocate( geo( self%n_Profiles )                         , &
             atm( self%n_Profiles )                         , &
             sfc( self%n_Profiles )                         , &
             rts( self%N_Channels, self%n_Profiles )        , &
             self%atm_K( self%N_Channels, self%n_Profiles ) , &
             self%sfc_K( self%N_Channels, self%n_Profiles ) , &
             rts_k( self%N_Channels, self%n_Profiles )      , &
             stat = alloc_stat                                )
   if ( alloc_stat /= 0 ) THEN
      message = 'Error allocating structure arrays (setTraj)'
      call display_message( program_name, message, failure )
      stop
   end if


   ! Create the input FORWARD structure (atm)
   ! ----------------------------------------
   call crtm_atmosphere_create( atm, self%n_Layers, self%rc%n_Absorbers, self%rc%n_Clouds, self%rc%n_Aerosols )
   if ( any(.NOT. crtm_atmosphere_associated(atm)) ) THEN
      message = 'Error allocating CRTM Forward Atmosphere structure (setTraj)'
      CALL display_message( program_name, message, failure )
      stop
   END IF


   ! Create the input FORWARD structure (sfc)
   ! ----------------------------------------
   call crtm_surface_create(sfc, self%N_Channels)
   IF ( any(.NOT. crtm_surface_associated(sfc)) ) THEN
      message = 'Error allocating CRTM Surface structure (setTraj)'
      CALL display_message( program_name, message, failure )
      stop
   END IF


   ! Create output K-MATRIX structure (atm)
   ! --------------------------------------
   call crtm_atmosphere_create( self%atm_K, self%n_Layers, self%rc%n_Absorbers, self%rc%n_Clouds, self%rc%n_Aerosols )
   if ( any(.NOT. crtm_atmosphere_associated(self%atm_K)) ) THEN
      message = 'Error allocating CRTM K-matrix Atmosphere structure (setTraj)'
      CALL display_message( program_name, message, failure )
      stop
   END IF


   ! Create output K-MATRIX structure (sfc)
   ! --------------------------------------
   call crtm_surface_create(self%sfc_K, self%N_Channels)
   IF ( any(.NOT. crtm_surface_associated(self%sfc_K)) ) THEN
      message = 'Error allocating CRTM K-matrix Surface structure (setTraj)'
      CALL display_message( program_name, message, failure )
      stop
   END IF


   !Assign the data from the GeoVaLs
   !--------------------------------
   call load_atm_data(self%N_PROFILES,self%N_LAYERS,geovals,atm)
   call load_sfc_data(self%N_PROFILES,self%N_LAYERS,self%N_Channels,geovals,sfc,chinfo,obss)
   call load_geom_data(obss,geo)


   ! Zero the K-matrix OUTPUT structures
   ! -----------------------------------
   call crtm_atmosphere_zero( self%atm_K )
   call crtm_surface_zero( self%sfc_K )


   ! Inintialize the K-matrix INPUT so that the results are dTb/dx
   ! -------------------------------------------------------------
   rts_k%Radiance               = zero
   rts_k%Brightness_Temperature = one


   ! Call the K-matrix model
   ! -----------------------
   err_stat = crtm_k_matrix( atm         , &  ! FORWARD  Input
                             sfc         , &  ! FORWARD  Input
                             rts_k       , &  ! K-MATRIX Input
                             geo         , &  ! Input
                             chinfo(n:n) , &  ! Input
                             self%atm_K  , &  ! K-MATRIX Output
                             self%sfc_K  , &  ! K-MATRIX Output
                             rts           )  ! FORWARD  Output
   if ( err_stat /= success ) THEN
      message = 'Error calling CRTM (setTraj) K-Matrix Model for '//trim(self%rc%SENSOR_ID(n))
      call display_message( program_name, message, failure )
      stop
   end if


   ! Deallocate the structures
   ! -------------------------
   call crtm_geometry_destroy(geo)
   call crtm_atmosphere_destroy(atm)
   call crtm_rtsolution_destroy(rts_k)
   call crtm_rtsolution_destroy(rts)
   call crtm_surface_destroy(sfc)


   ! Deallocate all arrays
   ! ---------------------
   deallocate(geo, atm, sfc, rts, rts_k, stat = alloc_stat)
   if ( alloc_stat /= 0 ) THEN
      message = 'Error deallocating structure arrays (setTraj)'
      call display_message( program_name, message, failure )
      stop
   end if

 end do sensor_loop


 ! Destroy CRTM instance
 ! ---------------------
 write( *, '( /5x, "Destroying the CRTM (setTraj)..." )' )
 err_stat = crtm_destroy( chinfo )
 if ( err_stat /= success ) THEN
    message = 'Error destroying CRTM (setTraj)'
    call display_message( program_name, message, failure )
    stop
 end if


 ! Set flag that the tracectory was set
 ! ------------------------------------
 self%ltraj = .true.

end subroutine ufo_radiance_tlad_settraj

! ------------------------------------------------------------------------------

subroutine ufo_radiance_simobs_tl(self, geovals, hofx, obss)

implicit none
class(ufo_radiance_tlad), intent(in) :: self
type(ufo_geovals),        intent(in) :: geovals
real(c_double),        intent(inout) :: hofx(:)
type(c_ptr), value,    intent(in)    :: obss

character(len=*), parameter :: myname_="ufo_radiance_simobs_tl"
character(max_string) :: err_msg
integer :: job, jprofile, jchannel, jlevel, ierr
type(ufo_geoval), pointer :: tv_d


 ! Initial checks
 ! --------------

 ! Check if trajectory was set
 if (.not. self%ltraj) then
   write(err_msg,*) myname_, ' trajectory wasnt set!'
   call abor1_ftn(err_msg)
 endif

 ! Check if nobs is consistent in geovals & hofx
 if (geovals%nobs /= self%n_Profiles) then
   write(err_msg,*) myname_, ' error: nobs inconsistent!'
   call abor1_ftn(err_msg)
 endif

 ! Initialize hofx
 ! ---------------
 hofx(:) = 0.0_kind_real


 ! Temperature
 ! -----------

 ! Check if variable is in geovals and get it
 call ufo_geovals_get_var(geovals, var_tv, tv_d, status=ierr )
 if (ierr/=0) then
   write(err_msg,*) myname_, trim(var_tv), ' doesnt exist'
   call abor1_ftn(err_msg)
 endif

 ! Check model levels is consistent in geovals & crtm
 if (tv_d%nval /= self%n_Layers) then
   write(err_msg,*) myname_, ' error: layers inconsistent!'
   call abor1_ftn(err_msg)
 endif

 ! Multiply by Jacobian and add to hofx
 job = 0
 do jprofile = 1, self%n_Profiles
   do jchannel = 1, self%n_Channels
     job = job + 1
     do jlevel = 1, tv_d%nval
       hofx(job) = hofx(job) + &
                    self%atm_K(jchannel,jprofile)%Temperature(jlevel) * tv_d%vals(jlevel,jprofile)
     enddo
   enddo
 enddo


end subroutine ufo_radiance_simobs_tl

! ------------------------------------------------------------------------------

subroutine ufo_radiance_simobs_ad(self, geovals, hofx, obss)

implicit none
class(ufo_radiance_tlad), intent(in) :: self
type(ufo_geovals),     intent(inout) :: geovals
real(c_double),           intent(in) :: hofx(:)
type(c_ptr), value,    intent(in)    :: obss

character(len=*), parameter :: myname_="ufo_radiance_simobs_ad"
character(max_string) :: err_msg
integer :: job, jprofile, jchannel, jlevel, ierr
type(ufo_geoval), pointer :: tv_d


 ! Initial checks
 ! --------------

 ! Check if trajectory was set
 if (.not. self%ltraj) then
   write(err_msg,*) myname_, ' trajectory wasnt set!'
   call abor1_ftn(err_msg)
 endif

 ! Check if nobs is consistent in geovals & hofx
 if (geovals%nobs /= self%n_Profiles) then
   write(err_msg,*) myname_, ' error: nobs inconsistent!'
   call abor1_ftn(err_msg)
 endif


 ! Temperature
 ! -----------

 ! Check if variable is in geovals and get it
 call ufo_geovals_get_var(geovals, var_tv, tv_d, status=ierr )
 if (ierr/=0) then
   write(err_msg,*) myname_, trim(var_tv), ' doesnt exist'
   call abor1_ftn(err_msg)
 endif

 ! allocate if not yet allocated
 if (.not. allocated(tv_d%vals)) then
    tv_d%nobs = self%n_Profiles
    tv_d%nval = self%n_Layers
    allocate(tv_d%vals(tv_d%nval,tv_d%nobs))
    tv_d%vals = 0.0_kind_real
 endif


 ! Multiply by Jacobian and add to hofx (adjoint)
 job = 0
 do jprofile = 1, self%n_Profiles
   do jchannel = 1, self%n_Channels
     job = job + 1
     do jlevel = 1, tv_d%nval
       tv_d%vals(jlevel,jprofile) = tv_d%vals(jlevel,jprofile) + &
                                    self%atm_K(jchannel,jprofile)%Temperature(jlevel) * hofx(job)
     enddo
   enddo
 enddo


 ! Once all geovals set replace flag
 ! ---------------------------------
 if (.not. geovals%linit ) geovals%linit=.true.


end subroutine ufo_radiance_simobs_ad

! ------------------------------------------------------------------------------

end module ufo_radiance_tlad_mod