FV3 Bundle
SOI_Module.f90
Go to the documentation of this file.
1 !
2 ! SOI_Module
3 !
4 ! Module containing the Successive Order of Interation (SOI) radiative
5 ! transfer solution procedures used in the CRTM.
6 !
7 !
8 ! CREATION HISTORY:
9 ! Written by: Tom Greenwald, CIMSS/SSEC; tom.greenwald@ssec.wisc.edu
10 ! Paul van Delst; paul.vandelst@noaa.gov
11 ! 20-Jan-2010
12 
13 MODULE soi_module
14 
15  ! ------------------
16  ! Environment set up
17  ! ------------------
18  ! Module use statements
19  USE type_kinds, ONLY: fp
20  USE message_handler, ONLY: success, failure, display_message
21  USE crtm_parameters, ONLY: set, zero, one, two, pi, &
22  max_n_layers, max_n_angles, max_n_legendre_terms, &
23  degrees_to_radians, &
24  secant_diffusivity, &
25  scattering_albedo_threshold, &
26  optical_depth_threshold
27  USE crtm_utility
28  USE rtv_define , ONLY: rtv_type, &
29  delta_optical_depth, &
30  max_n_doubling, &
31  max_n_soi_iterations
32  ! Disable all implicit typing
33  IMPLICIT NONE
34 
35 
36  ! --------------------
37  ! Default visibilities
38  ! --------------------
39  ! Everything private by default
40  PRIVATE
41  ! Procedures
42  PUBLIC :: crtm_soi
43  PUBLIC :: crtm_soi_tl
44  PUBLIC :: crtm_soi_ad
45  PUBLIC :: crtm_soi_version
46 
47 
48  ! -----------------
49  ! Module parameters
50  ! -----------------
51  ! Version Id for the module
52  CHARACTER(*), PARAMETER :: module_version_id = &
53  '$Id: SOI_Module.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'
54 
55 
56 CONTAINS
57 
58 
59 !################################################################################
60 !################################################################################
61 !## ##
62 !## ## PUBLIC MODULE ROUTINES ## ##
63 !## ##
64 !################################################################################
65 !################################################################################
66 
67  SUBROUTINE crtm_soi(n_Layers, & ! Input number of atmospheric layers
68  w, & ! Input layer scattering albedo
69  T_OD, & ! Input layer optical depth
70  cosmic_background, & ! Input cosmic background radiance
71  emissivity, & ! Input surface emissivity
72  reflectivity, & ! Input surface reflectivity matrix
73  Index_Sat_Angle, & ! Input satellite angle index
74  RTV) ! IN/Output upward radiance and others
75 ! ------------------------------------------------------------------------- !
76 ! !
77 ! FUNCTION: !
78 ! !
79 ! This subroutine calculates IR/MW radiance at the top of the atmosphere !
80 ! including multiple scattering using the SOI (Successive Order of !
81 ! Interaction) algorithm, which combines the Successive Order of !
82 ! Scattering and doubling methods. !
83 ! !
84 ! Written by Tom Greenwald tomg@ssec.wisc.edu !
85 ! ------------------------------------------------------------------------- !
86  IMPLICIT NONE
87  INTEGER, INTENT(IN) :: n_layers
88  REAL (fp), INTENT(IN), DIMENSION( : ) :: w, t_od
89  REAL (fp), INTENT(IN) :: cosmic_background
90  REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity
91  REAL (fp), INTENT(IN), DIMENSION( :,: ) :: reflectivity
92  INTEGER, INTENT( IN ) :: index_sat_angle
93  TYPE(rtv_type), INTENT( INOUT ) :: rtv
94 
95  ! -------------- internal variables --------------------------------- !
96 
97  INTEGER :: i, k, niter
98  REAL(fp) :: rad, rad_change
99  REAL(fp), PARAMETER :: initial_error = 1.e10
100  REAL(fp), PARAMETER :: small = 1.e-15
101  REAL(fp), PARAMETER :: sngl_scat_alb_thresh = 0.8
102  REAL(fp), PARAMETER :: opt_depth_thresh = 4.0
103  REAL(fp) :: radiance_thresh
104  REAL(fp), DIMENSION( MAX_N_ANGLES ) :: source
105 
106  ! Precompute layer R/T matrices and thermal sources
107  DO k = 1, n_layers
108  ! Precompute simple layer properties
109  rtv%e_Layer_Trans( 1 : rtv%n_Angles, k ) = exp( -t_od( k ) / rtv%COS_Angle( 1 : rtv%n_Angles ) )
110 
111  IF ( w( k ) > scattering_albedo_threshold ) THEN
112  IF ( w( k ) < sngl_scat_alb_thresh .AND. t_od( k ) < opt_depth_thresh ) THEN
113  CALL crtm_truncated_doubling( rtv%n_Streams, rtv%n_Angles, k, w( k ), t_od( k ), & ! Input
114  rtv%COS_Angle, rtv%COS_Weight, rtv%Pff( :, :, k ), & ! Input
115  rtv%Pbb( :, :, k ), rtv%Planck_Atmosphere( k ), & ! Input
116  rtv ) ! Output
117  ELSE
118  CALL crtm_doubling_layer( rtv%n_Streams, rtv%n_Angles, k, w( k ), t_od( k ), & ! Input
119  rtv%COS_Angle, rtv%COS_Weight, rtv%Pff( :, :, k ), & ! Input
120  rtv%Pbb( :, :, k ), rtv%Planck_Atmosphere( k ), & ! Input
121  rtv ) ! Output
122  END IF
123 
124 ! Subtract out direct transmission
125  DO i = 1, rtv%n_angles
126  rtv%s_Layer_Trans( i, i, k ) = rtv%s_Layer_Trans( i, i, k ) - rtv%e_layer_Trans( i, k )
127  END DO
128 
129  ELSE ! Thermal sources
130  DO i = 1, rtv%n_Angles
131  rtv%s_Layer_Source_UP( i, k ) = rtv%Planck_Atmosphere( k ) * ( one - rtv%e_Layer_Trans( i, k ) )
132  rtv%s_Layer_Source_DOWN( i, k ) = rtv%s_Layer_Source_UP( i, k )
133  END DO
134  END IF
135  END DO
136 
137  !------------------------------------
138  ! Arrays that *must* be zeroed out
139  !------------------------------------
140  rtv%s_level_Rad_UP( index_sat_angle, 0 ) = zero
141 
142  IF ( rtv%Number_SOI_Iter > 0 ) THEN
143  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, 0 : n_layers, 1 : rtv%Number_SOI_Iter ) = zero
144  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, 0 : n_layers, 1 : rtv%Number_SOI_Iter ) = zero
145  ELSE
146  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, 0 : n_layers, 1 ) = zero
147  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, 0 : n_layers, 1 ) = zero
148  END IF
149 
150  !---------------------------------
151  ! Set initial/boundary conditions
152  !---------------------------------
153  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, 0, 1 ) = cosmic_background
154  niter = 0 ! Set that we're on the zeroth order of scatter
155  rad_change = initial_error ! SOI loop uses this quantity to know when to stop iteration
156  rad = small
157 ! ACCEL = .FALSE. ! Turn on convergence acceleration
158 
159 ! IF ( ACCEL ) THEN
160 ! radiance_thresh = 5.E-7
161 ! radsum = 0.
162 ! q_old = 0.
163 ! ELSE
164  radiance_thresh = 1.e-4
165 ! END IF
166  !-----------------------------------------
167  ! This is the Order of Interaction loop
168  !-----------------------------------------
169  soi_loop: DO WHILE ( ( ( rad_change > radiance_thresh ) .AND. ( ( niter + 1 ) <= max_n_soi_iterations ) ) .OR. ( niter < 2 ) )
170 
171  niter = niter + 1 ! Note: niter = 1 is the zeroth order of interaction
172 
173  IF ( niter > 1 ) rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, 0, niter ) = zero ! No cosmic background contribution for
174  ! higher orders of interaction
175  !---------------------------------------------
176  ! Integrate down through atmospheric layers
177  !---------------------------------------------
178  layersdn_loop: DO k = 1, n_layers
179 
180  IF ( niter == 1 ) THEN
181  source( 1 : rtv%n_Angles ) = rtv%s_Layer_Source_DOWN( 1 : rtv%n_Angles, k )
182  ELSE
183  IF ( w( k ) > scattering_albedo_threshold .AND. t_od( k ) >= optical_depth_threshold ) THEN
184 ! DO j = 1, RTV%n_Angles
185 ! source( j ) = ZERO
186 ! DO i = 1, RTV%n_Angles ! Integrate over angle
187 ! source( j ) = source( j ) + RTV%s_Layer_Trans( j, i, k ) * RTV%s_Level_IterRad_DOWN( i, k - 1, niter - 1 ) &
188 ! + RTV%s_Layer_Refl( j, i, k ) * RTV%s_Level_IterRad_UP( i, k, niter - 1 )
189 ! END DO
190 ! END DO
191 
192  source( 1 : rtv%n_Angles ) = matmul( rtv%s_Layer_Trans( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
193  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, k - 1, niter - 1 ) ) + &
194  matmul( rtv%s_Layer_Refl( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
195  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, k, niter - 1 ) )
196 
197  ELSE
198  source( 1 : rtv%n_Angles ) = zero
199  END IF
200  END IF
201  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, k, niter ) = rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, k - 1, niter ) &
202  * rtv%e_Layer_Trans( 1 : rtv%n_Angles, k ) &
203  + source( 1 : rtv%n_Angles )
204 
205  END DO layersdn_loop
206 
207  !-----------------------------------------------
208  ! Account for surface reflection and emission
209  !-----------------------------------------------
210 
211  DO i = 1, rtv%n_Angles
212  rtv%s_Level_IterRad_UP( i, n_layers, niter ) = reflectivity( i, i ) * rtv%s_Level_IterRad_DOWN( i, n_layers, niter )
213  END DO
214  IF ( niter == 1 ) THEN
215  ! Add surface emission only for zeroth order of interaction
216 
217  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, n_layers, niter ) = &
218  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, n_layers, niter ) + &
219  emissivity( 1 : rtv%n_Angles ) * rtv%Planck_Surface
220 
221  END IF
222 
223  !-------------------------------------------------
224  ! Integrate back up through atmospheric layers
225  !-------------------------------------------------
226  layersup_loop: DO k = n_layers, 1, -1
227 
228  IF ( niter == 1 ) THEN
229  source( 1 :rtv%n_Angles ) = rtv%s_Layer_Source_UP( 1 : rtv%n_Angles, k )
230  ELSE
231  IF ( w( k ) > scattering_albedo_threshold .AND. t_od( k ) >= optical_depth_threshold ) THEN
232 ! DO j = 1, RTV%n_Angles
233 ! source( j ) = ZERO
234 ! DO i = 1, RTV%n_Angles ! Integrate over angle
235 ! source( j ) = source( j ) + RTV%s_Layer_Refl( j, i, k ) * RTV%s_Level_IterRad_DOWN( i, k - 1, niter - 1 ) &
236 ! + RTV%s_Layer_Trans( j, i, k ) * RTV%s_Level_IterRad_UP( i, k, niter - 1 )
237 ! END DO
238 ! END DO
239 
240  source( 1 : rtv%n_Angles ) = matmul( rtv%s_Layer_Refl( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
241  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, k - 1, niter - 1 ) ) + &
242  matmul( rtv%s_Layer_Trans( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
243  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, k, niter - 1 ) )
244 
245  ELSE
246  source( 1 : rtv%n_Angles ) = zero
247  END IF
248  END IF
249  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, k - 1, niter ) = rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, k, niter ) &
250  * rtv%e_Layer_Trans( 1 : rtv%n_Angles, k ) &
251  + source( 1 : rtv%n_Angles )
252  END DO layersup_loop
253 
254  !-------------------------------------------------------------------
255  ! Save solution at observation angle for this order of interaction
256  !-------------------------------------------------------------------
257  rad = rtv%s_Level_IterRad_UP( index_sat_angle, 0, niter )
258 
259  !---------------------------------
260  ! Add it to cumulative solution
261  !---------------------------------
262  rtv%s_level_Rad_UP( index_sat_angle, 0 ) = rtv%s_level_Rad_UP( index_sat_angle, 0 ) + rad
263 
264  !---------------------------------------
265  ! Accelerate covergence of iterations?
266  !---------------------------------------
267 ! IF ( ACCEL ) THEN
268 ! radsum = radsum + rad
269 ! IF ( niter > 1 ) THEN
270 ! diff = rad_old - rad
271 ! IF ( diff == 0. ) RETURN
272 ! q = radsum + rad ** 2 / diff
273 ! rad_change = ABS( q - q_old )
274 ! q_old = q
275 ! END IF
276 ! rad_old = rad
277 ! ELSE
278  ! Check how much cumulative solution has changed
279  rad_change = rad / rtv%s_level_Rad_UP( index_sat_angle, 0 )
280 ! END IF
281 
282  END DO soi_loop
283 
284  rtv%Number_SOI_Iter = niter
285 
286  RETURN
287  END SUBROUTINE crtm_soi
288 
289 
290  SUBROUTINE crtm_soi_tl(n_Layers, & ! Input number of atmospheric layers
291  w, & ! Input layer scattering albedo
292  T_OD, & ! Input layer optical depth
293  emissivity, & ! Input surface emissivity
294  reflectivity, & ! Input surface reflectivity matrix
295  Index_Sat_Angle, & ! Input satellite angle index
296  RTV, & ! Input structure containing forward part results
297  Planck_Atmosphere_TL, & ! Input tangent-linear atmospheric layer Planck radiance
298  Planck_Surface_TL, & ! Input TL surface Planck radiance
299  w_TL, & ! Input TL layer scattering albedo
300  T_OD_TL, & ! Input TL layer optical depth
301  emissivity_TL, & ! Input TL surface emissivity
302  reflectivity_TL, & ! Input TL reflectivity
303  Pff_TL, & ! Input TL forward phase matrix
304  Pbb_TL, & ! Input TL backward phase matrix
305  s_rad_up_TL) ! Output TL upward radiance
306 ! ------------------------------------------------------------------------- !
307 ! !
308 ! FUNCTION: !
309 ! !
310 ! This subroutine calculates the tangent linear of the IR/MW radiance at !
311 ! the top of the atmosphere including multiple scattering using the SOI !
312 ! (Successive Order of Interaction) algorithm, which combines the !
313 ! Successive Order of Scattering and doubling methods. !
314 ! !
315 ! Written by Tom Greenwald tomg@ssec.wisc.edu !
316 ! ------------------------------------------------------------------------- !
317  IMPLICIT NONE
318  INTEGER, INTENT(IN) :: n_layers
319  REAL (fp), INTENT(IN), DIMENSION( : ) :: w, t_od
320  REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity
321  REAL (fp), INTENT(IN), DIMENSION( :, : ) :: reflectivity
322  INTEGER, INTENT( IN ) :: index_sat_angle
323  TYPE(rtv_type), INTENT( IN ) :: rtv
324  REAL (fp), INTENT(IN), DIMENSION( 0: ) :: planck_atmosphere_tl
325  REAL (fp), INTENT(IN) :: planck_surface_tl
326  REAL (fp), INTENT(IN), DIMENSION( : ) :: w_tl, t_od_tl
327  REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity_tl
328  REAL (fp), INTENT(IN), DIMENSION( :, : ) :: reflectivity_tl
329  REAL (fp), INTENT(IN), DIMENSION( :, :, : ) :: pff_tl, pbb_tl
330  REAL (fp), INTENT(INOUT), DIMENSION( : ) :: s_rad_up_tl
331 
332  ! -------------- internal variables --------------------------------- !
333 
334  INTEGER :: i, k, iter
335  REAL(fp), PARAMETER :: sngl_scat_alb_thresh = 0.8
336  REAL(fp), PARAMETER :: opt_depth_thresh = 4.0
337  REAL(fp), DIMENSION( RTV%n_Angles ) :: source_tl
338  REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: e_trans_tl
339  REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: s_source_up_tl, s_source_down_tl
340  REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_trans_tl, s_refl_tl
341  REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_iterrad_up_tl, s_iterrad_down_tl
342 
343 
344  DO k = 1, n_layers
345 
346  e_trans_tl( 1 : rtv%n_Angles, k ) = -t_od_tl(k) * (exp( -t_od( k ) / rtv%COS_Angle( 1 : rtv%n_Angles ) ) ) / &
347  rtv%COS_Angle( 1 : rtv%n_Angles )
348 
349  IF ( w( k ) > scattering_albedo_threshold ) THEN
350  IF ( w( k ) < sngl_scat_alb_thresh .AND. t_od( k ) < opt_depth_thresh ) THEN
351  CALL crtm_truncated_doubling_tl(rtv%n_Streams, rtv%n_Angles, k, w( k ), t_od( k ), & !Input
352  rtv%COS_Angle( 1 : rtv%n_Angles ), rtv%COS_Weight( 1 : rtv%n_Angles ), & !Input
353  rtv%Pff( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), & !Input
354  rtv%Pbb( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), rtv%Planck_Atmosphere( k ), & !Input
355  w_tl( k ), t_od_tl( k ), pff_tl( :, :, k ), & !Input
356  pbb_tl( :, :, k ), planck_atmosphere_tl( k ), rtv, & !Input
357  s_trans_tl( :, :, k ), s_refl_tl( :, :, k ), s_source_up_tl( :, k ), & !Output
358  s_source_down_tl( :, k ) ) !Output
359  ELSE
360  CALL crtm_doubling_layer_tl(rtv%n_Streams, rtv%n_Angles, k, w( k ), t_od( k ), & !Input
361  rtv%COS_Angle( 1 : rtv%n_Angles ), rtv%COS_Weight( 1 : rtv%n_Angles ), & !Input
362  rtv%Pff( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), & !Input
363  rtv%Pbb( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), rtv%Planck_Atmosphere( k ), & !Input
364  w_tl( k ), t_od_tl( k ), pff_tl( :, :, k ), & !Input
365  pbb_tl( :, :, k ), planck_atmosphere_tl( k ), rtv, & !Input
366  s_trans_tl( :, :, k), s_refl_tl( :, :, k ), s_source_up_tl( :, k ), & !Output
367  s_source_down_tl( :, k ) ) !Output
368  END IF
369 
370 ! Subtract out direct transmission
371  DO i = 1, rtv%n_angles
372  s_trans_tl( i, i, k ) = s_trans_tl( i, i, k ) - e_trans_tl( i, k )
373  END DO
374 
375  ELSE ! Thermal sources
376  DO i = 1, rtv%n_Angles
377  s_source_up_tl( i, k ) = planck_atmosphere_tl( k ) * ( one - rtv%e_Layer_Trans( i, k ) ) - &
378  rtv%PLanck_Atmosphere( k ) * e_trans_tl( i, k )
379  s_source_down_tl( i, k ) = s_source_up_tl( i, k )
380  END DO
381  END IF
382 
383  END DO
384 
385  s_rad_up_tl( index_sat_angle ) = zero
386  s_iterrad_up_tl( 1 : rtv%n_Angles, 0 : n_layers, 1 : rtv%Number_SOI_Iter ) = zero
387  s_iterrad_down_tl( 1 : rtv%n_Angles, 0 : n_layers, 1 : rtv%Number_SOI_Iter ) = zero
388 
389  !-----------------------------------------
390  ! This is the Order of Interaction loop
391  !-----------------------------------------
392  DO iter = 1, rtv%Number_SOI_Iter
393 
394 
395  IF ( iter > 1 ) s_iterrad_down_tl( 1 : rtv%n_Angles, 0, iter ) = zero
396 
397  !---------------------------------------------
398  ! Integrate down through atmospheric layers
399  !---------------------------------------------
400  layersdn_loop: DO k = 1, n_layers
401 
402  IF ( iter == 1 ) THEN
403 
404  source_tl( 1 : rtv%n_Angles ) = s_source_down_tl( 1 : rtv%n_Angles, k )
405 
406  ELSE
407  IF ( w( k ) > scattering_albedo_threshold .AND. t_od( k ) >= optical_depth_threshold ) THEN
408  source_tl( 1 : rtv%n_Angles ) = matmul( rtv%s_Layer_Trans( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
409  s_iterrad_down_tl( 1 : rtv%n_Angles, k - 1, iter - 1 ) ) + &
410  matmul( s_trans_tl( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
411  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, k - 1, iter - 1 ) ) + &
412  matmul( rtv%s_Layer_Refl( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
413  s_iterrad_up_tl( 1 : rtv%n_Angles, k, iter - 1 ) ) + &
414  matmul( s_refl_tl( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
415  rtv%s_level_IterRad_UP( 1 : rtv%n_Angles, k, iter - 1 ) )
416  ELSE
417  source_tl( 1 : rtv%n_Angles ) = zero
418  END IF
419  END IF
420  s_iterrad_down_tl( 1 : rtv%n_Angles, k, iter ) = s_iterrad_down_tl( 1 : rtv%n_Angles, k - 1, iter ) &
421  * rtv%e_Layer_Trans( 1 : rtv%n_Angles, k ) + &
422  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, k - 1, iter ) &
423  * e_trans_tl( 1 : rtv%n_Angles, k ) &
424  + source_tl( 1 : rtv%n_Angles )
425  END DO layersdn_loop
426 
427  !-----------------------------------------------
428  ! Account for surface reflection and emission
429  !-----------------------------------------------
430 
431  DO i = 1, rtv%n_Angles
432  s_iterrad_up_tl( i, n_layers, iter ) = reflectivity_tl( i, i ) * rtv%s_Level_IterRad_DOWN( i, n_layers, iter ) + &
433  reflectivity( i, i ) * s_iterrad_down_tl( i, n_layers, iter )
434  END DO
435  IF ( iter == 1 ) THEN
436  ! Add surface emission only for zeroth order of interaction
437 
438  s_iterrad_up_tl( 1 : rtv%n_Angles, n_layers, iter ) = s_iterrad_up_tl( 1 : rtv%n_Angles, n_layers, iter ) + &
439  emissivity( 1 : rtv%n_Angles ) * planck_surface_tl + &
440  emissivity_tl( 1 : rtv%n_Angles ) * rtv%Planck_Surface
441  END IF
442 
443  !-------------------------------------------------
444  ! Integrate back up through atmospheric layers
445  !-------------------------------------------------
446  layersup_loop: DO k = n_layers, 1, -1
447 
448  IF ( iter == 1 ) THEN
449  source_tl( 1 :rtv%n_Angles ) = s_source_up_tl( 1 : rtv%n_Angles, k )
450  ELSE
451  IF ( w( k ) > scattering_albedo_threshold .AND. t_od( k ) >= optical_depth_threshold ) THEN
452  source_tl( 1 : rtv%n_Angles ) = matmul( rtv%s_Layer_Refl( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
453  s_iterrad_down_tl( 1 : rtv%n_Angles, k - 1, iter - 1 ) ) + &
454  matmul( s_refl_tl( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
455  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, k - 1, iter - 1 ) ) + &
456  matmul( rtv%s_Layer_Trans( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
457  s_iterrad_up_tl( 1 : rtv%n_Angles, k, iter - 1 ) ) + &
458  matmul( s_trans_tl( 1 : rtv%n_Angles, 1 : rtv%n_Angles, k ), &
459  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, k, iter - 1 ) )
460  ELSE
461  source_tl( 1 : rtv%n_Angles ) = zero
462  END IF
463  END IF
464  s_iterrad_up_tl( 1 : rtv%n_Angles, k - 1, iter ) = rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, k, iter ) &
465  * e_trans_tl( 1 : rtv%n_Angles, k ) + &
466  s_iterrad_up_tl( 1 : rtv%n_Angles, k, iter ) &
467  * rtv%e_Layer_Trans( 1 : rtv%n_Angles, k ) &
468  + source_tl( 1 : rtv%n_Angles )
469  END DO layersup_loop
470 
471  !----------------
472  ! Add up terms
473  !----------------
474  s_rad_up_tl( index_sat_angle ) = s_rad_up_tl( index_sat_angle ) + s_iterrad_up_tl( index_sat_angle, 0, iter )
475 
476  END DO
477 
478  RETURN
479  END SUBROUTINE crtm_soi_tl
480 
481  SUBROUTINE crtm_soi_ad(n_Layers, & ! Input number of atmospheric layers
482  w, & ! Input layer scattering albedo
483  T_OD, & ! Input layer optical depth
484  emissivity, & ! Input surface emissivity
485  reflectivity, & ! Input surface reflectivity matrix
486  Index_Sat_Angle, & ! Input satellite angle index
487  RTV, & ! Input structure containing forward results
488  s_rad_up_AD, & ! Input adjoint upward radiance
489  Planck_Atmosphere_AD, & ! Output AD atmospheric layer Planck radiance
490  Planck_Surface_AD, & ! Output AD surface Planck radiance
491  w_AD, & ! Output AD layer scattering albedo
492  T_OD_AD, & ! Output AD layer optical depth
493  emissivity_AD, & ! Output AD surface emissivity
494  reflectivity_AD, & ! Output AD surface reflectivity
495  Pff_AD, & ! Output AD forward phase matrix
496  Pbb_AD) ! Output AD backward phase matrix
497 ! ------------------------------------------------------------------------- !
498 ! FUNCTION: !
499 ! This subroutine calculates IR/MW adjoint radiance at the top of !
500 ! the atmosphere including atmospheric scattering. !
501 ! !
502 ! Tom Greenwald tomg@ssec.wisc.edu !
503 ! ------------------------------------------------------------------------- !
504  IMPLICIT NONE
505  INTEGER, INTENT(IN) :: n_layers
506  REAL (fp), INTENT(IN), DIMENSION( : ) :: w, t_od
507  REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity
508  REAL (fp), INTENT(IN), DIMENSION( :, : ) :: reflectivity
509  INTEGER, INTENT(IN) :: index_sat_angle
510  TYPE(rtv_type), INTENT(IN) :: rtv
511 
512  REAL (fp),INTENT(INOUT),DIMENSION( :, :, : ) :: pff_ad, pbb_ad
513  REAL (fp),INTENT(INOUT),DIMENSION( : ) :: w_ad,t_od_ad
514  REAL (fp),INTENT(INOUT),DIMENSION( 0: ) :: planck_atmosphere_ad
515  REAL (fp),INTENT(INOUT) :: planck_surface_ad
516  REAL (fp),INTENT(INOUT),DIMENSION( : ) :: emissivity_ad
517  REAL (fp),INTENT(INOUT),DIMENSION( :, : ) :: reflectivity_ad
518  REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_ad
519 
520 ! Local variables
521  REAL(fp), PARAMETER :: sngl_scat_alb_thresh = 0.8
522  REAL(fp), PARAMETER :: opt_depth_thresh = 4.0
523  INTEGER :: iter, k, i, j
524  REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_iterrad_up_ad
525  REAL(fp), DIMENSION( RTV%n_Angles, 0:n_Layers, RTV%Number_SOI_Iter ) :: s_iterrad_down_ad
526  REAL(fp), DIMENSION( RTV%n_Angles ) :: source_ad
527  REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: s_source_up_ad
528  REAL(fp), DIMENSION( RTV%n_Angles, n_layers ) :: s_source_down_ad
529  REAL(fp), DIMENSION( RTV%n_Angles, n_Layers ) :: e_trans_ad
530  REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_refl_ad
531  REAL(fp), DIMENSION( RTV%n_Angles, RTV%n_Angles, n_Layers ) :: s_trans_ad
532 
533 ! Zero out all local variables
534  s_iterrad_up_ad = zero
535  s_iterrad_down_ad = zero
536  source_ad = zero
537  s_source_up_ad = zero
538  s_source_down_ad = zero
539  e_trans_ad = zero
540  s_refl_ad = zero
541  s_trans_ad = zero
542  pff_ad = zero
543  pbb_ad = zero
544 
545 
546 !--------------------------------------------------------------------
547 ! Loop through each successive order of interaction in reverse order
548 !--------------------------------------------------------------------
549  DO iter = rtv%Number_SOI_Iter, 1, -1
550  s_iterrad_up_ad( index_sat_angle, 0, iter ) = s_rad_up_ad( index_sat_angle )
551 !---------------------------------------
552 ! Step down through upward integration
553 !---------------------------------------
554  DO k = 1, n_layers
555  source_ad( 1 : rtv%n_Angles ) = source_ad( 1 : rtv%n_Angles ) + s_iterrad_up_ad( 1 : rtv%n_Angles, k - 1, iter )
556  e_trans_ad( 1 : rtv%n_Angles, k ) = e_trans_ad( 1 : rtv%n_Angles, k ) + &
557  s_iterrad_up_ad( 1 : rtv%n_Angles, k - 1, iter ) * &
558  rtv%s_Level_IterRad_UP( 1 : rtv%n_Angles, k, iter )
559  s_iterrad_up_ad( 1 : rtv%n_Angles, k, iter ) = s_iterrad_up_ad( 1 : rtv%n_Angles, k, iter ) + &
560  s_iterrad_up_ad( 1 : rtv%n_Angles, k - 1, iter ) * &
561  rtv%e_Layer_Trans( 1 : rtv%n_Angles, k )
562  s_iterrad_up_ad( 1 : rtv%n_Angles, k - 1, iter ) = zero
563  IF ( iter == 1 ) THEN
564  s_source_up_ad( 1 : rtv%n_Angles, k ) = s_source_up_ad( 1 : rtv%n_Angles, k ) + source_ad( 1 : rtv%n_Angles )
565  source_ad( 1 : rtv%n_Angles ) = zero
566  ELSE
567  IF ( w( k ) > scattering_albedo_threshold .AND. t_od( k ) >= optical_depth_threshold ) THEN
568 
569  DO j = 1, rtv%n_Angles
570  DO i = 1, rtv%n_Angles ! Integrate over angle
571  s_refl_ad( j, i, k ) = s_refl_ad( j, i, k ) + source_ad( j ) * rtv%s_Level_IterRad_DOWN( i, k - 1, iter - 1 )
572  s_iterrad_down_ad( i, k - 1, iter - 1 ) = s_iterrad_down_ad( i, k - 1, iter - 1 ) + source_ad( j ) * &
573  rtv%s_Layer_Refl( j, i, k )
574  s_trans_ad( j, i, k ) = s_trans_ad( j, i, k ) + source_ad( j ) * rtv%s_Level_IterRad_UP( i, k, iter - 1 )
575  s_iterrad_up_ad( i, k, iter - 1 ) = s_iterrad_up_ad( i, k, iter - 1 ) + source_ad( j ) * &
576  rtv%s_Layer_Trans( j, i, k )
577  END DO
578  source_ad( j ) = zero
579  END DO
580  ELSE
581  source_ad( 1 : rtv%n_Angles ) = zero
582  END IF
583  END IF
584  END DO
585 
586  IF ( iter == 1 ) THEN
587  ! Add surface emission only for zeroth order of interaction
588 
589  emissivity_ad( 1 : rtv%n_Angles ) = emissivity_ad( 1 : rtv%n_Angles ) + &
590  s_iterrad_up_ad( 1 : rtv%n_Angles, n_layers, iter ) * rtv%Planck_Surface
591  planck_surface_ad = planck_surface_ad + sum( s_iterrad_up_ad( :, n_layers, iter ) * emissivity )
592  END IF
593  !-----------------------------------------------
594  ! Account for surface reflection and emission
595  !-----------------------------------------------
596  DO i = 1, rtv%n_Angles
597  reflectivity_ad( i, i ) = reflectivity_ad( i, i ) + s_iterrad_up_ad( i, n_layers, iter ) * &
598  rtv%s_Level_IterRad_DOWN( i, n_layers, iter )
599  s_iterrad_down_ad( i, n_layers, iter ) = s_iterrad_down_ad( i, n_layers, iter ) + &
600  s_iterrad_up_ad( i, n_layers, iter ) * reflectivity( i, i )
601  s_iterrad_up_ad( i, n_layers, iter ) = zero
602  END DO
603 
604 !---------------------------------------
605 ! Step up through downward integration
606 !---------------------------------------
607  DO k = n_layers, 1, -1
608  source_ad( 1 : rtv%n_Angles ) = source_ad( 1 : rtv%n_Angles ) + s_iterrad_down_ad( 1 : rtv%n_Angles, k, iter )
609  e_trans_ad( 1 : rtv%n_Angles, k ) = e_trans_ad( 1 : rtv%n_Angles, k ) + &
610  s_iterrad_down_ad( 1 : rtv%n_Angles, k, iter ) * &
611  rtv%s_Level_IterRad_DOWN( 1 : rtv%n_Angles, k - 1, iter )
612  s_iterrad_down_ad( 1 : rtv%n_Angles, k - 1, iter ) = s_iterrad_down_ad( 1 : rtv%n_Angles, k - 1, iter ) + &
613  s_iterrad_down_ad( 1 : rtv%n_Angles, k, iter ) * &
614  rtv%e_Layer_Trans( 1 : rtv%n_Angles, k )
615  s_iterrad_down_ad( 1 : rtv%n_Angles, k, iter ) = zero
616  IF ( iter == 1 ) THEN
617  s_source_down_ad( 1 : rtv%n_Angles, k ) = s_source_down_ad( 1 : rtv%n_Angles, k ) + source_ad( 1 : rtv%n_Angles )
618  source_ad( 1 : rtv%n_Angles ) = zero
619  ELSE
620  IF ( w( k ) > scattering_albedo_threshold .AND. t_od( k ) >= optical_depth_threshold ) THEN
621  DO j = 1, rtv%n_Angles
622  DO i = 1, rtv%n_Angles ! Integrate over angle
623  s_trans_ad( j, i, k ) = s_trans_ad( j, i, k ) + source_ad( j ) * rtv%s_Level_IterRad_DOWN( i, k - 1, iter - 1 )
624  s_iterrad_down_ad( i, k - 1, iter - 1 ) = s_iterrad_down_ad( i, k - 1, iter - 1 ) + source_ad( j ) * &
625  rtv%s_Layer_Trans( j, i, k )
626  s_refl_ad( j, i, k ) = s_refl_ad( j, i, k ) + source_ad( j ) * rtv%s_Level_IterRad_UP( i, k, iter - 1 )
627  s_iterrad_up_ad( i, k, iter - 1 ) = s_iterrad_up_ad( i, k, iter - 1 ) + source_ad( j ) * &
628  rtv%s_Layer_Refl( j, i, k )
629  END DO
630  source_ad( j ) = zero
631  END DO
632  ELSE
633  source_ad( 1 : rtv%n_Angles ) = zero
634  END IF
635  END IF
636  END DO
637 
638  IF ( iter > 1 ) s_iterrad_down_ad( 1 : rtv%n_Angles, 0, iter ) = zero ! No cosmic background contribution for
639  ! higher orders of interaction
640  END DO ! SOI iteration
641 
642  s_iterrad_down_ad( 1 : rtv%n_Angles, 0 : n_layers, 1 : rtv%Number_SOI_Iter ) = zero
643  s_iterrad_up_ad( 1 : rtv%n_Angles, 0 : n_layers, 1 : rtv%Number_SOI_Iter ) = zero
644  s_rad_up_ad( index_sat_angle ) = zero
645 
646  DO k = 1, n_layers
647 
648  IF ( w( k ) > scattering_albedo_threshold ) THEN
649 
650  DO i = 1, rtv%n_angles
651  e_trans_ad( i, k ) = e_trans_ad( i, k ) - s_trans_ad( i, i, k )
652  END DO
653 
654  IF ( w( k ) < sngl_scat_alb_thresh .AND. t_od( k ) < opt_depth_thresh ) THEN
655 
656  CALL crtm_truncated_doubling_ad(rtv%n_Streams, rtv%n_Angles, k, w( k ), t_od( k ), & !Input
657  rtv%COS_Angle, rtv%COS_Weight, rtv%Pff( :, :, k ), rtv%Pbb( :, :, k ), & ! Input
658  rtv%Planck_Atmosphere( k ), & !Input
659  s_trans_ad( :, :, k ), s_refl_ad( :, :, k ), s_source_up_ad( :, k ), &
660  s_source_down_ad( :, k ), rtv, w_ad( k ), t_od_ad( k ), pff_ad( :, :, k ), &
661  pbb_ad( :, :, k ), planck_atmosphere_ad( k ) ) !Output
662  ELSE
663 
664  CALL crtm_doubling_layer_ad(rtv%n_Streams, rtv%n_Angles, k, w( k ), t_od( k ), & !Input
665  rtv%COS_Angle, rtv%COS_Weight, rtv%Pff( :, :, k ), rtv%Pbb( :, :, k ), & ! Input
666  rtv%Planck_Atmosphere( k ), & !Input
667  s_trans_ad( :, :, k ), s_refl_ad( :, :, k ), s_source_up_ad( :, k ), &
668  s_source_down_ad( :, k ), rtv, w_ad( k ), t_od_ad( k ), pff_ad( :, :, k ), &
669  pbb_ad( :, :, k ), planck_atmosphere_ad( k ) ) !Output
670 
671  END IF
672 
673  ELSE ! Thermal sources
674 
675  DO i = 1, rtv%n_Angles
676  s_source_up_ad( i, k ) = s_source_up_ad( i, k ) + s_source_down_ad( i, k )
677  s_source_down_ad( i, k ) = zero
678 
679  planck_atmosphere_ad( k ) = planck_atmosphere_ad( k ) + s_source_up_ad( i, k ) * ( one - rtv%e_Layer_Trans( i, k ) )
680  e_trans_ad( i, k ) = e_trans_ad( i, k ) - rtv%Planck_Atmosphere( k ) * s_source_up_ad( i, k )
681  s_source_up_ad( i, k ) = zero
682  END DO
683  END IF
684 
685  t_od_ad( k ) = t_od_ad( k ) - sum( e_trans_ad( :, k ) * rtv%e_Layer_Trans( 1:rtv%n_Angles, k ) / &
686  rtv%COS_Angle(1:rtv%n_Angles) )
687  e_trans_ad( 1 : rtv%n_Angles, k ) = zero
688 
689  END DO
690 
691  END SUBROUTINE crtm_soi_ad
692 
693 !--------------------------------------------------------------------------------
694 !:sdoc+:
695 !
696 ! NAME:
697 ! CRTM_SOI_Version
698 !
699 ! PURPOSE:
700 ! Subroutine to return the module version information.
701 !
702 ! CALLING SEQUENCE:
703 ! CALL CRTM_SOI_Version( Id )
704 !
705 ! OUTPUT ARGUMENTS:
706 ! Id: Character string containing the version Id information
707 ! for the module.
708 ! UNITS: N/A
709 ! TYPE: CHARACTER(*)
710 ! DIMENSION: Scalar
711 ! ATTRIBUTES: INTENT(OUT)
712 !
713 !:sdoc-:
714 !--------------------------------------------------------------------------------
715 
716  SUBROUTINE crtm_soi_version( Id )
717  CHARACTER(*), INTENT(OUT) :: id
718  id = module_version_id
719  END SUBROUTINE crtm_soi_version
720 
721 
722 !################################################################################
723 !################################################################################
724 !## ##
725 !## ## PRIVATE MODULE ROUTINES ## ##
726 !## ##
727 !################################################################################
728 !################################################################################
729 
730  SUBROUTINE crtm_truncated_doubling( n_streams, & ! Input, number of streams
731  NANG, & ! Input, number of angles
732  KL, & ! Input, KL-th layer
733  single_albedo, & ! Input, single scattering albedo
734  optical_depth, & ! Input, layer optical depth
735  COS_Angle, & ! Input, COSINE of ANGLES
736  COS_Weight, & ! Input, GAUSSIAN Weights
737  ff, & ! Input, Phase matrix (forward part)
738  bb, & ! Input, Phase matrix (backward part)
739  planck_func, & ! Input, Planck for layer temperature
740  rtv) ! Output, layer transmittance, reflectance, and source
741 ! ---------------------------------------------------------------------------------------
742 ! FUNCTION
743 ! Compute layer transmission, reflection matrices and source function
744 ! at the top and bottom of the layer.
745 !
746 ! Method and References
747 ! It is a common doubling method and its theoretical basis is referred to
748 ! Hansen, J. E., 1971: Multiple scattering of polarized light in
749 ! planetary atmosphere. Part I. The doubling method, J. ATmos. Sci., 28, 120-125.
750 !
751 ! see also ADA method.
752 ! Quanhua Liu
753 ! Quanhua.Liu@noaa.gov
754 !
755 ! Modified from CRTM_Doubling_Layer by Tom Greenwald tomg@ssec.wisc.edu
756 ! ----------------------------------------------------------------------------------------
757  IMPLICIT NONE
758  INTEGER, INTENT(IN) :: n_streams, NANG, KL
759  REAL(fp), INTENT(IN) :: single_albedo, optical_depth, Planck_Func
760  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
761  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff, bb
762  TYPE(RTV_type), INTENT( INOUT ) :: RTV
763 
764  ! internal variables
765  REAL(fp), DIMENSION(NANG,NANG) :: term2,term3,trans,refl
766  REAL(fp), DIMENSION(NANG) :: C1, C2, source_up,source_down
767  REAL(fp) :: s, c
768  INTEGER :: i,j,k,L
769 
770 
771  ! Check for tiny optical depth
772 
773  IF ( optical_depth < optical_depth_threshold ) THEN
774  rtv%s_Layer_Trans(1:nang,1:nang,kl) = zero
775  DO i = 1, nang
776  rtv%s_Layer_Trans(i,i,kl) = one
777  ENDDO
778  rtv%s_Layer_Refl( 1 : nang, 1 : nang, kl ) = zero
779  rtv%s_Layer_Source_DOWN( 1 : nang, kl ) = zero
780  rtv%s_Layer_Source_UP( 1 : nang, kl ) = zero
781  RETURN
782  ENDIF
783 
784  ! -------------------------------------------------------------- !
785  ! Determine number of doubling steps and construct !
786  ! initial transmission and reflection matrix !
787  ! --------------------------------------------------------------!
788  rtv%Number_Doubling( kl ) = int( log( optical_depth / delta_optical_depth ) / log( two ) ) + 1
789  IF ( rtv%Number_Doubling( kl ) < 1 ) rtv%Number_Doubling( kl ) = 1
790  IF ( rtv%Number_Doubling( kl ) <= max_n_doubling ) THEN
791  rtv%Delta_Tau( kl ) = optical_depth / ( two**rtv%Number_Doubling( kl ) )
792  ELSE
793  rtv%Number_Doubling( kl ) = max_n_doubling
794  rtv%Delta_Tau( kl ) = delta_optical_depth
795  ENDIF
796  s = rtv%Delta_Tau( kl ) * single_albedo
797  DO i = 1, nang
798  c = s / cos_angle( i )
799  DO j = 1, nang
800  rtv%Refl( i, j, 0, kl ) = c * bb( i, j ) * cos_weight( j )
801  rtv%Trans( i, j, 0, kl ) = c * ff( i, j ) * cos_weight( j )
802  ENDDO
803  rtv%Trans( i, i, 0, kl ) = rtv%Trans( i, i, 0, kl ) + one - rtv%Delta_Tau( kl ) / cos_angle( i )
804  ENDDO
805 
806  ! -------------------------------------------------------------- !
807  ! Doubling divided sub-layers !
808  ! --------------------------------------------------------------!
809  DO l = 1, rtv%Number_Doubling( kl )
810  DO i = 1, nang
811  DO j = 1, nang
812  rtv%Inv_BeT( i, j, l, kl ) = zero
813  DO k = 1, nang
814 ! Add option to select two terms (i.e., RR+RR**2) if opd > 4 and ssa > 0.8
815  rtv%Inv_BeT( i, j, l, kl ) = rtv%Inv_BeT( i, j, l, kl ) + rtv%Refl( i, k, l - 1, kl ) * rtv%Refl( k, j, l - 1, kl )
816  ENDDO
817  ENDDO
818  rtv%Inv_BeT( i, i, l, kl ) = rtv%Inv_BeT( i, i, l, kl ) + one
819  ENDDO
820 
821  term2 = matmul( rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ), rtv%Inv_BeT( 1 : nang, 1 : nang, l, kl ) )
822  term3 = matmul( term2, rtv%Refl( 1 : nang, 1 : nang, l - 1, kl ) )
823  term3 = matmul( term3, rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) )
824 
825  rtv%Refl( 1 : nang, 1 : nang, l, kl ) = rtv%Refl( 1 : nang, 1 : nang, l - 1, kl ) + term3
826  rtv%Trans( 1 : nang, 1 : nang, l, kl ) = matmul( term2, rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) )
827 
828  ENDDO
829 
830  trans = rtv%Trans( 1 : nang, 1 : nang, rtv%Number_Doubling( kl ), kl )
831  refl = rtv%Refl( 1 : nang, 1 : nang, rtv%Number_Doubling( kl ), kl )
832 !
833 ! computing source function at up and downward directions.
834  DO i = 1, nang
835  c1( i ) = zero
836  c2( i ) = zero
837  DO j = 1, n_streams
838  c1( i ) = c1( i ) + trans( i, j )
839  c2( i ) = c2( i ) + refl( i, j )
840  ENDDO
841  IF ( i == nang .AND. nang == ( n_streams + 1 ) ) THEN
842  c1( i ) = c1( i ) + trans( nang, nang )
843  ENDIF
844  ENDDO
845 
846  DO i = 1, nang
847  source_up( i ) = ( one - c1( i ) - c2( i ) ) * planck_func
848  source_down( i ) = source_up( i )
849  ENDDO
850 
851  rtv%C1( 1 : nang, kl ) = c1
852  rtv%C2( 1 : nang, kl ) = c2
853  rtv%s_Layer_Trans( 1 : nang, 1 : nang, kl ) = trans
854  rtv%s_Layer_Refl( 1 : nang, 1 : nang, kl ) = refl
855  rtv%s_Layer_Source_DOWN( 1 : nang, kl ) = source_down
856  rtv%s_Layer_Source_UP( 1 : nang, kl ) = source_up
857 
858  RETURN
859 
860  END SUBROUTINE crtm_truncated_doubling
861 
862 
863  SUBROUTINE crtm_truncated_doubling_tl(n_streams, & ! Input, number of streams
864  NANG, & ! Input, number of angles
865  KL, & ! Input, number of angles
866  single_albedo, & ! Input, single scattering albedo
867  optical_depth, & ! Input, layer optical depth
868  COS_Angle, & ! Input, COSINE of ANGLES
869  COS_Weight, & ! Input, GAUSSIAN Weights
870  ff, & ! Input, Phase matrix (forward part)
871  bb, & ! Input, Phase matrix (backward part)
872  planck_func, & ! Input, Planck for layer temperature
873  single_albedo_tl, & ! Input, tangent-linear single albedo
874  optical_depth_tl, & ! Input, TL layer optical depth
875  ff_tl, & ! Input, TL forward Phase matrix
876  bb_tl, & ! Input, TL backward Phase matrix
877  planck_func_tl, & ! Input, TL Planck for layer temperature
878  rtv, & ! Input, structure containing forward results
879  trans_tl, & ! Output, layer tangent-linear trans
880  refl_tl, & ! Output, layer tangent-linear refl
881  source_up_tl, & ! Output, layer tangent-linear source_up
882  source_down_tl) ! Output, layer tangent-linear source_down
883 ! ---------------------------------------------------------------------------------------
884 ! FUNCTION
885 ! Compute tangent-linear layer transmission, reflection matrices and source function
886 ! at the top and bottom of the layer.
887 !
888 ! Tom Greenwald tomg@ssec.wisc.edu
889 ! ----------------------------------------------------------------------------------------
890  IMPLICIT NONE
891  INTEGER, INTENT(IN) :: n_streams, NANG, KL
892  TYPE(RTV_type), INTENT(IN) :: RTV
893  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff, bb
894  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
895  REAL(fp), INTENT(IN) :: single_albedo, optical_depth, Planck_Func
896 
897  ! Tangent-Linear Part
898  REAL(fp), INTENT(OUT), DIMENSION( :, : ) :: trans_TL, refl_TL
899  REAL(fp), INTENT(OUT), DIMENSION( : ) :: source_up_TL, source_down_TL
900  REAL(fp), INTENT(IN) :: single_albedo_TL
901  REAL(fp), INTENT(IN) :: optical_depth_TL, Planck_Func_TL
902  REAL(fp), INTENT(IN), DIMENSION( : ,: ) :: ff_TL, bb_TL
903 
904  ! internal variables
905  REAL(fp), DIMENSION( NANG, NANG ) :: term2, term3, term2_TL, term3_TL, ms_term_TL
906  REAL(fp) :: s, c
907  REAL(fp) :: s_TL, c_TL, Delta_Tau_TL
908  REAL(fp), DIMENSION( NANG ) :: C1_TL, C2_TL
909  INTEGER :: i, j, k, L
910 
911  ! Tangent-Linear Beginning
912 
913  IF( optical_depth < optical_depth_threshold ) THEN
914  trans_tl = zero
915  refl_tl = zero
916  source_up_tl = zero
917  source_down_tl = zero
918  RETURN
919  ENDIF
920 
921  s = rtv%Delta_Tau( kl ) * single_albedo
922  IF ( rtv%Number_Doubling( kl ) <= max_n_doubling ) THEN
923  delta_tau_tl = optical_depth_tl / ( two**rtv%Number_Doubling( kl ) )
924  ELSE
925  delta_tau_tl = zero
926  ENDIF
927  s_tl = delta_tau_tl * single_albedo + rtv%Delta_Tau( kl ) * single_albedo_tl
928  DO i = 1, nang
929  c = s/cos_angle( i )
930  c_tl = s_tl/cos_angle( i )
931  DO j = 1, nang
932  refl_tl( i, j ) = c_tl * bb( i, j ) * cos_weight( j ) + c * bb_tl( i, j ) * cos_weight( j )
933  trans_tl( i, j ) = c_tl * ff( i, j ) * cos_weight( j ) + c * ff_tl( i, j ) * cos_weight( j )
934  ENDDO
935  trans_tl( i, i ) = trans_tl( i, i ) - delta_tau_tl / cos_angle( i )
936  ENDDO
937 
938  DO l = 1, rtv%Number_Doubling( kl )
939  DO i = 1, nang
940  DO j = 1, nang
941  ms_term_tl( i, j ) = zero
942  DO k = 1, nang
943  ms_term_tl( i, j ) = ms_term_tl( i, j ) + rtv%Refl( i, k, l - 1, kl ) * refl_tl( k, j ) &
944  + refl_tl( i, k ) * rtv%Refl( k, j, l - 1, kl )
945  END DO
946  END DO
947  END DO
948 
949  term2 = matmul( rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ), rtv%Inv_BeT( 1 : nang, 1 : nang, l, kl ) )
950  term2_tl = matmul( rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ), ms_term_tl ) + &
951  matmul( trans_tl, rtv%Inv_BeT( 1 : nang, 1 : nang, l, kl ) )
952  term3 = matmul( term2, rtv%Refl( 1 : nang, 1 : nang, l - 1, kl ) )
953  term3_tl = matmul( term2, refl_tl ) + matmul( term2_tl, rtv%Refl( 1 : nang, 1 : nang, l - 1, kl ) )
954  term3 = matmul( term3, rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) )
955  term3_tl = matmul( term3, trans_tl ) + matmul( term3_tl, rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) )
956  refl_tl = refl_tl + term3_tl
957  trans_tl = matmul( term2, trans_tl ) + matmul( term2_tl, rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) )
958 
959  ENDDO
960 
961 ! compute TL of source function at up and downward directions.
962 
963  DO i = 1, nang
964  c1_tl( i ) = zero
965  c2_tl( i ) = zero
966  DO j = 1, n_streams
967  c1_tl( i ) = c1_tl( i ) + trans_tl( i, j )
968  c2_tl( i ) = c2_tl( i ) + refl_tl( i, j )
969  ENDDO
970  IF ( i == nang .AND. nang == ( n_streams + 1 ) ) THEN
971  c1_tl( i ) = c1_tl( i ) + trans_tl( nang, nang )
972  ENDIF
973  ENDDO
974 
975  DO i = 1, nang
976  source_up_tl( i ) = - ( c1_tl( i ) + c2_tl( i ) ) * planck_func &
977  + ( one - rtv%C1( i, kl ) - rtv%C2( i, kl ) ) * planck_func_tl
978  source_down_tl( i ) = source_up_tl( i )
979  END DO
980 
981  RETURN
982 
983  END SUBROUTINE crtm_truncated_doubling_tl
984 
985  SUBROUTINE crtm_truncated_doubling_ad(n_streams, & ! Input, number of streams
986  NANG, & ! Input, number of angles
987  KL, & ! Input, number of angles
988  single_albedo, & ! Input, single scattering albedo
989  optical_depth, & ! Input, layer optical depth
990  COS_Angle, & ! Input, COSINE of ANGLES
991  COS_Weight, & ! Input, GAUSSIAN Weights
992  ff, & ! Input, Phase matrix (forward part)
993  bb, & ! Input, Phase matrix (backward part)
994  planck_func, & ! Input, Planck for layer temperature
995  trans_ad, & ! Input, layer tangent-linear trans
996  refl_ad, & ! Input, layer tangent-linear refl
997  source_up_ad, & ! Input, layer tangent-linear source_up
998  source_down_ad, & ! Input, layer tangent-linear source_down
999  rtv, & ! Input, structure containing forward results
1000  single_albedo_ad, & ! Output adjoint single scattering albedo
1001  optical_depth_ad, & ! Output AD layer optical depth
1002  ff_ad, & ! Output AD forward Phase matrix
1003  bb_ad, & ! Output AD backward Phase matrix
1004  planck_func_ad) ! Output AD Planck for layer temperature
1005 
1006 
1007  INTEGER, INTENT(IN) :: n_streams,NANG,KL
1008  TYPE(RTV_type), INTENT(IN) :: RTV
1009  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
1010  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
1011  REAL(fp), INTENT(IN) :: single_albedo,optical_depth,Planck_Func
1012 
1013  ! Tangent-Linear Part
1014  REAL(fp), INTENT( INOUT ), DIMENSION( :,: ) :: trans_AD,refl_AD
1015  REAL(fp), INTENT( INOUT ), DIMENSION( : ) :: source_up_AD,source_down_AD
1016  REAL(fp), INTENT( INOUT ) :: single_albedo_AD
1017  REAL(fp), INTENT( INOUT ) :: optical_depth_AD,Planck_Func_AD
1018  REAL(fp), INTENT(INOUT), DIMENSION(:,:) :: ff_AD,bb_AD
1019 
1020  ! internal variables
1021  REAL(fp), DIMENSION( NANG, NANG ) :: term2, term3, term2_AD, term3_AD, ms_term_AD
1022  REAL(fp) :: s, c
1023  REAL(fp) :: s_AD, c_AD, Delta_Tau_AD
1024  REAL(fp), DIMENSION(NANG) :: C1_AD, C2_AD
1025  INTEGER :: i, j, L
1026 
1027  ! Adjoint Beginning
1028 
1029  IF ( optical_depth < optical_depth_threshold ) THEN
1030  trans_ad = zero
1031  refl_ad = zero
1032  source_up_ad = zero
1033  source_down_ad = zero
1034  RETURN
1035  ENDIF
1036 
1037 ! Remember, all output adjoint variables should be set to zero
1038 ! single_albedo_AD = 0, etc.
1039 
1040  DO i = nang, 1, -1
1041  source_up_ad( i ) = source_up_ad( i ) + source_down_ad( i )
1042  source_down_ad( i ) = zero
1043  c2_ad( i ) = -source_up_ad( i ) * planck_func
1044  c1_ad( i ) = -source_up_ad( i ) * planck_func
1045  planck_func_ad = planck_func_ad + ( one - rtv%C1( i, kl ) - rtv%C2( i, kl ) ) * source_up_ad( i )
1046  END DO
1047 
1048  ! Compute the source function in the up and downward directions.
1049  DO i = nang, 1, -1
1050 
1051  IF(i == nang .AND. nang == ( n_streams + 1 ) ) THEN
1052  trans_ad( nang, nang ) = trans_ad( nang, nang ) + c1_ad( i )
1053  ENDIF
1054 
1055  DO j = n_streams, 1, -1
1056  refl_ad( i, j ) = refl_ad( i, j ) + c2_ad( i )
1057  trans_ad( i, j ) = trans_ad( i, j ) + c1_ad( i )
1058  END DO
1059 
1060  END DO
1061 
1062  term2_ad = zero
1063  term3_ad = zero
1064  ms_term_ad = zero
1065  DO l = rtv%Number_Doubling(kl), 1, -1
1066 
1067 ! Forward calculations
1068  term2 = matmul( rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ), rtv%Inv_BeT( 1 : nang, 1 : nang, l, kl ) )
1069  term3 = matmul( term2, rtv%Refl( 1 : nang, 1 : nang, l - 1, kl ) )
1070  term3 = matmul( term3, rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) )
1071 
1072 ! Back to adjoint
1073  term2_ad = term2_ad + matmul( trans_ad, transpose( rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) ) )
1074  trans_ad = matmul( transpose( term2 ), trans_ad )
1075 
1076  term3_ad = term3_ad + refl_ad
1077 
1078  trans_ad = trans_ad + matmul( transpose( term3 ), term3_ad )
1079  term3_ad = matmul( term3_ad, transpose( rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) ) )
1080 
1081  term2_ad = term2_ad + matmul( term3_ad, transpose( rtv%Refl( 1 : nang, 1 : nang, l - 1, kl ) ) )
1082  refl_ad = refl_ad + matmul( transpose( term2 ), term3_ad )
1083  term3_ad = zero
1084 
1085  ms_term_ad = ms_term_ad + matmul( transpose( rtv%Trans( 1 : nang, 1 : nang, l - 1, kl ) ), term2_ad )
1086  trans_ad = trans_ad + matmul( term2_ad, transpose( rtv%Inv_BeT( 1 : nang, 1 : nang, l, kl ) ) )
1087  term2_ad = zero
1088 
1089  refl_ad = refl_ad + matmul( ms_term_ad, transpose( rtv%Refl( 1 : nang, 1 : nang, l - 1, kl ) ) ) + &
1090  matmul( transpose( rtv%Refl( 1 : nang, 1 : nang, l - 1, kl ) ), ms_term_ad )
1091  ms_term_ad = zero
1092 
1093  ENDDO
1094 
1095  s = rtv%Delta_Tau( kl ) * single_albedo
1096  c_ad = zero
1097  s_ad = zero
1098  delta_tau_ad = zero
1099 
1100  DO i = nang, 1, -1
1101 
1102  c = s / cos_angle( i )
1103  delta_tau_ad = delta_tau_ad - trans_ad( i, i ) / cos_angle( i )
1104 
1105  DO j = nang, 1, -1
1106  c_ad = c_ad + trans_ad( i, j ) * ff( i, j ) * cos_weight( j )
1107  ff_ad( i, j ) = ff_ad( i, j ) + trans_ad( i, j ) * c * cos_weight( j )
1108  c_ad = c_ad + refl_ad( i, j ) * bb( i, j ) * cos_weight( j )
1109  bb_ad( i, j ) = bb_ad( i, j ) + refl_ad( i, j ) * c * cos_weight( j )
1110  END DO
1111 
1112  s_ad = s_ad + c_ad / cos_angle( i )
1113  c_ad = zero
1114 
1115  ENDDO
1116 
1117  delta_tau_ad = delta_tau_ad + s_ad * single_albedo
1118  single_albedo_ad = single_albedo_ad + rtv%Delta_Tau( kl ) * s_ad
1119  optical_depth_ad = optical_depth_ad + delta_tau_ad / ( two**rtv%Number_Doubling( kl ) )
1120 
1121  END SUBROUTINE crtm_truncated_doubling_ad
1122 
1123 
1124  SUBROUTINE crtm_doubling_layer(n_streams, & ! Input, number of streams
1125  NANG, & ! Input, number of angles
1126  KL, & ! Input, KL-th layer
1127  single_albedo, & ! Input, single scattering albedo
1128  optical_depth, & ! Input, layer optical depth
1129  COS_Angle, & ! Input, COSINE of ANGLES
1130  COS_Weight, & ! Input, GAUSSIAN Weights
1131  ff, & ! Input, Phase matrix (forward part)
1132  bb, & ! Input, Phase matrix (backward part)
1133  planck_func, & ! Input, Planck for layer temperature
1134  rtv) ! Output, layer transmittance, reflectance, and source
1135 ! ---------------------------------------------------------------------------------------
1136 ! FUNCTION
1137 ! Compute layer transmission, reflection matrices and source function
1138 ! at the top and bottom of the layer.
1139 !
1140 ! Method and References
1141 ! It is a common doubling method and its theoretical basis is referred to
1142 ! Hansen, J. E., 1971: Multiple scattering of polarized light in
1143 ! planetary atmosphere. Part I. The doubling method, J. ATmos. Sci., 28, 120-125.
1144 !
1145 ! see also ADA method.
1146 ! Quanhua Liu
1147 ! Quanhua.Liu@noaa.gov
1148 ! ----------------------------------------------------------------------------------------
1149  IMPLICIT NONE
1150  INTEGER, INTENT(IN) :: n_streams,NANG,KL
1151  TYPE(RTV_type), INTENT( INOUT ) :: RTV
1152  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
1153  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
1154  REAL(fp), INTENT(IN) :: single_albedo,optical_depth,Planck_Func
1155 
1156  ! internal variables
1157  REAL(fp), DIMENSION(NANG,NANG) :: term2,term3,term4,trans,refl
1158  REAL(fp), DIMENSION(NANG) :: C1, C2, source_up,source_down
1159  REAL(fp) :: s, c
1160  INTEGER :: i,j,k,L
1161  INTEGER :: Error_Status
1162 !
1163 
1164  ! Forward part beginning
1165 
1166  IF( optical_depth < optical_depth_threshold ) THEN
1167  rtv%s_Layer_Trans(1:nang,1:nang,kl) = zero
1168  DO i = 1, nang
1169  rtv%s_Layer_Trans(i,i,kl) = one
1170  ENDDO
1171  rtv%s_Layer_Refl(1:nang,1:nang,kl) = zero
1172  rtv%s_Layer_Source_DOWN(1:nang,kl) = zero
1173  rtv%s_Layer_Source_UP(1:nang,kl) = zero
1174  RETURN
1175  ENDIF
1176 
1177  ! -------------------------------------------------------------- !
1178  ! Determining number of doubling processes and constructing !
1179  ! initial transmission and reflection matrix
1180  ! --------------------------------------------------------------!
1181  rtv%Number_Doubling(kl)=int(log(optical_depth/delta_optical_depth)/log(two))+1
1182  IF( rtv%Number_Doubling(kl) < 1 ) rtv%Number_Doubling(kl) = 1
1183  IF( rtv%Number_Doubling(kl) <= max_n_doubling ) THEN
1184  rtv%Delta_Tau(kl) = optical_depth/(two**rtv%Number_Doubling(kl))
1185  ELSE
1186  rtv%Number_Doubling(kl)=max_n_doubling
1187  rtv%Delta_Tau(kl) = delta_optical_depth
1188  ENDIF
1189  s = rtv%Delta_Tau(kl) * single_albedo
1190  DO i = 1, nang
1191  c = s/cos_angle(i)
1192  DO j = 1, nang
1193  rtv%Refl(i,j,0,kl) = c * bb(i,j) * cos_weight(j)
1194  rtv%Trans(i,j,0,kl) = c * ff(i,j) * cos_weight(j)
1195  ENDDO
1196  rtv%Trans(i,i,0,kl) = rtv%Trans(i,i,0,kl) + one - rtv%Delta_Tau(kl)/cos_angle(i)
1197  ENDDO
1198 
1199  ! -------------------------------------------------------------- !
1200  ! Doubling divided sub-layers !
1201  ! --------------------------------------------------------------!
1202  DO l = 1, rtv%Number_Doubling(kl)
1203  DO i = 1, nang
1204  DO j = 1, nang
1205  term4(i,j) = zero
1206  DO k = 1, nang
1207  term4(i,j) = term4(i,j) - rtv%Refl(i,k,l-1,kl)*rtv%Refl(k,j,l-1,kl)
1208  ENDDO
1209  ENDDO
1210  term4(i,i) = term4(i,i) + one
1211  ENDDO
1212 
1213  rtv%Inv_BeT(1:nang,1:nang,l,kl) = matinv(term4, error_status)
1214  IF( error_status /= success ) THEN
1215  print *,' error at matinv in CRTM_Doubling_layer '
1216  rtv%s_Layer_Trans(1:nang,1:nang,kl) = zero
1217  DO i = 1, nang
1218  rtv%s_Layer_Trans(i,i,kl) = exp(-optical_depth/cos_angle(i))
1219  ENDDO
1220  rtv%s_Layer_Refl(1:nang,1:nang,kl) = zero
1221  rtv%s_Layer_Source_DOWN(1:nang,kl) = zero
1222  rtv%s_Layer_Source_UP(1:nang,kl) = zero
1223  RETURN
1224  ENDIF
1225 
1226  term2 = matmul(rtv%Trans(1:nang,1:nang,l-1,kl), rtv%Inv_BeT(1:nang,1:nang,l,kl))
1227  term3 = matmul(term2, rtv%Refl(1:nang,1:nang,l-1,kl))
1228  term3 = matmul(term3, rtv%Trans(1:nang,1:nang,l-1,kl))
1229 
1230  rtv%Refl(1:nang,1:nang,l,kl) = rtv%Refl(1:nang,1:nang,l-1,kl) + term3
1231  rtv%Trans(1:nang,1:nang,l,kl) = matmul(term2, rtv%Trans(1:nang,1:nang,l-1,kl))
1232  ENDDO
1233 
1234  trans = rtv%Trans(1:nang,1:nang,rtv%Number_Doubling(kl),kl)
1235  refl = rtv%Refl(1:nang,1:nang,rtv%Number_Doubling(kl),kl)
1236 !
1237 ! computing source function at up and downward directions.
1238  DO i = 1, nang
1239  c1(i) = zero
1240  c2(i) = zero
1241  DO j = 1, n_streams
1242  c1(i) = c1(i) + trans(i,j)
1243  c2(i) = c2(i) + refl(i,j)
1244  ENDDO
1245  IF(i == nang .AND. nang == (n_streams+1)) THEN
1246  c1(i) = c1(i)+trans(nang,nang)
1247  ENDIF
1248  ENDDO
1249 
1250  DO i = 1, nang
1251  source_up(i) = (one-c1(i)-c2(i))*planck_func
1252  source_down(i) = source_up(i)
1253  ENDDO
1254 
1255  rtv%C1( 1:nang,kl ) = c1
1256  rtv%C2( 1:nang,kl ) = c2
1257  rtv%s_Layer_Trans(1:nang,1:nang,kl) = trans
1258  rtv%s_Layer_Refl(1:nang,1:nang,kl) = refl
1259  rtv%s_Layer_Source_DOWN(1:nang,kl) = source_down
1260  rtv%s_Layer_Source_UP(1:nang,kl) = source_up
1261 
1262  RETURN
1263 
1264  END SUBROUTINE crtm_doubling_layer
1265 !
1266 !
1267  SUBROUTINE crtm_doubling_layer_tl(n_streams, & ! Input, number of streams
1268  NANG, & ! Input, number of angles
1269  KL, & ! Input, number of angles
1270  single_albedo, & ! Input, single scattering albedo
1271  optical_depth, & ! Input, layer optical depth
1272  COS_Angle, & ! Input, COSINE of ANGLES
1273  COS_Weight, & ! Input, GAUSSIAN Weights
1274  ff, & ! Input, Phase matrix (forward part)
1275  bb, & ! Input, Phase matrix (backward part)
1276  planck_func, & ! Input, Planck for layer temperature
1277  single_albedo_tl, & ! Input, tangent-linear single albedo
1278  optical_depth_tl, & ! Input, TL layer optical depth
1279  ff_tl, & ! Input, TL forward Phase matrix
1280  bb_tl, & ! Input, TL backward Phase matrix
1281  planck_func_tl, & ! Input, TL Planck for layer temperature
1282  rtv, & ! Input, structure containing forward results
1283  trans_tl, & ! Output, layer tangent-linear trans
1284  refl_tl, & ! Output, layer tangent-linear refl
1285  source_up_tl, & ! Output, layer tangent-linear source_up
1286  source_down_tl) ! Output, layer tangent-linear source_down
1287 ! ---------------------------------------------------------------------------------------
1288 ! FUNCTION
1289 ! Compute tangent-linear layer transmission, reflection matrices and source function
1290 ! at the top and bottom of the layer.
1291 !
1292 ! Quanhua Liu
1293 ! Quanhua.Liu@noaa.gov
1294 ! ----------------------------------------------------------------------------------------
1295  IMPLICIT NONE
1296  INTEGER, INTENT(IN) :: n_streams,NANG,KL
1297  TYPE(RTV_type), INTENT(IN) :: RTV
1298  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
1299  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
1300  REAL(fp), INTENT(IN) :: single_albedo,optical_depth,Planck_Func
1301 
1302  ! Tangent-Linear Part
1303  REAL(fp), INTENT(OUT), DIMENSION( :,: ) :: trans_TL,refl_TL
1304  REAL(fp), INTENT(OUT), DIMENSION( : ) :: source_up_TL,source_down_TL
1305  REAL(fp), INTENT(IN) :: single_albedo_TL
1306  REAL(fp), INTENT(IN) :: optical_depth_TL,Planck_Func_TL
1307  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff_TL,bb_TL
1308 
1309  ! internal variables
1310  REAL(fp), DIMENSION(NANG,NANG) :: term1,term2,term3,term4,term5_TL
1311  REAL(fp) :: s, c
1312  REAL(fp) :: s_TL, c_TL, Delta_Tau_TL
1313  REAL(fp), DIMENSION(NANG) :: C1_TL, C2_TL
1314  INTEGER :: i,j,L
1315 !
1316  ! Tangent-Linear Beginning
1317 
1318  IF( optical_depth < optical_depth_threshold ) THEN
1319  trans_tl = zero
1320  refl_tl = zero
1321  source_up_tl = zero
1322  source_down_tl = zero
1323  RETURN
1324  ENDIF
1325 
1326  s = rtv%Delta_Tau(kl) * single_albedo
1327  IF ( rtv%Number_Doubling( kl ) <= max_n_doubling ) THEN
1328  delta_tau_tl = optical_depth_tl / (two**rtv%Number_Doubling( kl ) )
1329  ELSE
1330  delta_tau_tl = zero
1331  ENDIF
1332  s_tl = delta_tau_tl * single_albedo + rtv%Delta_Tau(kl) * single_albedo_tl
1333  DO i = 1, nang
1334  c = s/cos_angle(i)
1335  c_tl = s_tl/cos_angle(i)
1336  DO j = 1, nang
1337  refl_tl(i,j) = c_tl*bb(i,j)*cos_weight(j)+c*bb_tl(i,j)*cos_weight(j)
1338  trans_tl(i,j) =c_tl*ff(i,j)*cos_weight(j)+c*ff_tl(i,j)*cos_weight(j)
1339  ENDDO
1340  trans_tl(i,i) =trans_tl(i,i) - delta_tau_tl/cos_angle(i)
1341  ENDDO
1342 
1343  DO l = 1, rtv%Number_Doubling(kl)
1344 
1345  term1 = matmul(rtv%Trans(1:nang,1:nang,l-1,kl),rtv%Inv_BeT(1:nang,1:nang,l,kl))
1346  term2 = matmul(rtv%Inv_BeT(1:nang,1:nang,l,kl),rtv%Refl(1:nang,1:nang,l-1,kl))
1347  term3 = matmul(rtv%Inv_BeT(1:nang,1:nang,l,kl),rtv%Trans(1:nang,1:nang,l-1,kl))
1348  term4 = matmul(term2,rtv%Trans(1:nang,1:nang,l-1,kl))
1349  term5_tl = matmul(refl_tl,rtv%Refl(1:nang,1:nang,l-1,kl)) &
1350  + matmul(rtv%Refl(1:nang,1:nang,l-1,kl),refl_tl)
1351 
1352  refl_tl=refl_tl+matmul(matmul(term1,term5_tl),term4)+matmul(trans_tl,term4) &
1353  +matmul(matmul(term1,refl_tl),rtv%Trans(1:nang,1:nang,l-1,kl))
1354  refl_tl=refl_tl+matmul(matmul(term1,rtv%Refl(1:nang,1:nang,l-1,kl)),trans_tl)
1355 
1356  trans_tl=matmul(trans_tl,term3) &
1357  +matmul(matmul(term1,term5_tl),term3)+matmul(term1,trans_tl)
1358 
1359  ENDDO
1360 
1361 ! computing source function at up and downward directions.
1362 
1363  DO i = 1, nang
1364  c1_tl(i) = zero
1365  c2_tl(i) = zero
1366  DO j = 1, n_streams
1367  c1_tl(i) = c1_tl(i) + trans_tl(i,j)
1368  c2_tl(i) = c2_tl(i) + refl_tl(i,j)
1369  ENDDO
1370  IF(i == nang .AND. nang == (n_streams+1)) THEN
1371  c1_tl(i) = c1_tl(i)+trans_tl(nang,nang)
1372  ENDIF
1373  ENDDO
1374 
1375  DO i = 1, nang
1376  source_up_tl(i) = -(c1_tl(i)+c2_tl(i))*planck_func &
1377  + (one-rtv%C1(i,kl)-rtv%C2(i,kl))*planck_func_tl
1378  source_down_tl(i) = source_up_tl(i)
1379  END DO
1380 
1381  RETURN
1382 
1383  END SUBROUTINE crtm_doubling_layer_tl
1384 
1385 
1386 ! ---------------------------------------------------------------------------------------
1387 ! FUNCTION
1388 ! Compute layer adjoint transmission, reflection matrices and source function
1389 ! at the top and bottom of the layer.
1390 !
1391 ! Quanhua Liu
1392 ! Quanhua.Liu@noaa.gov
1393 ! ----------------------------------------------------------------------------------------
1394 
1395  SUBROUTINE crtm_doubling_layer_ad(n_streams, & ! Input, number of streams
1396  NANG, & ! Input, number of angles
1397  KL, & ! Input, number of angles
1398  single_albedo, & ! Input, single scattering albedo
1399  optical_depth, & ! Input, layer optical depth
1400  COS_Angle, & ! Input, COSINE of ANGLES
1401  COS_Weight, & ! Input, GAUSSIAN Weights
1402  ff, & ! Input, Phase matrix (forward part)
1403  bb, & ! Input, Phase matrix (backward part)
1404  planck_func, & ! Input, Planck for layer temperature
1405  trans_ad, & ! Input, layer tangent-linear trans
1406  refl_ad, & ! Input, layer tangent-linear refl
1407  source_up_ad, & ! Input, layer tangent-linear source_up
1408  source_down_ad, & ! Input, layer tangent-linear source_down
1409  rtv, & ! Input, structure containing forward results
1410  single_albedo_ad, & ! Output adjoint single scattering albedo
1411  optical_depth_ad, & ! Output AD layer optical depth
1412  ff_ad, & ! Output AD forward Phase matrix
1413  bb_ad, & ! Output AD backward Phase matrix
1414  planck_func_ad) ! Output AD Planck for layer temperature
1415 
1416 
1417  INTEGER, INTENT(IN) :: n_streams,NANG,KL
1418  TYPE(RTV_type), INTENT(IN) :: RTV
1419  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
1420  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
1421  REAL(fp), INTENT(IN) :: single_albedo,optical_depth,Planck_Func
1422 
1423  ! Tangent-Linear Part
1424  REAL(fp), INTENT( INOUT ), DIMENSION( :,: ) :: trans_AD,refl_AD
1425  REAL(fp), INTENT( INOUT ), DIMENSION( : ) :: source_up_AD,source_down_AD
1426  REAL(fp), INTENT( INOUT ) :: single_albedo_AD
1427  REAL(fp), INTENT( INOUT ) :: optical_depth_AD,Planck_Func_AD
1428  REAL(fp), INTENT(INOUT), DIMENSION(:,:) :: ff_AD,bb_AD
1429 
1430  ! internal variables
1431  REAL(fp), DIMENSION(NANG,NANG) :: term1,term2,term3,term4,term5_AD
1432  REAL(fp) :: s, c
1433  REAL(fp) :: s_AD, c_AD, Delta_Tau_AD
1434  REAL(fp), DIMENSION(NANG) :: C1_AD, C2_AD
1435  INTEGER :: i,j,L
1436 
1437  ! Tangent-Linear Beginning
1438 
1439  IF( optical_depth < optical_depth_threshold ) THEN
1440  trans_ad = zero
1441  refl_ad = zero
1442  source_up_ad = zero
1443  source_down_ad = zero
1444  RETURN
1445  ENDIF
1446 
1447  DO i = nang, 1, -1
1448  source_up_ad(i) = source_up_ad(i) + source_down_ad(i)
1449  source_down_ad(i) = zero
1450  c2_ad(i) = -source_up_ad(i)*planck_func
1451  c1_ad(i) = -source_up_ad(i)*planck_func
1452  planck_func_ad = planck_func_ad + (one-rtv%C1(i,kl)-rtv%C2(i,kl))*source_up_ad(i)
1453  END DO
1454 
1455  ! Compute the source function in the up and downward directions.
1456  DO i = nang, 1, -1
1457 
1458  IF(i == nang .AND. nang == (n_streams+1)) THEN
1459  trans_ad(nang,nang)=trans_ad(nang,nang)+c1_ad(i)
1460  ENDIF
1461 
1462  DO j = n_streams, 1, -1
1463  refl_ad(i,j)=refl_ad(i,j)+c2_ad(i)
1464  trans_ad(i,j)=trans_ad(i,j)+c1_ad(i)
1465  END DO
1466 
1467  END DO
1468 
1469  DO l = rtv%Number_Doubling(kl), 1, -1
1470 
1471  term1 = matmul(rtv%Trans(1:nang,1:nang,l-1,kl),rtv%Inv_BeT(1:nang,1:nang,l,kl))
1472  term2 = matmul(rtv%Inv_BeT(1:nang,1:nang,l,kl),rtv%Refl(1:nang,1:nang,l-1,kl))
1473  term3 = matmul(rtv%Inv_BeT(1:nang,1:nang,l,kl),rtv%Trans(1:nang,1:nang,l-1,kl))
1474  term4 = matmul(term2,rtv%Trans(1:nang,1:nang,l-1,kl))
1475 
1476  term5_ad = matmul(matmul(transpose(term1),trans_ad),transpose(term3))
1477  trans_ad = matmul(trans_ad,transpose(term3))+matmul(transpose(term1),trans_ad)
1478 
1479  trans_ad=trans_ad+matmul(transpose(matmul(term1,rtv%Refl(1:nang,1:nang,l-1,kl))),refl_ad)
1480 
1481  term5_ad =term5_ad+matmul(matmul(transpose(term1),refl_ad),transpose(term4))
1482  trans_ad = trans_ad+matmul(refl_ad,transpose(term4))
1483  refl_ad = refl_ad+matmul(matmul(transpose(term1),refl_ad),transpose(rtv%Trans(1:nang,1:nang,l-1,kl)))
1484  refl_ad = refl_ad+matmul(term5_ad,transpose(rtv%Refl(1:nang,1:nang,l-1,kl)))
1485  refl_ad = refl_ad+matmul(transpose(rtv%Refl(1:nang,1:nang,l-1,kl)),term5_ad)
1486 
1487  ENDDO
1488 
1489  s = rtv%Delta_Tau(kl) * single_albedo
1490  c_ad = zero
1491  s_ad = zero
1492  delta_tau_ad=zero
1493 
1494  DO i = nang, 1, -1
1495 
1496  c = s/cos_angle(i)
1497  delta_tau_ad = delta_tau_ad - trans_ad(i,i)/cos_angle(i)
1498 
1499  DO j = nang, 1, -1
1500  c_ad = c_ad + trans_ad(i,j)*ff(i,j)*cos_weight(j)
1501  ff_ad(i,j)=ff_ad(i,j)+trans_ad(i,j)*c*cos_weight(j)
1502  c_ad = c_ad + refl_ad(i,j)*bb(i,j)*cos_weight(j)
1503  bb_ad(i,j)=bb_ad(i,j) + refl_ad(i,j)*c*cos_weight(j)
1504  END DO
1505 
1506  s_ad = s_ad + c_ad/cos_angle(i)
1507  c_ad = zero
1508 
1509  ENDDO
1510 
1511  delta_tau_ad = delta_tau_ad + s_ad* single_albedo
1512  single_albedo_ad = single_albedo_ad+rtv%Delta_Tau(kl) * s_ad
1513  optical_depth_ad = optical_depth_ad + delta_tau_ad/(two**rtv%Number_Doubling(kl))
1514 
1515  END SUBROUTINE crtm_doubling_layer_ad
1516 
1517 END MODULE soi_module
1518 
soi_module::crtm_doubling_layer_tl
subroutine crtm_doubling_layer_tl(n_streams, NANG, KL, single_albedo, optical_depth, COS_Angle, COS_Weight, ff, bb, Planck_Func, single_albedo_TL, optical_depth_TL, ff_TL, bb_TL, Planck_Func_TL, RTV, trans_TL, refl_TL, source_up_TL, source_down_TL)
Definition: SOI_Module.f90:1287
message_handler::failure
integer, parameter, public failure
Definition: Message_Handler.f90:30
crtm_utility
Definition: CRTM_Utility.f90:17
crtm_parameters::set
integer, parameter, public set
Definition: CRTM_Parameters.f90:113
crtm_parameters::zero
real(fp), parameter, public zero
Definition: CRTM_Parameters.f90:70
type_kinds::fp
integer, parameter, public fp
Definition: Type_Kinds.f90:124
crtm_parameters::max_n_angles
integer, parameter, public max_n_angles
Definition: CRTM_Parameters.f90:261
crtm_parameters::scattering_albedo_threshold
real(fp), parameter, public scattering_albedo_threshold
Definition: CRTM_Parameters.f90:255
rtv_define::delta_optical_depth
real(fp), parameter, public delta_optical_depth
Definition: RTV_Define.f90:78
crtm_parameters
Definition: CRTM_Parameters.f90:12
crtm_parameters::optical_depth_threshold
real(fp), parameter, public optical_depth_threshold
Definition: CRTM_Parameters.f90:252
soi_module::crtm_soi
subroutine, public crtm_soi(n_Layers, w, T_OD, cosmic_background, emissivity, reflectivity, Index_Sat_Angle, RTV)
Definition: SOI_Module.f90:75
crtm_parameters::max_n_legendre_terms
integer, parameter, public max_n_legendre_terms
Definition: CRTM_Parameters.f90:258
rtv_define::max_n_soi_iterations
integer, parameter, public max_n_soi_iterations
Definition: RTV_Define.f90:89
message_handler
Definition: Message_Handler.f90:4
crtm_parameters::secant_diffusivity
real(fp), parameter, public secant_diffusivity
Definition: CRTM_Parameters.f90:227
crtm_utility::module_version_id
character(*), parameter module_version_id
Definition: CRTM_Utility.f90:58
soi_module::crtm_soi_version
subroutine, public crtm_soi_version(Id)
Definition: SOI_Module.f90:717
crtm_parameters::one
real(fp), parameter, public one
Definition: CRTM_Parameters.f90:71
message_handler::display_message
recursive subroutine, public display_message(Routine_Name, Message, Error_State, Message_Log)
Definition: Message_Handler.f90:118
soi_module::crtm_doubling_layer
subroutine crtm_doubling_layer(n_streams, NANG, KL, single_albedo, optical_depth, COS_Angle, COS_Weight, ff, bb, Planck_Func, RTV)
Definition: SOI_Module.f90:1135
soi_module::crtm_truncated_doubling_ad
subroutine crtm_truncated_doubling_ad(n_streams, NANG, KL, single_albedo, optical_depth, COS_Angle, COS_Weight, ff, bb, Planck_Func, trans_AD, refl_AD, source_up_AD, source_down_AD, RTV, single_albedo_AD, optical_depth_AD, ff_AD, bb_AD, Planck_Func_AD)
Definition: SOI_Module.f90:1005
crtm_parameters::two
real(fp), parameter, public two
Definition: CRTM_Parameters.f90:72
soi_module::crtm_truncated_doubling_tl
subroutine crtm_truncated_doubling_tl(n_streams, NANG, KL, single_albedo, optical_depth, COS_Angle, COS_Weight, ff, bb, Planck_Func, single_albedo_TL, optical_depth_TL, ff_TL, bb_TL, Planck_Func_TL, RTV, trans_TL, refl_TL, source_up_TL, source_down_TL)
Definition: SOI_Module.f90:883
crtm_parameters::degrees_to_radians
real(fp), parameter, public degrees_to_radians
Definition: CRTM_Parameters.f90:88
rtv_define::rtv_type
Definition: RTV_Define.f90:106
soi_module::crtm_soi_tl
subroutine, public crtm_soi_tl(n_Layers, w, T_OD, emissivity, reflectivity, Index_Sat_Angle, RTV, Planck_Atmosphere_TL, Planck_Surface_TL, w_TL, T_OD_TL, emissivity_TL, reflectivity_TL, Pff_TL, Pbb_TL, s_rad_up_TL)
Definition: SOI_Module.f90:306
crtm_parameters::max_n_layers
integer, parameter, public max_n_layers
Definition: CRTM_Parameters.f90:64
rtv_define::max_n_doubling
integer, parameter, public max_n_doubling
Definition: RTV_Define.f90:86
soi_module::crtm_soi_ad
subroutine, public crtm_soi_ad(n_Layers, w, T_OD, emissivity, reflectivity, Index_Sat_Angle, RTV, s_rad_up_AD, Planck_Atmosphere_AD, Planck_Surface_AD, w_AD, T_OD_AD, emissivity_AD, reflectivity_AD, Pff_AD, Pbb_AD)
Definition: SOI_Module.f90:497
type_kinds
Definition: Type_Kinds.f90:12
soi_module::crtm_truncated_doubling
subroutine crtm_truncated_doubling(n_streams, NANG, KL, single_albedo, optical_depth, COS_Angle, COS_Weight, ff, bb, Planck_Func, RTV)
Definition: SOI_Module.f90:741
soi_module
Definition: SOI_Module.f90:13
soi_module::crtm_doubling_layer_ad
subroutine crtm_doubling_layer_ad(n_streams, NANG, KL, single_albedo, optical_depth, COS_Angle, COS_Weight, ff, bb, Planck_Func, trans_AD, refl_AD, source_up_AD, source_down_AD, RTV, single_albedo_AD, optical_depth_AD, ff_AD, bb_AD, Planck_Func_AD)
Definition: SOI_Module.f90:1415
message_handler::success
integer, parameter, public success
Definition: Message_Handler.f90:27
crtm_utility::matinv
real(fp) function, dimension(size(a, 1), size(a, 2)), public matinv(A, Error_Status)
Definition: CRTM_Utility.f90:106
rtv_define
Definition: RTV_Define.f90:21
crtm_parameters::pi
real(fp), parameter, public pi
Definition: CRTM_Parameters.f90:86