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