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