UFO
ufo_radiancecrtm_tlad_mod.F90
Go to the documentation of this file.
1 ! (C) Copyright 2017-2018 UCAR
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 
6 !> Fortran module to handle tl/ad for radiancecrtm observations
7 
8 module ufo_radiancecrtm_tlad_mod
9 
10  use crtm_module
11 
12  use fckit_configuration_module, only: fckit_configuration
13  use iso_c_binding
14  use kinds
15  use missing_values_mod
16 
17  use obsspace_mod
18 
19  use ufo_geovals_mod, only: ufo_geovals, ufo_geoval, ufo_geovals_get_var
20  use ufo_vars_mod
21  use ufo_crtm_utils_mod
22 
23  use ufo_constants_mod, only: deg2rad
24 
25  implicit none
26  private
27 
28  !> Fortran derived type for radiancecrtm trajectory
29  type, public :: ufo_radiancecrtm_tlad
30  private
31  character(len=MAXVARLEN), public, allocatable :: varin(:) ! variables requested from the model
32  integer, allocatable :: channels(:)
33  type(crtm_conf) :: conf
34  type(crtm_conf) :: conf_traj
35  integer :: n_profiles
36  integer :: n_layers
37  integer :: n_channels
38  type(crtm_atmosphere_type), allocatable :: atm_k(:,:)
39  type(crtm_surface_type), allocatable :: sfc_k(:,:)
40  logical :: ltraj
41  logical, allocatable :: skip_profiles(:)
42  contains
43  procedure :: setup => ufo_radiancecrtm_tlad_setup
44  procedure :: delete => ufo_radiancecrtm_tlad_delete
45  procedure :: settraj => ufo_radiancecrtm_tlad_settraj
46  procedure :: simobs_tl => ufo_radiancecrtm_simobs_tl
47  procedure :: simobs_ad => ufo_radiancecrtm_simobs_ad
48  end type ufo_radiancecrtm_tlad
49 
50  character(len=maxvarlen), dimension(1), parameter :: varin_default = &
51  (/var_ts/)
52 
53 contains
54 
55 ! ------------------------------------------------------------------------------
56 
57 subroutine ufo_radiancecrtm_tlad_setup(self, f_confOper, channels)
58 
59 implicit none
60 class(ufo_radiancecrtm_tlad), intent(inout) :: self
61 type(fckit_configuration), intent(in) :: f_confOper
62 integer(c_int), intent(in) :: channels(:) !List of channels to use
63 
64 integer :: nvars_in
65 integer :: ind, jspec
66 type(fckit_configuration) :: f_confOpts,f_confLinOper
67 
68  call f_confoper%get_or_die("obs options",f_confopts)
69  call crtm_conf_setup(self%conf_traj, f_confopts, f_confoper)
70 
71  if ( f_confoper%has("linear obs operator") ) then
72  call f_confoper%get_or_die("linear obs operator",f_conflinoper)
73  call crtm_conf_setup(self%conf, f_confopts, f_conflinoper)
74  else
75  call crtm_conf_setup(self%conf, f_confopts, f_confoper)
76  end if
77 
78  ! request from the model var_ts +
79  ! 1 * n_Absorbers
80  ! 1 * n_Clouds (mass content only)
81  nvars_in = size(varin_default) + self%conf%n_Absorbers + self%conf%n_Clouds + self%conf%n_Surfaces
82  allocate(self%varin(nvars_in))
83  self%varin(1:size(varin_default)) = varin_default
84  ind = size(varin_default) + 1
85 
86  !Use list of Absorbers and Clouds from conf
87  do jspec = 1, self%conf%n_Absorbers
88  self%varin(ind) = self%conf%Absorbers(jspec)
89  ind = ind + 1
90  end do
91  do jspec = 1, self%conf%n_Clouds
92  self%varin(ind) = self%conf%Clouds(jspec,1)
93  ind = ind + 1
94  end do
95  do jspec = 1, self%conf%n_Surfaces
96  self%varin(ind) = self%conf%Surfaces(jspec)
97  ind = ind + 1
98  end do
99 
100  ! save channels
101  allocate(self%channels(size(channels)))
102  self%channels(:) = channels(:)
103 
104 end subroutine ufo_radiancecrtm_tlad_setup
105 
106 ! ------------------------------------------------------------------------------
107 
108 subroutine ufo_radiancecrtm_tlad_delete(self)
109 
110 implicit none
111 class(ufo_radiancecrtm_tlad), intent(inout) :: self
112 
113  self%ltraj = .false.
114 
115  call crtm_conf_delete(self%conf)
116  call crtm_conf_delete(self%conf_traj)
117 
118  if (allocated(self%atm_k)) then
119  call crtm_atmosphere_destroy(self%atm_K)
120  deallocate(self%atm_k)
121  endif
122 
123  if (allocated(self%sfc_k)) then
124  call crtm_surface_destroy(self%sfc_K)
125  deallocate(self%sfc_k)
126  endif
127 
128  if (allocated(self%Skip_Profiles)) deallocate(self%Skip_Profiles)
129 
130 end subroutine ufo_radiancecrtm_tlad_delete
131 
132 ! ------------------------------------------------------------------------------
133 
134 subroutine ufo_radiancecrtm_tlad_settraj(self, geovals, obss, hofxdiags)
135 use fckit_mpi_module, only: fckit_mpi_comm
136 
137 implicit none
138 
139 class(ufo_radiancecrtm_tlad), intent(inout) :: self
140 type(ufo_geovals), intent(in) :: geovals
141 type(c_ptr), value, intent(in) :: obss
142 type(ufo_geovals), intent(inout) :: hofxdiags !non-h(x) diagnostics
143 
144 ! Local Variables
145 character(*), parameter :: PROGRAM_NAME = 'ufo_radiancecrtm_tlad_settraj'
146 character(255) :: message, version
147 character(max_string) :: err_msg
148 integer :: err_stat, alloc_stat
149 integer :: n
150 type(ufo_geoval), pointer :: temp
151 integer :: jvar, ichannel, jchannel, jprofile, jlevel, jspec
152 real(c_double) :: missing
153 type(fckit_mpi_comm) :: f_comm
154 
155 
156 ! Define the "non-demoninational" arguments
157 type(crtm_channelinfo_type) :: chinfo(self%conf_traj%n_Sensors)
158 type(crtm_geometry_type), allocatable :: geo(:)
159 
160 ! Define the FORWARD variables
161 type(crtm_atmosphere_type), allocatable :: atm(:)
162 type(crtm_surface_type), allocatable :: sfc(:)
163 type(crtm_rtsolution_type), allocatable :: rts(:,:)
164 type(crtm_options_type), allocatable :: Options(:)
165 
166 ! Define the K-MATRIX variables
167 type(crtm_rtsolution_type), allocatable :: rts_K(:,:)
168 
169 ! Used to parse hofxdiags
170 character(len=MAXVARLEN) :: varstr
171 character(len=MAXVARLEN), dimension(hofxdiags%nvar) :: &
172  ystr_diags, xstr_diags
173 character(10), parameter :: jacobianstr = "_jacobian_"
174 integer :: str_pos(4), ch_diags(hofxdiags%nvar)
175 
176 real(kind_real) :: total_od, secant_term, wfunc_max
177 real(kind_real), allocatable :: TmpVar(:)
178 real(kind_real), allocatable :: Tao(:)
179 real(kind_real), allocatable :: Wfunc(:)
180 
181  call obsspace_get_comm(obss, f_comm)
182 
183  ! Get number of profile and layers from geovals
184  ! ---------------------------------------------
185  self%n_Profiles = geovals%nlocs
186  call ufo_geovals_get_var(geovals, var_ts, temp)
187  self%n_Layers = temp%nval
188  nullify(temp)
189 
190  ! Program header
191  ! --------------
192  ! call CRTM_Version( Version )
193  ! call Program_Message( PROGRAM_NAME, &
194  ! 'UFO interface for the CRTM Forward and K-Matrix (setTraj) functions using '//&
195  ! trim(self%conf_traj%ENDIAN_type)//' coefficient datafiles', &
196  ! 'CRTM Version: '//TRIM(Version) )
197 
198 
199  ! Initialise all the sensors at once
200  ! ----------------------------------
201  !** NOTE: CRTM_Init points to the various binary files needed for CRTM. See the
202  !** CRTM_Lifecycle.f90 for more details.
203 
204  ! write( *,'(/5x,"Initializing the CRTM (setTraj) ...")' )
205  err_stat = crtm_init( self%conf_traj%SENSOR_ID, chinfo, &
206  file_path=trim(self%conf_traj%COEFFICIENT_PATH), &
207  irwatercoeff_file=trim(self%conf_traj%IRwaterCoeff_File), &
208  irlandcoeff_file=trim(self%conf_traj%IRlandCoeff_File), &
209  irsnowcoeff_file=trim(self%conf_traj%IRsnowCoeff_File), &
210  iricecoeff_file=trim(self%conf_traj%IRiceCoeff_File), &
211  viswatercoeff_file=trim(self%conf_traj%VISwaterCoeff_File), &
212  vislandcoeff_file=trim(self%conf_traj%VISlandCoeff_File), &
213  vissnowcoeff_file=trim(self%conf_traj%VISsnowCoeff_File), &
214  visicecoeff_file=trim(self%conf_traj%VISiceCoeff_File), &
215  mwwatercoeff_file=trim(self%conf_traj%MWwaterCoeff_File), &
216  quiet=.true.)
217  message = 'Error initializing CRTM (setTraj)'
218  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
219 
220  ! Loop over all sensors. Not necessary if we're calling CRTM for each sensor
221  ! ----------------------------------------------------------------------------
222  sensor_loop:do n = 1, self%conf_traj%n_Sensors
223 
224 
225  ! Pass channel list to CRTM
226  ! -------------------------
227  err_stat = crtm_channelinfo_subset(chinfo(n), self%channels, reset=.false.)
228  message = 'Error subsetting channels'
229  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
230 
231  ! Determine the number of channels for the current sensor
232  ! -------------------------------------------------------
233  self%n_Channels = crtm_channelinfo_n_channels(chinfo(n))
234 
235 
236  ! Allocate the ARRAYS
237  ! -------------------
238  allocate( geo( self%n_Profiles ) , &
239  atm( self%n_Profiles ) , &
240  sfc( self%n_Profiles ) , &
241  rts( self%n_Channels, self%n_Profiles ) , &
242  self%atm_K( self%n_Channels, self%n_Profiles ) , &
243  self%sfc_K( self%n_Channels, self%n_Profiles ) , &
244  rts_k( self%n_Channels, self%n_Profiles ) , &
245  options( self%n_Profiles ) , &
246  stat = alloc_stat )
247  message = 'Error allocating structure arrays (setTraj)'
248  call crtm_comm_stat_check(alloc_stat, program_name, message, f_comm)
249 
250  ! Create the input FORWARD structure (atm)
251  ! ----------------------------------------
252  call crtm_atmosphere_create( atm, self%n_Layers, self%conf_traj%n_Absorbers, self%conf_traj%n_Clouds, self%conf_traj%n_Aerosols )
253  if ( any(.NOT. crtm_atmosphere_associated(atm)) ) THEN
254  message = 'Error allocating CRTM Forward Atmosphere structure (setTraj)'
255  CALL display_message( program_name, message, failure )
256  stop
257  END IF
258 
259  if (self%n_Layers > 0) CALL crtm_rtsolution_create(rts, self%n_Layers )
260 
261  ! Create the input FORWARD structure (sfc)
262  ! ----------------------------------------
263  call crtm_surface_create(sfc, self%n_Channels)
264  IF ( any(.NOT. crtm_surface_associated(sfc)) ) THEN
265  message = 'Error allocating CRTM Surface structure (setTraj)'
266  CALL display_message( program_name, message, failure )
267  stop
268  END IF
269 
270 
271  ! Create output K-MATRIX structure (atm)
272  ! --------------------------------------
273  call crtm_atmosphere_create( self%atm_K, self%n_Layers, self%conf_traj%n_Absorbers, &
274  self%conf_traj%n_Clouds, self%conf_traj%n_Aerosols )
275  if ( any(.NOT. crtm_atmosphere_associated(self%atm_K)) ) THEN
276  message = 'Error allocating CRTM K-matrix Atmosphere structure (setTraj)'
277  CALL display_message( program_name, message, failure )
278  stop
279  END IF
280 
281 
282  ! Create output K-MATRIX structure (sfc)
283  ! --------------------------------------
284  call crtm_surface_create(self%sfc_K, self%n_Channels)
285  IF ( any(.NOT. crtm_surface_associated(self%sfc_K)) ) THEN
286  message = 'Error allocating CRTM K-matrix Surface structure (setTraj)'
287  CALL display_message( program_name, message, failure )
288  stop
289  END IF
290 
291  !Assign the data from the GeoVaLs
292  !--------------------------------
293  call load_atm_data(self%N_PROFILES,self%N_LAYERS,geovals,atm,self%conf_traj)
294  call load_sfc_data(self%N_PROFILES,self%n_Channels,self%channels,geovals,sfc,chinfo,obss,self%conf_traj)
295  call load_geom_data(obss,geo)
296 
297  ! Zero the K-matrix OUTPUT structures
298  ! -----------------------------------
299  call crtm_atmosphere_zero( self%atm_K )
300  call crtm_surface_zero( self%sfc_K )
301 
302 
303  ! Inintialize the K-matrix INPUT so that the results are dTb/dx
304  ! -------------------------------------------------------------
305  rts_k%Radiance = zero
306  rts_k%Brightness_Temperature = one
307 
308  if (allocated(self%Skip_Profiles)) deallocate(self%Skip_Profiles)
309  allocate(self%Skip_Profiles(self%n_Profiles))
310  call ufo_crtm_skip_profiles(self%n_Profiles,self%n_Channels,self%channels,obss,self%Skip_Profiles)
311  profile_loop: do jprofile = 1, self%n_Profiles
312  options(jprofile)%Skip_Profile = self%Skip_Profiles(jprofile)
313  ! check for pressure monotonicity
314  do jlevel = atm(jprofile)%n_layers, 1, -1
315  if ( atm(jprofile)%level_pressure(jlevel) <= atm(jprofile)%level_pressure(jlevel-1) ) then
316  options(jprofile)%Skip_Profile = .true.
317  cycle profile_loop
318  end if
319  end do
320  end do profile_loop
321 
322  ! Call the K-matrix model
323  ! -----------------------
324  err_stat = crtm_k_matrix( atm , & ! FORWARD Input
325  sfc , & ! FORWARD Input
326  rts_k , & ! K-MATRIX Input
327  geo , & ! Input
328  chinfo(n:n) , & ! Input
329  self%atm_K , & ! K-MATRIX Output
330  self%sfc_K , & ! K-MATRIX Output
331  rts , & ! FORWARD Output
332  options ) ! Input
333  message = 'Error calling CRTM (setTraj) K-Matrix Model for '//trim(self%conf_traj%SENSOR_ID(n))
334  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
335 
336  !call CRTM_RTSolution_Inspect(rts)
337 
338 
339  !! Parse hofxdiags%variables into independent/dependent variables and channel
340  !! assumed formats:
341  !! jacobian var --> <ystr>_jacobian_<xstr>_<chstr>
342  !! non-jacobian var --> <ystr>_<chstr>
343 
344  ch_diags = -9999
345  do jvar = 1, hofxdiags%nvar
346  varstr = hofxdiags%variables(jvar)
347  str_pos(4) = len_trim(varstr)
348  if (str_pos(4) < 1) cycle
349  str_pos(3) = index(varstr,"_",back=.true.) !final "_" before channel
350  read(varstr(str_pos(3)+1:str_pos(4)),*, err=999) ch_diags(jvar)
351  999 str_pos(1) = index(varstr,jacobianstr) - 1 !position before jacobianstr
352  if (str_pos(1) == 0) then
353  write(err_msg,*) 'ufo_radiancecrtm_simobs: _jacobian_ must be // &
354  & preceded by dependent variable in config: ', &
355  & hofxdiags%variables(jvar)
356  call abor1_ftn(err_msg)
357  else if (str_pos(1) > 0) then
358  !Diagnostic is a Jacobian member (dy/dx)
359  ystr_diags(jvar) = varstr(1:str_pos(1))
360  str_pos(2) = str_pos(1) + len(jacobianstr) + 1 !begin xstr_diags
361  str_pos(4) = str_pos(3) - str_pos(2)
362  xstr_diags(jvar)(1:str_pos(4)) = varstr(str_pos(2):str_pos(3)-1)
363  xstr_diags(jvar)(str_pos(4)+1:) = ""
364  else !null
365  !Diagnostic is a dependent variable (y)
366  xstr_diags(jvar) = ""
367  ystr_diags(jvar)(1:str_pos(3)-1) = varstr(1:str_pos(3)-1)
368  ystr_diags(jvar)(str_pos(3):) = ""
369  if (ch_diags(jvar) < 0) ystr_diags(jvar) = varstr
370  end if
371  end do
372 
373  ! Set missing value
374  missing = missing_value(missing)
375 
376  ! Put simulated diagnostics into hofxdiags
377  ! ----------------------------------------------
378  do jvar = 1, hofxdiags%nvar
379 
380  if (len(trim(hofxdiags%variables(jvar))) < 1) cycle
381 
382  if (ch_diags(jvar) > 0) then
383  if (size(pack(self%channels,self%channels==ch_diags(jvar))) /= 1) then
384  write(err_msg,*) 'ufo_radiancecrtm_simobs: mismatch between// &
385  & h(x) channels(', self%channels,') and// &
386  & ch_diags(jvar) = ', ch_diags(jvar)
387  call abor1_ftn(err_msg)
388  end if
389  end if
390 
391  do ichannel = 1, size(self%channels)
392  if (ch_diags(jvar) == self%channels(ichannel)) then
393  jchannel = ichannel
394  exit
395  end if
396  end do
397 
398  if (allocated(hofxdiags%geovals(jvar)%vals)) &
399  deallocate(hofxdiags%geovals(jvar)%vals)
400 
401  !============================================
402  ! Diagnostics used for QC and bias correction
403  !============================================
404  if (trim(xstr_diags(jvar)) == "") then
405  ! forward h(x) diags
406  select case(ystr_diags(jvar))
407  ! variable: optical_thickness_of_atmosphere_layer_CH
408  case (var_opt_depth)
409  hofxdiags%geovals(jvar)%nval = self%n_Layers
410  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,self%n_Profiles))
411  hofxdiags%geovals(jvar)%vals = missing
412  do jprofile = 1, self%n_Profiles
413  if (.not.self%Skip_Profiles(jprofile)) then
414  do jlevel = 1, hofxdiags%geovals(jvar)%nval
415  hofxdiags%geovals(jvar)%vals(jlevel,jprofile) = &
416  rts(jchannel,jprofile) % layer_optical_depth(jlevel)
417  end do
418  end if
419  end do
420 
421  ! variable: brightness_temperature_CH
422  case (var_tb)
423  hofxdiags%geovals(jvar)%nval = 1
424  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,self%n_Profiles))
425  hofxdiags%geovals(jvar)%vals = missing
426  do jprofile = 1, self%n_Profiles
427  if (.not.self%Skip_Profiles(jprofile)) then
428  hofxdiags%geovals(jvar)%vals(1,jprofile) = &
429  rts(jchannel,jprofile) % Brightness_Temperature
430  end if
431  end do
432 
433  ! variable: transmittances_of_atmosphere_layer_CH
434  case (var_lvl_transmit)
435  hofxdiags%geovals(jvar)%nval = self%n_Layers
436  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,self%n_Profiles))
437  hofxdiags%geovals(jvar)%vals = missing
438  allocate(tmpvar(self%n_Profiles))
439  call obsspace_get_db(obss, "MetaData", "sensor_zenith_angle", tmpvar)
440  do jprofile = 1, self%n_Profiles
441  if (.not.self%Skip_Profiles(jprofile)) then
442  secant_term = one/cos(tmpvar(jprofile)*deg2rad)
443  total_od = 0.0
444  do jlevel = 1, self%n_Layers
445  total_od = total_od + rts(jchannel,jprofile) % layer_optical_depth(jlevel)
446  hofxdiags%geovals(jvar)%vals(jlevel,jprofile) = &
447  exp(-min(limit_exp,total_od*secant_term))
448  end do
449  end if
450  end do
451  deallocate(tmpvar)
452 
453  ! variable: weightingfunction_of_atmosphere_layer_CH
454  case (var_lvl_weightfunc)
455  hofxdiags%geovals(jvar)%nval = self%n_Layers
456  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,self%n_Profiles))
457  hofxdiags%geovals(jvar)%vals = missing
458  allocate(tmpvar(self%n_Profiles))
459  allocate(tao(self%n_Layers))
460  call obsspace_get_db(obss, "MetaData", "sensor_zenith_angle", tmpvar)
461  do jprofile = 1, self%n_Profiles
462  if (.not.self%Skip_Profiles(jprofile)) then
463  ! get layer-to-space transmittance
464  secant_term = one/cos(tmpvar(jprofile)*deg2rad)
465  total_od = 0.0
466  do jlevel = 1, self%n_Layers
467  total_od = total_od + rts(jchannel,jprofile) % layer_optical_depth(jlevel)
468  tao(jlevel) = exp(-min(limit_exp,total_od*secant_term))
469  end do
470  ! get weighting function
471  do jlevel = self%n_Layers-1, 1, -1
472  hofxdiags%geovals(jvar)%vals(jlevel,jprofile) = &
473  abs( (tao(jlevel+1)-tao(jlevel))/ &
474  (log(atm(jprofile)%pressure(jlevel+1))- &
475  log(atm(jprofile)%pressure(jlevel))) )
476  end do
477  hofxdiags%geovals(jvar)%vals(self%n_Layers,jprofile) = &
478  hofxdiags%geovals(jvar)%vals(self%n_Layers-1,jprofile)
479  end if
480  end do
481  deallocate(tmpvar)
482  deallocate(tao)
483 
484  ! variable: pressure_level_at_peak_of_weightingfunction_CH
485  case (var_pmaxlev_weightfunc)
486  hofxdiags%geovals(jvar)%nval = 1
487  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,self%n_Profiles))
488  hofxdiags%geovals(jvar)%vals = missing
489  allocate(tmpvar(self%n_Profiles))
490  allocate(tao(self%n_Layers))
491  allocate(wfunc(self%n_Layers))
492  call obsspace_get_db(obss, "MetaData", "sensor_zenith_angle", tmpvar)
493  do jprofile = 1, self%n_Profiles
494  if (.not.self%Skip_Profiles(jprofile)) then
495  ! get layer-to-space transmittance
496  secant_term = one/cos(tmpvar(jprofile)*deg2rad)
497  total_od = 0.0
498  do jlevel = 1, self%n_Layers
499  total_od = total_od + rts(jchannel,jprofile) % layer_optical_depth(jlevel)
500  tao(jlevel) = exp(-min(limit_exp,total_od*secant_term))
501  end do
502  ! get weighting function
503  do jlevel = self%n_Layers-1, 1, -1
504  wfunc(jlevel) = &
505  abs( (tao(jlevel+1)-tao(jlevel))/ &
506  (log(atm(jprofile)%pressure(jlevel+1))- &
507  log(atm(jprofile)%pressure(jlevel))) )
508  end do
509  wfunc(self%n_Layers) = wfunc(self%n_Layers-1)
510  ! get pressure level at the peak of the weighting function
511  wfunc_max = -999.0
512  do jlevel = self%n_Layers-1, 1, -1
513  if (wfunc(jlevel) > wfunc_max) then
514  wfunc_max = wfunc(jlevel)
515  hofxdiags%geovals(jvar)%vals(1,jprofile) = jlevel
516  endif
517  enddo
518  end if
519  end do
520  deallocate(tmpvar)
521  deallocate(tao)
522  deallocate(wfunc)
523 
524  case default
525  write(err_msg,*) 'ufo_radiancecrtm_simobs: //&
526  & ObsDiagnostic is unsupported, ', &
527  & hofxdiags%variables(jvar)
528  ! call abor1_ftn(err_msg)
529  hofxdiags%geovals(jvar)%nval = 1
530  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,self%n_Profiles))
531  hofxdiags%geovals(jvar)%vals = missing
532  end select
533  else if (ystr_diags(jvar) == var_tb) then
534  ! var_tb jacobians
535  select case (xstr_diags(jvar))
536  ! variable: brightness_temperature_jacobian_surface_emissivity_CH (nval=1)
537  case (var_sfc_emiss)
538  hofxdiags%geovals(jvar)%nval = 1
539  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,self%n_Profiles))
540  hofxdiags%geovals(jvar)%vals = missing
541  do jprofile = 1, self%n_Profiles
542  if (.not.self%Skip_Profiles(jprofile)) then
543  hofxdiags%geovals(jvar)%vals(1,jprofile) = &
544  rts_k(jchannel,jprofile) % surface_emissivity
545  end if
546  end do
547  case default
548  write(err_msg,*) 'ufo_radiancecrtm_simobs: //&
549  & ObsDiagnostic is unsupported, ', &
550  & hofxdiags%variables(jvar)
551  call abor1_ftn(err_msg)
552  end select
553  else
554  write(err_msg,*) 'ufo_radiancecrtm_simobs: //&
555  & ObsDiagnostic is unsupported, ', &
556  & hofxdiags%variables(jvar)
557  call abor1_ftn(err_msg)
558  end if
559  end do
560 
561  ! Deallocate the structures
562  ! -------------------------
563  call crtm_geometry_destroy(geo)
564  call crtm_atmosphere_destroy(atm)
565  call crtm_rtsolution_destroy(rts_k)
566  call crtm_rtsolution_destroy(rts)
567  call crtm_surface_destroy(sfc)
568 
569 
570  ! Deallocate all arrays
571  ! ---------------------
572  deallocate(geo, atm, sfc, rts, rts_k, options, stat = alloc_stat)
573  message = 'Error deallocating structure arrays (setTraj)'
574  call crtm_comm_stat_check(alloc_stat, program_name, message, f_comm)
575 
576  end do sensor_loop
577 
578 
579  ! Destroy CRTM instance
580  ! ---------------------
581  ! write( *, '( /5x, "Destroying the CRTM (setTraj)..." )' )
582  err_stat = crtm_destroy( chinfo )
583  message = 'Error destroying CRTM (setTraj)'
584  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
585 
586  ! Set flag that the tracectory was set
587  ! ------------------------------------
588  self%ltraj = .true.
589 
590 end subroutine ufo_radiancecrtm_tlad_settraj
591 
592 ! ------------------------------------------------------------------------------
593 
594 subroutine ufo_radiancecrtm_simobs_tl(self, geovals, obss, nvars, nlocs, hofx)
595 
596 implicit none
597 class(ufo_radiancecrtm_tlad), intent(in) :: self
598 type(ufo_geovals), intent(in) :: geovals
599 type(c_ptr), value, intent(in) :: obss
600 integer, intent(in) :: nvars, nlocs
601 real(c_double), intent(inout) :: hofx(nvars, nlocs)
602 
603 character(len=*), parameter :: myname_="ufo_radiancecrtm_simobs_tl"
604 character(max_string) :: err_msg
605 integer :: jprofile, jchannel, jlevel, jspec, ispec
606 type(ufo_geoval), pointer :: geoval_d
607 
608  ! Initial checks
609  ! --------------
610 
611  ! Check if trajectory was set
612  if (.not. self%ltraj) then
613  write(err_msg,*) myname_, ' trajectory wasnt set!'
614  call abor1_ftn(err_msg)
615  endif
616 
617  ! Check if nlocs is consistent in geovals & hofx
618  if (geovals%nlocs /= self%n_Profiles) then
619  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
620  call abor1_ftn(err_msg)
621  endif
622 
623  ! Initialize hofx
624  ! ---------------
625  hofx(:,:) = 0.0_kind_real
626 
627  ! Temperature
628  ! -----------
629 
630  ! Get t from geovals
631  call ufo_geovals_get_var(geovals, var_ts, geoval_d)
632 
633  ! Check model levels is consistent in geovals & crtm
634  if (geoval_d%nval /= self%n_Layers) then
635  write(err_msg,*) myname_, ' error: layers inconsistent!'
636  call abor1_ftn(err_msg)
637  endif
638 
639  ! Multiply by Jacobian and add to hofx
640  do jprofile = 1, self%n_Profiles
641  if (.not.self%Skip_Profiles(jprofile)) then
642  do jchannel = 1, size(self%channels)
643  do jlevel = 1, geoval_d%nval
644  hofx(jchannel, jprofile) = hofx(jchannel, jprofile) + &
645  self%atm_K(jchannel,jprofile)%Temperature(jlevel) * &
646  geoval_d%vals(jlevel,jprofile)
647  enddo
648  enddo
649  end if
650  enddo
651 
652  ! Absorbers
653  ! ---------
654 
655  do jspec = 1, self%conf%n_Absorbers
656  ! Get Absorber from geovals
657  call ufo_geovals_get_var(geovals, self%conf%Absorbers(jspec), geoval_d)
658 
659  ispec = ufo_vars_getindex(self%conf_traj%Absorbers, self%conf%Absorbers(jspec))
660 
661  ! Multiply by Jacobian and add to hofx
662  do jprofile = 1, self%n_Profiles
663  if (.not.self%Skip_Profiles(jprofile)) then
664  do jchannel = 1, size(self%channels)
665  do jlevel = 1, geoval_d%nval
666  hofx(jchannel, jprofile) = hofx(jchannel, jprofile) + &
667  self%atm_K(jchannel,jprofile)%Absorber(jlevel,ispec) * &
668  geoval_d%vals(jlevel,jprofile)
669  enddo
670  enddo
671  end if
672  enddo
673  end do
674 
675  ! Clouds (mass content only)
676  ! --------------------------
677 
678  do jspec = 1, self%conf%n_Clouds
679  ! Get Cloud from geovals
680  call ufo_geovals_get_var(geovals, self%conf%Clouds(jspec,1), geoval_d)
681 
682  ispec = ufo_vars_getindex(self%conf_traj%Clouds(:,1), self%conf%Clouds(jspec,1))
683 
684  ! Multiply by Jacobian and add to hofx
685  do jprofile = 1, self%n_Profiles
686  if (.not.self%Skip_Profiles(jprofile)) then
687  do jchannel = 1, size(self%channels)
688  do jlevel = 1, geoval_d%nval
689  hofx(jchannel, jprofile) = hofx(jchannel, jprofile) + &
690  self%atm_K(jchannel,jprofile)%Cloud(ispec)%Water_Content(jlevel) * &
691  geoval_d%vals(jlevel,jprofile)
692  enddo
693  enddo
694  end if
695  enddo
696  end do
697 
698  ! Surface Variables
699  ! --------------------------
700 
701  do jspec = 1, self%conf%n_Surfaces
702  ! Get Surface from geovals
703  call ufo_geovals_get_var(geovals, self%conf%Surfaces(jspec), geoval_d)
704 
705  select case(self%conf%Surfaces(jspec))
706 
707  case(var_sfc_wtmp)
708 
709  ! Multiply by Jacobian and add to hofx
710  do jprofile = 1, self%n_Profiles
711  if (.not.self%Skip_Profiles(jprofile)) then
712  do jchannel = 1, size(self%channels)
713  jlevel = 1
714  hofx(jchannel, jprofile) = hofx(jchannel, jprofile) + &
715  self%sfc_K(jchannel,jprofile)%water_temperature * &
716  geoval_d%vals(jlevel,jprofile)
717  enddo
718  end if
719  enddo
720  case(var_sfc_wspeed)
721  ! Multiply by Jacobian and add to hofx
722  do jprofile = 1, self%n_Profiles
723  if (.not.self%Skip_Profiles(jprofile)) then
724  do jchannel = 1, size(self%channels)
725  jlevel = 1
726  hofx(jchannel, jprofile) = hofx(jchannel, jprofile) + &
727  self%sfc_K(jchannel,jprofile)%wind_speed * &
728  geoval_d%vals(jlevel,jprofile)
729  enddo
730  end if
731  enddo
732  case(var_sfc_wdir)
733  ! Multiply by Jacobian and add to hofx
734  do jprofile = 1, self%n_Profiles
735  if (.not.self%Skip_Profiles(jprofile)) then
736  do jchannel = 1, size(self%channels)
737  jlevel = 1
738  hofx(jchannel, jprofile) = hofx(jchannel, jprofile) + &
739  self%sfc_K(jchannel,jprofile)%wind_direction * &
740  geoval_d%vals(jlevel,jprofile)
741  enddo
742  end if
743  enddo
744  case(var_sfc_sss)
745  ! Multiply by Jacobian and add to hofx
746  do jprofile = 1, self%n_Profiles
747  if (.not.self%Skip_Profiles(jprofile)) then
748  do jchannel = 1, size(self%channels)
749  jlevel = 1
750  hofx(jchannel, jprofile) = hofx(jchannel, jprofile) + &
751  self%sfc_K(jchannel,jprofile)%salinity * &
752  geoval_d%vals(jlevel,jprofile)
753  enddo
754  end if
755  enddo
756 
757  end select
758  end do
759 
760 
761 end subroutine ufo_radiancecrtm_simobs_tl
762 
763 
764 ! ------------------------------------------------------------------------------
765 
766 subroutine ufo_radiancecrtm_simobs_ad(self, geovals, obss, nvars, nlocs, hofx)
767 
768 implicit none
769 class(ufo_radiancecrtm_tlad), intent(in) :: self
770 type(ufo_geovals), intent(inout) :: geovals
771 type(c_ptr), value, intent(in) :: obss
772 integer, intent(in) :: nvars, nlocs
773 real(c_double), intent(in) :: hofx(nvars, nlocs)
774 
775 character(len=*), parameter :: myname_="ufo_radiancecrtm_simobs_ad"
776 character(max_string) :: err_msg
777 integer :: jprofile, jchannel, jlevel, jspec, ispec
778 type(ufo_geoval), pointer :: geoval_d
779 real(c_double) :: missing
780 
781 
782  ! Initial checks
783  ! --------------
784 
785  ! Check if trajectory was set
786  if (.not. self%ltraj) then
787  write(err_msg,*) myname_, ' trajectory wasnt set!'
788  call abor1_ftn(err_msg)
789  endif
790 
791  ! Check if nlocs is consistent in geovals & hofx
792  if (geovals%nlocs /= self%n_Profiles) then
793  write(err_msg,*) myname_, ' error: nlocs inconsistent!'
794  call abor1_ftn(err_msg)
795  endif
796 
797  ! Set missing value
798  missing = missing_value(missing)
799 
800  ! Temperature
801  ! -----------
802 
803  ! Get t from geovals
804  call ufo_geovals_get_var(geovals, var_ts, geoval_d)
805 
806  ! allocate if not yet allocated
807  if (.not. allocated(geoval_d%vals)) then
808  geoval_d%nlocs = self%n_Profiles
809  geoval_d%nval = self%n_Layers
810  allocate(geoval_d%vals(geoval_d%nval,geoval_d%nlocs))
811  geoval_d%vals = 0.0_kind_real
812  endif
813 
814  ! Multiply by Jacobian and add to hofx (adjoint)
815  do jprofile = 1, self%n_Profiles
816  if (.not.self%Skip_Profiles(jprofile)) then
817  do jchannel = 1, size(self%channels)
818  if (hofx(jchannel, jprofile) /= missing) then
819  do jlevel = 1, geoval_d%nval
820  geoval_d%vals(jlevel,jprofile) = geoval_d%vals(jlevel,jprofile) + &
821  self%atm_K(jchannel,jprofile)%Temperature(jlevel) * &
822  hofx(jchannel, jprofile)
823  enddo
824  endif
825  enddo
826  end if
827  enddo
828 
829  ! Absorbers
830  ! ---------
831 
832  do jspec = 1, self%conf%n_Absorbers
833  ! Get Absorber from geovals
834  call ufo_geovals_get_var(geovals, self%conf%Absorbers(jspec), geoval_d)
835 
836  ! allocate if not yet allocated
837  if (.not. allocated(geoval_d%vals)) then
838  geoval_d%nlocs = self%n_Profiles
839  geoval_d%nval = self%n_Layers
840  allocate(geoval_d%vals(geoval_d%nval,geoval_d%nlocs))
841  geoval_d%vals = 0.0_kind_real
842  endif
843 
844  ispec = ufo_vars_getindex(self%conf_traj%Absorbers, self%conf%Absorbers(jspec))
845 
846  ! Multiply by Jacobian and add to hofx (adjoint)
847  do jprofile = 1, self%n_Profiles
848  if (.not.self%Skip_Profiles(jprofile)) then
849  do jchannel = 1, size(self%channels)
850  if (hofx(jchannel, jprofile) /= missing) then
851  do jlevel = 1, geoval_d%nval
852  geoval_d%vals(jlevel,jprofile) = geoval_d%vals(jlevel,jprofile) + &
853  self%atm_K(jchannel,jprofile)%Absorber(jlevel,ispec) * &
854  hofx(jchannel, jprofile)
855  enddo
856  endif
857  enddo
858  end if
859  enddo
860  end do
861 
862  ! Clouds (mass content only)
863  ! --------------------------
864 
865  do jspec = 1, self%conf%n_Clouds
866  ! Get Cloud from geovals
867  call ufo_geovals_get_var(geovals, self%conf%Clouds(jspec,1), geoval_d)
868 
869  ! allocate if not yet allocated
870  if (.not. allocated(geoval_d%vals)) then
871  geoval_d%nlocs = self%n_Profiles
872  geoval_d%nval = self%n_Layers
873  allocate(geoval_d%vals(geoval_d%nval,geoval_d%nlocs))
874  geoval_d%vals = 0.0_kind_real
875  endif
876 
877  ispec = ufo_vars_getindex(self%conf_traj%Clouds(:,1), self%conf%Clouds(jspec,1))
878 
879  ! Multiply by Jacobian and add to hofx (adjoint)
880  do jprofile = 1, self%n_Profiles
881  if (.not.self%Skip_Profiles(jprofile)) then
882  do jchannel = 1, size(self%channels)
883  if (hofx(jchannel, jprofile) /= missing) then
884  do jlevel = 1, geoval_d%nval
885  geoval_d%vals(jlevel,jprofile) = geoval_d%vals(jlevel,jprofile) + &
886  self%atm_K(jchannel,jprofile)%Cloud(ispec)%Water_Content(jlevel) * &
887  hofx(jchannel, jprofile)
888  enddo
889  endif
890  enddo
891  end if
892  enddo
893  end do
894 
895 
896  ! Surface Variables
897  ! --------------------------
898  do jspec = 1, self%conf%n_Surfaces
899  ! Get Cloud from geovals
900  call ufo_geovals_get_var(geovals, self%conf%Surfaces(jspec), geoval_d)
901 
902  ! allocate if not yet allocated
903  if (.not. allocated(geoval_d%vals)) then
904  geoval_d%nlocs = self%n_Profiles
905  geoval_d%nval = self%n_Layers
906  allocate(geoval_d%vals(geoval_d%nval,geoval_d%nlocs))
907  geoval_d%vals = 0.0_kind_real
908  endif
909 
910  select case(self%conf%Surfaces(jspec))
911 
912  case(var_sfc_wtmp)
913  ! Multiply by Jacobian and add to hofx
914  do jprofile = 1, self%n_Profiles
915  if (.not.self%Skip_Profiles(jprofile)) then
916  do jchannel = 1, size(self%channels)
917  if (hofx(jchannel, jprofile) /= missing) then
918  jlevel = 1
919  geoval_d%vals(jlevel, jprofile) = geoval_d%vals(jlevel,jprofile) + &
920  self%sfc_K(jchannel,jprofile)%water_temperature * &
921  hofx(jchannel,jprofile)
922  endif
923  enddo
924  end if
925  enddo
926  case(var_sfc_wspeed)
927  ! Multiply by Jacobian and add to hofx
928  do jprofile = 1, self%n_Profiles
929  if (.not.self%Skip_Profiles(jprofile)) then
930  do jchannel = 1, size(self%channels)
931  if (hofx(jchannel, jprofile) /= missing) then
932  jlevel = 1
933  geoval_d%vals(jlevel, jprofile) = geoval_d%vals(jlevel,jprofile) + &
934  self%sfc_K(jchannel,jprofile)%wind_speed * &
935  hofx(jchannel,jprofile)
936  endif
937  enddo
938  end if
939  enddo
940  case(var_sfc_wdir)
941  ! Multiply by Jacobian and add to hofx
942  do jprofile = 1, self%n_Profiles
943  if (.not.self%Skip_Profiles(jprofile)) then
944  do jchannel = 1, size(self%channels)
945  if (hofx(jchannel, jprofile) /= missing) then
946  jlevel = 1
947  geoval_d%vals(jlevel, jprofile) = geoval_d%vals(jlevel,jprofile) + &
948  self%sfc_K(jchannel,jprofile)%wind_direction * &
949  hofx(jchannel,jprofile)
950  endif
951  enddo
952  end if
953  enddo
954  case(var_sfc_sss)
955  ! Multiply by Jacobian and add to hofx
956  do jprofile = 1, self%n_Profiles
957  if (.not.self%Skip_Profiles(jprofile)) then
958  do jchannel = 1, size(self%channels)
959  if (hofx(jchannel, jprofile) /= missing) then
960  jlevel = 1
961  geoval_d%vals(jlevel, jprofile) = geoval_d%vals(jlevel,jprofile) + &
962  self%sfc_K(jchannel,jprofile)%salinity * &
963  hofx(jchannel,jprofile)
964  endif
965  enddo
966  end if
967  enddo
968 
969  end select
970 
971  enddo
972 
973  ! Once all geovals set replace flag
974  ! ---------------------------------
975  if (.not. geovals%linit ) geovals%linit=.true.
976 
977 
978 end subroutine ufo_radiancecrtm_simobs_ad
979 
980 ! ------------------------------------------------------------------------------
981 
982 end module ufo_radiancecrtm_tlad_mod
ufo_crtm_utils_mod::load_atm_data
subroutine, public load_atm_data(n_Profiles, n_Layers, geovals, atm, conf)
Definition: ufo_crtm_utils_mod.F90:497
ufo_radiancecrtm_tlad_mod::ufo_radiancecrtm_tlad_settraj
subroutine ufo_radiancecrtm_tlad_settraj(self, geovals, obss, hofxdiags)
Definition: ufo_radiancecrtm_tlad_mod.F90:135
ufo_vars_mod::var_sfc_wdir
character(len=maxvarlen), parameter, public var_sfc_wdir
Definition: ufo_variables_mod.F90:62
ufo_crtm_utils_mod::crtm_comm_stat_check
subroutine, public crtm_comm_stat_check(stat, PROGRAM_NAME, message, f_comm)
Definition: ufo_crtm_utils_mod.F90:428
ufo_radiancecrtm_tlad_mod::ufo_radiancecrtm_tlad_delete
subroutine ufo_radiancecrtm_tlad_delete(self)
Definition: ufo_radiancecrtm_tlad_mod.F90:109
ufo_crtm_utils_mod::load_geom_data
subroutine, public load_geom_data(obss, geo)
Definition: ufo_crtm_utils_mod.F90:770
ufo_vars_mod::var_lvl_weightfunc
character(len=maxvarlen), parameter, public var_lvl_weightfunc
Definition: ufo_variables_mod.F90:81
ufo_radiancecrtm_tlad_mod::ufo_radiancecrtm_tlad
Fortran derived type for radiancecrtm trajectory.
Definition: ufo_radiancecrtm_tlad_mod.F90:29
ufo_radiancecrtm_tlad_mod
Fortran module to handle tl/ad for radiancecrtm observations.
Definition: ufo_radiancecrtm_tlad_mod.F90:8
ufo_vars_mod::var_sfc_sss
character(len=maxvarlen), parameter, public var_sfc_sss
Definition: ufo_variables_mod.F90:75
ufo_vars_mod::var_sfc_wspeed
character(len=maxvarlen), parameter, public var_sfc_wspeed
Definition: ufo_variables_mod.F90:61
ufo_vars_mod::var_sfc_emiss
character(len=maxvarlen), parameter, public var_sfc_emiss
Definition: ufo_variables_mod.F90:74
ufo_radiancecrtm_tlad_mod::ufo_radiancecrtm_simobs_ad
subroutine ufo_radiancecrtm_simobs_ad(self, geovals, obss, nvars, nlocs, hofx)
Definition: ufo_radiancecrtm_tlad_mod.F90:767
ufo_crtm_utils_mod::load_sfc_data
subroutine, public load_sfc_data(n_Profiles, n_Channels, channels, geovals, sfc, chinfo, obss, conf)
Definition: ufo_crtm_utils_mod.F90:589
ufo_vars_mod::var_lvl_transmit
character(len=maxvarlen), parameter, public var_lvl_transmit
Definition: ufo_variables_mod.F90:80
ufo_vars_mod::ufo_vars_getindex
integer function, public ufo_vars_getindex(vars, varname)
Definition: ufo_variables_mod.F90:204
ufo_crtm_utils_mod
Fortran module to provide code shared between nonlinear and tlm/adm radiance calculations.
Definition: ufo_crtm_utils_mod.F90:8
ufo_geovals_mod
Definition: ufo_geovals_mod.F90:7
ufo_radiancecrtm_tlad_mod::ufo_radiancecrtm_tlad_setup
subroutine ufo_radiancecrtm_tlad_setup(self, f_confOper, channels)
Definition: ufo_radiancecrtm_tlad_mod.F90:58
ufo_crtm_utils_mod::crtm_conf_delete
subroutine, public crtm_conf_delete(conf)
Definition: ufo_crtm_utils_mod.F90:411
ufo_constants_mod
Definition: ufo_constants_mod.F90:2
ufo_vars_mod::var_sfc_wtmp
character(len=maxvarlen), parameter, public var_sfc_wtmp
Definition: ufo_variables_mod.F90:55
ufo_crtm_utils_mod::ufo_crtm_skip_profiles
subroutine, public ufo_crtm_skip_profiles(n_Profiles, n_Channels, channels, obss, Skip_Profiles)
Definition: ufo_crtm_utils_mod.F90:450
ufo_aodcrtm_mod::varin_default
character(len=maxvarlen), dimension(5), parameter varin_default
Definition: ufo_aodcrtm_mod.F90:35
ufo_vars_mod
Definition: ufo_variables_mod.F90:8
ufo_crtm_utils_mod::crtm_conf
Definition: ufo_crtm_utils_mod.F90:57
ufo_vars_mod::var_tb
character(len=maxvarlen), parameter, public var_tb
Definition: ufo_variables_mod.F90:78
ufo_geovals_mod::ufo_geovals_get_var
subroutine, public ufo_geovals_get_var(self, varname, geoval)
Definition: ufo_geovals_mod.F90:128
ufo_constants_mod::deg2rad
real(kind_real), parameter, public deg2rad
Definition: ufo_constants_mod.F90:50
ufo_vars_mod::var_ts
character(len=maxvarlen), parameter, public var_ts
Definition: ufo_variables_mod.F90:19
ufo_vars_mod::var_pmaxlev_weightfunc
character(len=maxvarlen), parameter, public var_pmaxlev_weightfunc
Definition: ufo_variables_mod.F90:82
ufo_radiancecrtm_tlad_mod::ufo_radiancecrtm_simobs_tl
subroutine ufo_radiancecrtm_simobs_tl(self, geovals, obss, nvars, nlocs, hofx)
Definition: ufo_radiancecrtm_tlad_mod.F90:595
ufo_crtm_utils_mod::crtm_conf_setup
subroutine, public crtm_conf_setup(conf, f_confOpts, f_confOper)
Definition: ufo_crtm_utils_mod.F90:161
ufo_geovals_mod::ufo_geovals
type to hold interpolated fields required by the obs operators
Definition: ufo_geovals_mod.F90:47
ufo_vars_mod::var_opt_depth
character(len=maxvarlen), parameter, public var_opt_depth
Definition: ufo_variables_mod.F90:76
ufo_geovals_mod::ufo_geoval
type to hold interpolated field for one variable, one observation
Definition: ufo_geovals_mod.F90:40
conf
Definition: conf.py:1