UFO
ufo_radiancecrtm_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 radiancecrtm observations
7 
8 module ufo_radiancecrtm_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
30  private
31  character(len=MAXVARLEN), public, allocatable :: varin(:) ! variables requested from the model
32  integer, allocatable :: channels(:)
33  type(crtm_conf) :: conf
34  contains
35  procedure :: setup => ufo_radiancecrtm_setup
36  procedure :: delete => ufo_radiancecrtm_delete
37  procedure :: simobs => ufo_radiancecrtm_simobs
38  end type ufo_radiancecrtm
39 
40  character(len=maxvarlen), dimension(21), parameter :: varin_default = &
41  (/var_ts, var_prs, var_prsi, &
42  var_sfc_wfrac, var_sfc_lfrac, var_sfc_ifrac, var_sfc_sfrac, &
43  var_sfc_wtmp, var_sfc_ltmp, var_sfc_itmp, var_sfc_stmp, &
44  var_sfc_vegfrac, var_sfc_wspeed, var_sfc_wdir, var_sfc_lai, &
45  var_sfc_soilm, var_sfc_soilt, var_sfc_landtyp, &
46  var_sfc_vegtyp, var_sfc_soiltyp, var_sfc_sdepth/)
47 
48 contains
49 
50 ! ------------------------------------------------------------------------------
51 
52 subroutine ufo_radiancecrtm_setup(self, f_confOper, channels)
53 
54 implicit none
55 class(ufo_radiancecrtm), intent(inout) :: self
56 type(fckit_configuration), intent(in) :: f_confOper
57 integer(c_int), intent(in) :: channels(:) !List of channels to use
58 
59 integer :: nvars_in
60 integer :: ind, jspec
61 character(len=max_string) :: err_msg
62 type(fckit_configuration) :: f_confOpts
63 logical :: request_cldfrac
64 
65  call f_confoper%get_or_die("obs options",f_confopts)
66 
67  call crtm_conf_setup(self%conf,f_confopts,f_confoper)
68  if ( ufo_vars_getindex(self%conf%Absorbers, var_mixr) < 1 ) then
69  write(err_msg,*) 'ufo_radiancecrtm_setup error: H2O must be included in CRTM Absorbers'
70  call abor1_ftn(err_msg)
71  end if
72  if ( ufo_vars_getindex(self%conf%Absorbers, var_oz) < 1 ) then
73  write(err_msg,*) 'ufo_radiancecrtm_setup error: O3 must be included in CRTM Absorbers'
74  call abor1_ftn(err_msg)
75  end if
76 
77  ! request from the model all the hardcoded atmospheric & surface variables +
78  ! 1 * n_Absorbers
79  ! 2 * n_Clouds (mass content and effective radius)
80  nvars_in = size(varin_default) + self%conf%n_Absorbers + 2 * self%conf%n_Clouds
81  request_cldfrac = &
82  self%conf%n_Clouds > 0 .and. &
83  self%conf%Cloud_Fraction < 0.0
84  if ( request_cldfrac ) then
85  nvars_in = nvars_in + 1
86  end if
87  ! if sss is in obs options + sss
88  if (trim(self%conf%salinity_option) == "on") then
89  nvars_in = nvars_in + 1
90  end if
91 
92  allocate(self%varin(nvars_in))
93  self%varin(1:size(varin_default)) = varin_default
94  ind = size(varin_default) + 1
95  !Use list of Absorbers and Clouds from conf
96  do jspec = 1, self%conf%n_Absorbers
97  self%varin(ind) = self%conf%Absorbers(jspec)
98  ind = ind + 1
99  end do
100  do jspec = 1, self%conf%n_Clouds
101  self%varin(ind) = self%conf%Clouds(jspec,1)
102  ind = ind + 1
103  self%varin(ind) = self%conf%Clouds(jspec,2)
104  ind = ind + 1
105  end do
106  if ( request_cldfrac ) then
107  self%varin(ind) = var_cldfrac
108  ind = ind + 1
109  end if
110  if (trim(self%conf%salinity_option) == "on") then
111  self%varin(ind) = var_sfc_sss
112  ind = ind + 1
113  end if
114 
115  ! save channels
116  allocate(self%channels(size(channels)))
117  self%channels(:) = channels(:)
118 
119 end subroutine ufo_radiancecrtm_setup
120 
121 ! ------------------------------------------------------------------------------
122 
123 subroutine ufo_radiancecrtm_delete(self)
124 
125 implicit none
126 class(ufo_radiancecrtm), intent(inout) :: self
127 
128  call crtm_conf_delete(self%conf)
129 
130 end subroutine ufo_radiancecrtm_delete
131 
132 ! ------------------------------------------------------------------------------
133 
134 subroutine ufo_radiancecrtm_simobs(self, geovals, obss, nvars, nlocs, hofx, hofxdiags)
135 use fckit_mpi_module, only: fckit_mpi_comm
136 
137 implicit none
138 
139 class(ufo_radiancecrtm), intent(in) :: self !Radiance object
140 type(ufo_geovals), intent(in) :: geovals !Inputs from the model
141 integer, intent(in) :: nvars, nlocs
142 real(c_double), intent(inout) :: hofx(nvars, nlocs) !h(x) to return
143 type(ufo_geovals), intent(inout) :: hofxdiags !non-h(x) diagnostics
144 type(c_ptr), value, intent(in) :: obss !ObsSpace
145 
146 ! Local Variables
147 character(*), parameter :: PROGRAM_NAME = 'ufo_radiancecrtm_simobs'
148 character(255) :: message, version
149 character(max_string) :: err_msg
150 integer :: err_stat, alloc_stat
151 integer :: l, m, n
152 type(ufo_geoval), pointer :: temp
153 integer :: jvar, jprofile, jlevel, jchannel, ichannel, jspec
154 real(c_double) :: missing
155 type(fckit_mpi_comm) :: f_comm
156 real(kind_real) :: total_od, secant_term, wfunc_max
157 real(kind_real), allocatable :: TmpVar(:)
158 real(kind_real), allocatable :: Tao(:)
159 real(kind_real), allocatable :: Wfunc(:)
160 
161 integer :: n_Profiles, n_Layers, n_Channels
162 logical, allocatable :: Skip_Profiles(:)
163 
164 ! Define the "non-demoninational" arguments
165 type(crtm_channelinfo_type) :: chinfo(self%conf%n_Sensors)
166 type(crtm_geometry_type), allocatable :: geo(:)
167 
168 ! Define the FORWARD variables
169 type(crtm_atmosphere_type), allocatable :: atm(:)
170 type(crtm_surface_type), allocatable :: sfc(:)
171 type(crtm_rtsolution_type), allocatable :: rts(:,:)
172 
173 ! Define the K-MATRIX variables for hofxdiags
174 type(crtm_atmosphere_type), allocatable :: atm_K(:,:)
175 type(crtm_surface_type), allocatable :: sfc_K(:,:)
176 type(crtm_rtsolution_type), allocatable :: rts_K(:,:)
177 type(crtm_options_type), allocatable :: Options(:)
178 
179 ! Used to parse hofxdiags
180 character(len=MAXVARLEN) :: varstr
181 character(len=MAXVARLEN), dimension(hofxdiags%nvar) :: &
182  ystr_diags, xstr_diags
183 character(10), parameter :: jacobianstr = "_jacobian_"
184 integer :: str_pos(4), ch_diags(hofxdiags%nvar)
185 logical :: jacobian_needed
186 
187  call obsspace_get_comm(obss, f_comm)
188 
189  ! Get number of profile and layers from geovals
190  ! ---------------------------------------------
191  n_profiles = geovals%nlocs
192  call ufo_geovals_get_var(geovals, var_ts, temp)
193  n_layers = temp%nval
194  nullify(temp)
195 
196 
197  ! Program header
198  ! --------------
199  ! call CRTM_Version( Version )
200  ! call Program_Message( PROGRAM_NAME, &
201  ! 'UFO interface for the CRTM Forward and K-Matrix functions using '//&
202  ! trim(self%conf%ENDIAN_type)//' coefficient datafiles', &
203  ! 'CRTM Version: '//TRIM(Version) )
204 
205 
206  ! Initialise all the sensors at once
207  ! ----------------------------------
208  !** NOTE: CRTM_Init points to the various binary files needed for CRTM. See the
209  !** CRTM_Lifecycle.f90 for more details.
210 
211  ! write( *,'(/5x,"Initializing the CRTM...")' )
212  err_stat = crtm_init( self%conf%SENSOR_ID, chinfo, &
213  file_path=trim(self%conf%COEFFICIENT_PATH), &
214  irwatercoeff_file=trim(self%conf%IRwaterCoeff_File), &
215  irlandcoeff_file=trim(self%conf%IRlandCoeff_File), &
216  irsnowcoeff_file=trim(self%conf%IRsnowCoeff_File), &
217  iricecoeff_file=trim(self%conf%IRiceCoeff_File), &
218  viswatercoeff_file=trim(self%conf%VISwaterCoeff_File), &
219  vislandcoeff_file=trim(self%conf%VISlandCoeff_File), &
220  vissnowcoeff_file=trim(self%conf%VISsnowCoeff_File), &
221  visicecoeff_file=trim(self%conf%VISiceCoeff_File), &
222  mwwatercoeff_file=trim(self%conf%MWwaterCoeff_File), &
223  quiet=.true.)
224  message = 'Error initializing CRTM'
225  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
226 
227  ! Loop over all sensors. Not necessary if we're calling CRTM for each sensor
228  ! ----------------------------------------------------------------------------
229  sensor_loop:do n = 1, self%conf%n_Sensors
230 
231 
232  ! Pass channel list to CRTM
233  ! -------------------------
234  err_stat = crtm_channelinfo_subset(chinfo(n), self%channels, reset=.false.)
235  message = 'Error subsetting channels!'
236  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
237 
238  ! Determine the number of channels for the current sensor
239  ! -------------------------------------------------------
240  n_channels = crtm_channelinfo_n_channels(chinfo(n))
241 
242  ! Allocate the ARRAYS (for CRTM_Forward)
243  ! --------------------------------------
244  allocate( geo( n_profiles ), &
245  atm( n_profiles ), &
246  sfc( n_profiles ), &
247  rts( n_channels, n_profiles ), &
248  options( n_profiles ), &
249  stat = alloc_stat )
250  message = 'Error allocating structure arrays'
251  call crtm_comm_stat_check(alloc_stat, program_name, message, f_comm)
252 
253  if (n_layers > 0) call crtm_rtsolution_create (rts, n_layers)
254 
255  ! Create the input FORWARD structure (atm)
256  ! ----------------------------------------
257  call crtm_atmosphere_create( atm, n_layers, self%conf%n_Absorbers, self%conf%n_Clouds, self%conf%n_Aerosols )
258  if ( any(.NOT. crtm_atmosphere_associated(atm)) ) THEN
259  message = 'Error allocating CRTM Forward Atmosphere structure'
260  CALL display_message( program_name, message, failure )
261  stop
262  END IF
263 
264 
265  ! Create the input FORWARD structure (sfc)
266  ! ----------------------------------------
267  call crtm_surface_create(sfc, n_channels)
268  IF ( any(.NOT. crtm_surface_associated(sfc)) ) THEN
269  message = 'Error allocating CRTM Surface structure'
270  CALL display_message( program_name, message, failure )
271  stop
272  END IF
273 
274  CALL crtm_rtsolution_create(rts, n_layers )
275 
276  !Assign the data from the GeoVaLs
277  !--------------------------------
278  call load_atm_data(n_profiles,n_layers,geovals,atm,self%conf)
279  call load_sfc_data(n_profiles,n_channels,self%channels,geovals,sfc,chinfo,obss,self%conf)
280  call load_geom_data(obss,geo)
281 
282  ! Call THE CRTM inspection
283  ! ------------------------
284  if (self%conf%inspect > 0) then
285  call crtm_atmosphere_inspect(atm(self%conf%inspect))
286  call crtm_surface_inspect(sfc(self%conf%inspect))
287  call crtm_geometry_inspect(geo(self%conf%inspect))
288  call crtm_channelinfo_inspect(chinfo(n))
289  endif
290 
291  !! Parse hofxdiags%variables into independent/dependent variables and channel
292  !! assumed formats:
293  !! jacobian var --> <ystr>_jacobian_<xstr>_<chstr>
294  !! non-jacobian var --> <ystr>_<chstr>
295 
296  jacobian_needed = .false.
297  ch_diags = -9999
298  do jvar = 1, hofxdiags%nvar
299  varstr = hofxdiags%variables(jvar)
300  str_pos(4) = len_trim(varstr)
301  if (str_pos(4) < 1) cycle
302  str_pos(3) = index(varstr,"_",back=.true.) !final "_" before channel
303  read(varstr(str_pos(3)+1:str_pos(4)),*, err=999) ch_diags(jvar)
304  999 str_pos(1) = index(varstr,jacobianstr) - 1 !position before jacobianstr
305  if (str_pos(1) == 0) then
306  write(err_msg,*) 'ufo_radiancecrtm_simobs: _jacobian_ must be // &
307  & preceded by dependent variable in config: ', &
308  & hofxdiags%variables(jvar)
309  call abor1_ftn(err_msg)
310  else if (str_pos(1) > 0) then
311  !Diagnostic is a Jacobian member (dy/dx)
312  ystr_diags(jvar) = varstr(1:str_pos(1))
313  str_pos(2) = str_pos(1) + len(jacobianstr) + 1 !begin xstr_diags
314  jacobian_needed = .true.
315  str_pos(4) = str_pos(3) - str_pos(2)
316  xstr_diags(jvar)(1:str_pos(4)) = varstr(str_pos(2):str_pos(3)-1)
317  xstr_diags(jvar)(str_pos(4)+1:) = ""
318  else !null
319  !Diagnostic is a dependent variable (y)
320  xstr_diags(jvar) = ""
321  ystr_diags(jvar)(1:str_pos(3)-1) = varstr(1:str_pos(3)-1)
322  ystr_diags(jvar)(str_pos(3):) = ""
323  if (ch_diags(jvar) < 0) ystr_diags(jvar) = varstr
324  end if
325  end do
326 
327  allocate(skip_profiles(n_profiles))
328  call ufo_crtm_skip_profiles(n_profiles,n_channels,self%channels,obss,skip_profiles)
329  profile_loop: do jprofile = 1, n_profiles
330  options(jprofile)%Skip_Profile = skip_profiles(jprofile)
331  ! check for pressure monotonicity
332  do jlevel = atm(jprofile)%n_layers, 1, -1
333  if ( atm(jprofile)%level_pressure(jlevel) <= atm(jprofile)%level_pressure(jlevel-1) ) then
334  options(jprofile)%Skip_Profile = .true.
335  cycle profile_loop
336  end if
337  end do
338  end do profile_loop
339 
340  if (jacobian_needed) then
341  ! Allocate the ARRAYS (for CRTM_K_Matrix)
342  ! --------------------------------------
343  allocate( atm_k( n_channels, n_profiles ), &
344  sfc_k( n_channels, n_profiles ), &
345  rts_k( n_channels, n_profiles ), &
346  stat = alloc_stat )
347  message = 'Error allocating K structure arrays'
348  call crtm_comm_stat_check(alloc_stat, program_name, message, f_comm)
349 
350  ! Create output K-MATRIX structure (atm)
351  ! --------------------------------------
352  call crtm_atmosphere_create( atm_k, n_layers, self%conf%n_Absorbers, self%conf%n_Clouds, self%conf%n_Aerosols )
353  if ( any(.NOT. crtm_atmosphere_associated(atm_k)) ) THEN
354  message = 'Error allocating CRTM K-matrix Atmosphere structure (setTraj)'
355  CALL display_message( program_name, message, failure )
356  stop
357  END IF
358 
359  ! Create output K-MATRIX structure (sfc)
360  ! --------------------------------------
361  call crtm_surface_create( sfc_k, n_channels)
362  IF ( any(.NOT. crtm_surface_associated(sfc_k)) ) THEN
363  message = 'Error allocating CRTM K-matrix Surface structure (setTraj)'
364  CALL display_message( program_name, message, failure )
365  stop
366  END IF
367 
368  ! Zero the K-matrix OUTPUT structures
369  ! -----------------------------------
370  call crtm_atmosphere_zero( atm_k )
371  call crtm_surface_zero( sfc_k )
372 
373  ! Inintialize the K-matrix INPUT so that the results are dTb/dx
374  ! -------------------------------------------------------------
375  rts_k%Radiance = zero
376  rts_k%Brightness_Temperature = one
377 
378 
379  ! Call the K-matrix model
380  ! -----------------------
381  err_stat = crtm_k_matrix( atm , & ! FORWARD Input
382  sfc , & ! FORWARD Input
383  rts_k , & ! K-MATRIX Input
384  geo , & ! Input
385  chinfo(n:n) , & ! Input
386  atm_k , & ! K-MATRIX Output
387  sfc_k , & ! K-MATRIX Output
388  rts , & ! FORWARD Output
389  options ) ! Input
390  message = 'Error calling CRTM (setTraj) K-Matrix Model for '//trim(self%conf%SENSOR_ID(n))
391  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
392  else
393  ! Call the forward model call for each sensor
394  ! -------------------------------------------
395  err_stat = crtm_forward( atm , & ! Input
396  sfc , & ! Input
397  geo , & ! Input
398  chinfo(n:n) , & ! Input
399  rts , & ! Output
400  options ) ! Input
401  message = 'Error calling CRTM Forward Model for '//trim(self%conf%SENSOR_ID(n))
402  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
403  end if ! jacobian_needed
404 
405  !call CRTM_RTSolution_Inspect(rts)
406 
407  ! Put simulated brightness temperature into hofx
408  ! ----------------------------------------------
409 
410  ! Set missing value
411  missing = missing_value(missing)
412 
413  !Set to missing, then retrieve non-missing profiles
414  hofx = missing
415  do m = 1, n_profiles
416  if (.not.skip_profiles(m)) then
417  do l = 1, size(self%channels)
418  hofx(l,m) = rts(l,m)%Brightness_Temperature
419  end do
420  end if
421  end do
422 
423  ! Put simulated diagnostics into hofxdiags
424  ! ----------------------------------------------
425  do jvar = 1, hofxdiags%nvar
426  if (len(trim(hofxdiags%variables(jvar))) < 1) cycle
427 
428  if (ch_diags(jvar) > 0) then
429  if (size(pack(self%channels,self%channels==ch_diags(jvar))) /= 1) then
430  write(err_msg,*) 'ufo_radiancecrtm_simobs: mismatch between// &
431  & h(x) channels(', self%channels,') and// &
432  & ch_diags(jvar) = ', ch_diags(jvar)
433  call abor1_ftn(err_msg)
434  end if
435  end if
436 
437  jchannel = -1
438  do ichannel = 1, size(self%channels)
439  if (ch_diags(jvar) == self%channels(ichannel)) then
440  jchannel = ichannel
441  exit
442  end if
443  end do
444 
445  if (allocated(hofxdiags%geovals(jvar)%vals)) &
446  deallocate(hofxdiags%geovals(jvar)%vals)
447 
448  !============================================
449  ! Diagnostics used for QC and bias correction
450  !============================================
451  if (trim(xstr_diags(jvar)) == "") then
452  ! forward h(x) diags
453  select case(ystr_diags(jvar))
454  ! variable: optical_thickness_of_atmosphere_layer_CH
455  case (var_opt_depth)
456  hofxdiags%geovals(jvar)%nval = n_layers
457  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
458  hofxdiags%geovals(jvar)%vals = missing
459  do jprofile = 1, n_profiles
460  if (.not.skip_profiles(jprofile)) then
461  do jlevel = 1, hofxdiags%geovals(jvar)%nval
462  hofxdiags%geovals(jvar)%vals(jlevel,jprofile) = &
463  rts(jchannel,jprofile) % layer_optical_depth(jlevel)
464  end do
465  end if
466  end do
467 
468  ! variable: toa_outgoing_radiance_per_unit_wavenumber_CH [mW / (m^2 sr cm^-1)] (nval=1)
469  case (var_radiance)
470  hofxdiags%geovals(jvar)%nval = 1
471  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
472  hofxdiags%geovals(jvar)%vals = missing
473  do jprofile = 1, n_profiles
474  if (.not.skip_profiles(jprofile)) then
475  hofxdiags%geovals(jvar)%vals(1,jprofile) = &
476  rts(jchannel,jprofile) % Radiance
477  end if
478  end do
479 
480  ! variable: brightness_temperature_assuming_clear_sky_CH
481  case (var_tb_clr)
482  hofxdiags%geovals(jvar)%nval = 1
483  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
484  hofxdiags%geovals(jvar)%vals = missing
485  do jprofile = 1, n_profiles
486  if (.not.skip_profiles(jprofile)) then
487  ! Note: Using Tb_Clear requires CRTM_Atmosphere_IsFractional(cloud_coverage_flag)
488  ! to be true. For CRTM v2.3.0, that happens when
489  ! atm(jprofile)%Cloud_Fraction > MIN_COVERAGE_THRESHOLD (1e.-6)
490  hofxdiags%geovals(jvar)%vals(1,jprofile) = &
491  rts(jchannel,jprofile) % Tb_Clear
492  end if
493  end do
494 
495  ! variable: brightness_temperature_CH
496  case (var_tb)
497  hofxdiags%geovals(jvar)%nval = 1
498  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
499  hofxdiags%geovals(jvar)%vals = missing
500  do jprofile = 1, n_profiles
501  if (.not.skip_profiles(jprofile)) then
502  hofxdiags%geovals(jvar)%vals(1,jprofile) = &
503  rts(jchannel,jprofile) % Brightness_Temperature
504  end if
505  end do
506 
507  ! variable: transmittances_of_atmosphere_layer_CH
508  case (var_lvl_transmit)
509  hofxdiags%geovals(jvar)%nval = n_layers
510  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
511  hofxdiags%geovals(jvar)%vals = missing
512  allocate(tmpvar(n_profiles))
513  call obsspace_get_db(obss, "MetaData", "sensor_zenith_angle", tmpvar)
514  do jprofile = 1, n_profiles
515  if (.not.skip_profiles(jprofile)) then
516  secant_term = one/cos(tmpvar(jprofile)*deg2rad)
517  total_od = 0.0
518  do jlevel = 1, n_layers
519  total_od = total_od + rts(jchannel,jprofile) % layer_optical_depth(jlevel)
520  hofxdiags%geovals(jvar)%vals(jlevel,jprofile) = &
521  exp(-min(limit_exp,total_od*secant_term))
522  end do
523  end if
524  end do
525  deallocate(tmpvar)
526 
527  ! variable: weightingfunction_of_atmosphere_layer_CH
528  case (var_lvl_weightfunc)
529  hofxdiags%geovals(jvar)%nval = n_layers
530  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
531  hofxdiags%geovals(jvar)%vals = missing
532  allocate(tmpvar(n_profiles))
533  allocate(tao(n_layers))
534  call obsspace_get_db(obss, "MetaData", "sensor_zenith_angle", tmpvar)
535  do jprofile = 1, n_profiles
536  if (.not.skip_profiles(jprofile)) then
537  ! get layer-to-space transmittance
538  secant_term = one/cos(tmpvar(jprofile)*deg2rad)
539  total_od = 0.0
540  do jlevel = 1, n_layers
541  total_od = total_od + rts(jchannel,jprofile) % layer_optical_depth(jlevel)
542  tao(jlevel) = exp(-min(limit_exp,total_od*secant_term))
543  end do
544  ! get weighting function
545  do jlevel = n_layers-1, 1, -1
546  hofxdiags%geovals(jvar)%vals(jlevel,jprofile) = &
547  abs( (tao(jlevel+1)-tao(jlevel))/ &
548  (log(atm(jprofile)%pressure(jlevel+1))- &
549  log(atm(jprofile)%pressure(jlevel))) )
550  end do
551  hofxdiags%geovals(jvar)%vals(n_layers,jprofile) = &
552  hofxdiags%geovals(jvar)%vals(n_layers-1,jprofile)
553  end if
554  end do
555  deallocate(tmpvar)
556  deallocate(tao)
557 
558  ! variable: pressure_level_at_peak_of_weightingfunction_CH
559  case (var_pmaxlev_weightfunc)
560  hofxdiags%geovals(jvar)%nval = 1
561  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
562  hofxdiags%geovals(jvar)%vals = missing
563  allocate(tmpvar(n_profiles))
564  allocate(tao(n_layers))
565  allocate(wfunc(n_layers))
566  call obsspace_get_db(obss, "MetaData", "sensor_zenith_angle", tmpvar)
567  do jprofile = 1, n_profiles
568  if (.not.skip_profiles(jprofile)) then
569  ! get layer-to-space transmittance
570  secant_term = one/cos(tmpvar(jprofile)*deg2rad)
571  total_od = 0.0
572  do jlevel = 1, n_layers
573  total_od = total_od + rts(jchannel,jprofile) % layer_optical_depth(jlevel)
574  tao(jlevel) = exp(-min(limit_exp,total_od*secant_term))
575  end do
576  ! get weighting function
577  do jlevel = n_layers-1, 1, -1
578  wfunc(jlevel) = &
579  abs( (tao(jlevel+1)-tao(jlevel))/ &
580  (log(atm(jprofile)%pressure(jlevel+1))- &
581  log(atm(jprofile)%pressure(jlevel))) )
582  end do
583  wfunc(n_layers) = wfunc(n_layers-1)
584  ! get pressure level at the peak of the weighting function
585  wfunc_max = -999.0
586  do jlevel = n_layers-1, 1, -1
587  if (wfunc(jlevel) > wfunc_max) then
588  wfunc_max = wfunc(jlevel)
589  hofxdiags%geovals(jvar)%vals(1,jprofile) = jlevel
590  endif
591  enddo
592  end if
593  end do
594  deallocate(tmpvar)
595  deallocate(tao)
596  deallocate(wfunc)
597 
598  case default
599  write(err_msg,*) 'ufo_radiancecrtm_simobs: //&
600  & ObsDiagnostic is unsupported, ', &
601  & hofxdiags%variables(jvar)
602  ! call abor1_ftn(err_msg)
603  hofxdiags%geovals(jvar)%nval = 1
604  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
605  hofxdiags%geovals(jvar)%vals = missing
606  end select
607  else if (ystr_diags(jvar) == var_tb) then
608  ! var_tb jacobians
609  select case (xstr_diags(jvar))
610  ! variable: brightness_temperature_jacobian_air_temperature_CH
611  case (var_ts)
612  hofxdiags%geovals(jvar)%nval = n_layers
613  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
614  hofxdiags%geovals(jvar)%vals = missing
615  do jprofile = 1, n_profiles
616  if (.not.skip_profiles(jprofile)) then
617  do jlevel = 1, hofxdiags%geovals(jvar)%nval
618  hofxdiags%geovals(jvar)%vals(jlevel,jprofile) = &
619  atm_k(jchannel,jprofile) % Temperature(jlevel)
620  end do
621  end if
622  end do
623  ! variable: brightness_temperature_jacobian_humidity_mixing_ratio_CH (nval==n_Layers) --> requires MAXVARLEN=58
624  case (var_mixr)
625  hofxdiags%geovals(jvar)%nval = n_layers
626  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
627  hofxdiags%geovals(jvar)%vals = missing
628  jspec = ufo_vars_getindex(self%conf%Absorbers, var_mixr)
629  do jprofile = 1, n_profiles
630  if (.not.skip_profiles(jprofile)) then
631  do jlevel = 1, hofxdiags%geovals(jvar)%nval
632  hofxdiags%geovals(jvar)%vals(jlevel,jprofile) = &
633  atm_k(jchannel,jprofile) % Absorber(jlevel,jspec)
634  end do
635  end if
636  end do
637 
638  ! variable: brightness_temperature_jacobian_surface_temperature_CH (nval=1)
639  case (var_sfc_t)
640  hofxdiags%geovals(jvar)%nval = 1
641  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
642  hofxdiags%geovals(jvar)%vals = missing
643  do jprofile = 1, n_profiles
644  if (.not.skip_profiles(jprofile)) then
645  hofxdiags%geovals(jvar)%vals(1,jprofile) = &
646  sfc_k(jchannel,jprofile) % water_temperature &
647  + sfc_k(jchannel,jprofile) % land_temperature &
648  + sfc_k(jchannel,jprofile) % ice_temperature &
649  + sfc_k(jchannel,jprofile) % snow_temperature
650  end if
651  end do
652 
653  ! variable: brightness_temperature_jacobian_surface_emissivity_CH (nval=1)
654  case (var_sfc_emiss)
655  hofxdiags%geovals(jvar)%nval = 1
656  allocate(hofxdiags%geovals(jvar)%vals(hofxdiags%geovals(jvar)%nval,n_profiles))
657  hofxdiags%geovals(jvar)%vals = missing
658  do jprofile = 1, n_profiles
659  if (.not.skip_profiles(jprofile)) then
660  hofxdiags%geovals(jvar)%vals(1,jprofile) = &
661  rts_k(jchannel,jprofile) % surface_emissivity
662  end if
663  end do
664  case default
665  write(err_msg,*) 'ufo_radiancecrtm_simobs: //&
666  & ObsDiagnostic is unsupported, ', &
667  & hofxdiags%variables(jvar)
668  call abor1_ftn(err_msg)
669  end select
670  else
671  write(err_msg,*) 'ufo_radiancecrtm_simobs: //&
672  & ObsDiagnostic is unsupported, ', &
673  & hofxdiags%variables(jvar)
674  call abor1_ftn(err_msg)
675  end if
676  end do
677 
678  ! Deallocate the structures
679  ! -------------------------
680  call crtm_geometry_destroy(geo)
681  call crtm_atmosphere_destroy(atm)
682  call crtm_surface_destroy(sfc)
683  call crtm_rtsolution_destroy(rts)
684 
685  ! Deallocate all arrays
686  ! ---------------------
687  deallocate(geo, atm, sfc, rts, options, skip_profiles, stat = alloc_stat)
688  message = 'Error deallocating structure arrays'
689  call crtm_comm_stat_check(alloc_stat, program_name, message, f_comm)
690 
691  if (jacobian_needed) then
692  ! Deallocate the K structures
693  ! ---------------------------
694  call crtm_atmosphere_destroy(atm_k)
695  call crtm_surface_destroy(sfc_k)
696  call crtm_rtsolution_destroy(rts_k)
697 
698  ! Deallocate all K arrays
699  ! -----------------------
700  deallocate(atm_k, sfc_k, rts_k, stat = alloc_stat)
701  message = 'Error deallocating K structure arrays'
702  call crtm_comm_stat_check(alloc_stat, program_name, message, f_comm)
703  end if
704 
705  end do sensor_loop
706 
707 
708  ! Destroy CRTM instance
709  ! ---------------------
710  ! write( *, '( /5x, "Destroying the CRTM..." )' )
711  err_stat = crtm_destroy( chinfo )
712  message = 'Error destroying CRTM'
713  call crtm_comm_stat_check(err_stat, program_name, message, f_comm)
714 
715 end subroutine ufo_radiancecrtm_simobs
716 
717 ! ------------------------------------------------------------------------------
718 
719 end module ufo_radiancecrtm_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_vars_mod::var_cldfrac
character(len=maxvarlen), parameter, public var_cldfrac
Definition: ufo_variables_mod.F90:46
ufo_vars_mod::var_sfc_vegtyp
character(len=maxvarlen), parameter, public var_sfc_vegtyp
Definition: ufo_variables_mod.F90:67
ufo_vars_mod::var_sfc_vegfrac
character(len=maxvarlen), parameter, public var_sfc_vegfrac
Definition: ufo_variables_mod.F90:60
ufo_vars_mod::var_radiance
character(len=maxvarlen), parameter, public var_radiance
Definition: ufo_variables_mod.F90:77
ufo_vars_mod::var_sfc_soilm
character(len=maxvarlen), parameter, public var_sfc_soilm
Definition: ufo_variables_mod.F90:64
ufo_radiancecrtm_mod::ufo_radiancecrtm_simobs
subroutine ufo_radiancecrtm_simobs(self, geovals, obss, nvars, nlocs, hofx, hofxdiags)
Definition: ufo_radiancecrtm_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_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_vars_mod::var_sfc_ifrac
character(len=maxvarlen), parameter, public var_sfc_ifrac
Definition: ufo_variables_mod.F90:53
ufo_vars_mod::var_sfc_t
character(len=maxvarlen), parameter, public var_sfc_t
Definition: ufo_variables_mod.F90:72
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_soilt
character(len=maxvarlen), parameter, public var_sfc_soilt
Definition: ufo_variables_mod.F90:65
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_mod::ufo_radiancecrtm_setup
subroutine ufo_radiancecrtm_setup(self, f_confOper, channels)
Definition: ufo_radiancecrtm_mod.F90:53
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::var_sfc_wfrac
character(len=maxvarlen), parameter, public var_sfc_wfrac
Definition: ufo_variables_mod.F90:51
ufo_vars_mod::ufo_vars_getindex
integer function, public ufo_vars_getindex(vars, varname)
Definition: ufo_variables_mod.F90:204
ufo_vars_mod::var_mixr
character(len=maxvarlen), parameter, public var_mixr
Definition: ufo_variables_mod.F90:21
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_vars_mod::var_prsi
character(len=maxvarlen), parameter, public var_prsi
Definition: ufo_variables_mod.F90:26
ufo_vars_mod::var_sfc_itmp
character(len=maxvarlen), parameter, public var_sfc_itmp
Definition: ufo_variables_mod.F90:57
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_vars_mod::var_sfc_soiltyp
character(len=maxvarlen), parameter, public var_sfc_soiltyp
Definition: ufo_variables_mod.F90:68
ufo_vars_mod::var_sfc_sfrac
character(len=maxvarlen), parameter, public var_sfc_sfrac
Definition: ufo_variables_mod.F90:54
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_vars_mod::var_sfc_lfrac
character(len=maxvarlen), parameter, public var_sfc_lfrac
Definition: ufo_variables_mod.F90:52
ufo_aodcrtm_mod::varin_default
character(len=maxvarlen), dimension(5), parameter varin_default
Definition: ufo_aodcrtm_mod.F90:35
ufo_vars_mod::var_sfc_landtyp
character(len=maxvarlen), parameter, public var_sfc_landtyp
Definition: ufo_variables_mod.F90:66
ufo_vars_mod
Definition: ufo_variables_mod.F90:8
ufo_crtm_utils_mod::crtm_conf
Definition: ufo_crtm_utils_mod.F90:57
ufo_radiancecrtm_mod::ufo_radiancecrtm
Fortran derived type for radiancecrtm trajectory.
Definition: ufo_radiancecrtm_mod.F90:29
ufo_vars_mod::var_tb
character(len=maxvarlen), parameter, public var_tb
Definition: ufo_variables_mod.F90:78
ufo_vars_mod::var_sfc_lai
character(len=maxvarlen), parameter, public var_sfc_lai
Definition: ufo_variables_mod.F90:63
ufo_geovals_mod::ufo_geovals_get_var
subroutine, public ufo_geovals_get_var(self, varname, geoval)
Definition: ufo_geovals_mod.F90:128
ufo_vars_mod::var_tb_clr
character(len=maxvarlen), parameter, public var_tb_clr
Definition: ufo_variables_mod.F90:79
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_vars_mod::var_sfc_stmp
character(len=maxvarlen), parameter, public var_sfc_stmp
Definition: ufo_variables_mod.F90:58
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_vars_mod::var_oz
character(len=maxvarlen), parameter, public var_oz
Definition: ufo_variables_mod.F90:32
ufo_vars_mod::var_sfc_sdepth
character(len=maxvarlen), parameter, public var_sfc_sdepth
Definition: ufo_variables_mod.F90:59
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_radiancecrtm_mod
Fortran module to handle radiancecrtm observations.
Definition: ufo_radiancecrtm_mod.F90:8
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
ufo_radiancecrtm_mod::ufo_radiancecrtm_delete
subroutine ufo_radiancecrtm_delete(self)
Definition: ufo_radiancecrtm_mod.F90:124
ufo_vars_mod::var_prs
character(len=maxvarlen), parameter, public var_prs
Definition: ufo_variables_mod.F90:25
ufo_vars_mod::var_sfc_ltmp
character(len=maxvarlen), parameter, public var_sfc_ltmp
Definition: ufo_variables_mod.F90:56