FV3 Bundle
ADA_Module.f90
Go to the documentation of this file.
1 !
2 ! ADA_Module
3 !
4 ! Module containing the Adding Doubling Adding (ADA) radiative
5 ! transfer solution procedures used in the CRTM.
6 !
7 !
8 ! CREATION HISTORY:
9 ! Written by: Quanhua Liu, QSS at JCSDA; quanhua.liu@noaa.gov
10 ! Yong Han, NOAA/NESDIS; yong.han@noaa.gov
11 ! Paul van Delst; CIMMS/SSEC; paul.vandelst@noaa.gov
12 ! 08-Jun-2004
13 
14 MODULE ada_module
15 
16  ! ------------------
17  ! Environemnt set up
18  ! ------------------
19  ! Module use statements
20  USE rtv_define
21  USE crtm_parameters
22  USE type_kinds
23  USE message_handler
24  USE crtm_utility
25 
26  IMPLICIT NONE
27 
28  ! --------------------
29  ! Default visibilities
30  ! --------------------
31  ! Everything private by default
32  PRIVATE
33 
34  PUBLIC crtm_ada
35  PUBLIC crtm_ada_tl
36  PUBLIC crtm_ada_ad
37 
38  ! -----------------
39  ! Module parameters
40  ! -----------------
41  ! Version Id for the module
42  CHARACTER(*), PARAMETER :: module_version_id = &
43  '$Id: $'
44 
45 CONTAINS
46 
47 !################################################################################
48 !################################################################################
49 !## ##
50 !## ## PUBLIC MODULE ROUTINES ## ##
51 !## ##
52 !################################################################################
53 !################################################################################
54 
55 
56  SUBROUTINE crtm_ada(n_Layers, & ! Input number of atmospheric layers
57  w, & ! Input layer scattering albedo
58  T_OD, & ! Input layer optical depth
59  cosmic_background, & ! Input cosmic background radiance
60  emissivity, & ! Input surface emissivity
61  reflectivity, & ! Input surface reflectivity matrix
62  direct_reflectivity, & ! Input surface direct reflectivity
63  RTV) ! IN/Output upward radiance and others
64 ! ------------------------------------------------------------------------- !
65 ! FUNCTION: !
66 ! This subroutine calculates IR/MW radiance at the top of the atmosphere !
67 ! including atmospheric scattering. The scheme will include solar part. !
68 ! The ADA algorithm computes layer reflectance and transmittance as well !
69 ! as source function by the subroutine CRTM_Doubling_layer, then uses !
70 ! an adding method to integrate the layer and surface components. !
71 ! !
72 ! Quanhua Liu Quanhua.Liu@noaa.gov !
73 ! ------------------------------------------------------------------------- !
74  IMPLICIT NONE
75  INTEGER, INTENT(IN) :: n_layers
76  INTEGER nz
77  TYPE(rtv_type), INTENT( INOUT ) :: rtv
78  REAL (fp), INTENT(IN), DIMENSION( : ) :: w,t_od
79  REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity, direct_reflectivity
80  REAL (fp), INTENT(IN), DIMENSION( :,: ) :: reflectivity
81  REAL (fp), INTENT(IN) :: cosmic_background
82 
83  ! -------------- internal variables --------------------------------- !
84  ! Abbreviations: !
85  ! s: scattering, rad: radiance, trans: transmission, !
86  ! refl: reflection, up: upward, down: downward !
87  ! --------------------------------------------------------------------!
88  REAL (fp), DIMENSION(RTV%n_Angles, RTV%n_Angles) :: temporal_matrix
89  REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down
90  REAL (fp), DIMENSION(0:n_Layers) :: total_opt
91  INTEGER :: i, j, k, error_status
92  CHARACTER(*), PARAMETER :: routine_name = 'CRTM_ADA'
93  CHARACTER(256) :: message
94 
95  total_opt(0) = zero
96  DO k = 1, n_layers
97  total_opt(k) = total_opt(k-1) + t_od(k)
98  END DO
99 
100  nz = rtv%n_Angles
101  rtv%s_Layer_Trans = zero
102  rtv%s_Layer_Refl = zero
103  rtv%s_Level_Refl_UP = zero
104  rtv%s_Level_Rad_UP = zero
105  rtv%s_Layer_Source_UP = zero
106  rtv%s_Layer_Source_DOWN = zero
107 
108  rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,n_layers)=reflectivity(1:rtv%n_Angles,1:rtv%n_Angles)
109 
110  IF( rtv%mth_Azi == 0 ) THEN
111  rtv%s_Level_Rad_UP(1:rtv%n_Angles,n_layers ) = emissivity(1:rtv%n_Angles)*rtv%Planck_Surface
112  END IF
113 
114  IF( rtv%Solar_Flag_true ) THEN
115  rtv%s_Level_Rad_UP(1:nz,n_layers ) = rtv%s_Level_Rad_UP(1:nz,n_layers )+direct_reflectivity(1:nz)* &
116  rtv%COS_SUN*rtv%Solar_irradiance/pi*exp(-total_opt(n_layers)/rtv%COS_SUN)
117  END IF
118 
119  ! UPWARD ADDING LOOP STARTS FROM BOTTOM LAYER TO ATMOSPHERIC TOP LAYER.
120  DO 10 k = n_layers, 1, -1
121 
122  ! Compute tranmission and reflection matrices for a layer
123  IF(w(k) > scattering_albedo_threshold) THEN
124 
125  ! ----------------------------------------------------------- !
126  ! CALL multiple-stream algorithm for computing layer !
127  ! transmission, reflection, and source functions. !
128  ! ----------------------------------------------------------- !
129 
130  CALL crtm_amom_layer( &
131  rtv%n_Streams, &
132  rtv%n_Angles,k,w(k), &
133  t_od(k), &
134  total_opt(k-1), &
135  rtv%COS_Angle, & ! Input
136  rtv%COS_Weight, &
137  rtv%Pff(:,:,k), &
138  rtv%Pbb(:,:,k), & ! Input
139  rtv%Planck_Atmosphere(k), & ! Input
140  rtv ) ! Internal variable
141  ! ----------------------------------------------------------- !
142  ! Adding method to add the layer to the present level !
143  ! to compute upward radiances and reflection matrix !
144  ! at new level. !
145  ! ----------------------------------------------------------- !
146 
147  temporal_matrix = -matmul(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k), &
148  rtv%s_Layer_Refl(1:rtv%n_Angles,1:rtv%n_Angles,k))
149  DO i = 1, rtv%n_Angles
150  temporal_matrix(i,i) = one + temporal_matrix(i,i)
151  END DO
152 
153  rtv%Inv_Gamma(1:rtv%n_Angles,1:rtv%n_Angles,k) = matinv(temporal_matrix, error_status)
154  IF( error_status /= success ) THEN
155  WRITE( message,'("Error in matrix inversion matinv(temporal_matrix, Error_Status) ")' )
156  CALL display_message( routine_name, &
157  trim(message), &
158  error_status )
159  RETURN
160  END IF
161 
162  rtv%Inv_GammaT(1:rtv%n_Angles,1:rtv%n_Angles,k) = &
163  matmul(rtv%s_Layer_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k), rtv%Inv_Gamma(1:rtv%n_Angles,1:rtv%n_Angles,k))
164  refl_down(1:rtv%n_Angles,k) = matmul(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k), &
165  rtv%s_Layer_Source_DOWN(1:rtv%n_Angles,k))
166 
167  rtv%s_Level_Rad_UP(1:rtv%n_Angles,k-1 )=rtv%s_Layer_Source_UP(1:rtv%n_Angles,k)+ &
168  matmul(rtv%Inv_GammaT(1:rtv%n_Angles,1:rtv%n_Angles,k),refl_down(1:rtv%n_Angles,k) &
169  +rtv%s_Level_Rad_UP(1:rtv%n_Angles,k ))
170  rtv%Refl_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k) = matmul(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k), &
171  rtv%s_Layer_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k))
172  rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k-1)=rtv%s_Layer_Refl(1:rtv%n_Angles,1:rtv%n_Angles,k) + &
173  matmul(rtv%Inv_GammaT(1:rtv%n_Angles,1:rtv%n_Angles,k),rtv%Refl_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k))
174 
175  ELSE
176  DO i = 1, rtv%n_Angles
177  rtv%s_Layer_Trans(i,i,k) = exp(-t_od(k)/rtv%COS_Angle(i))
178  rtv%s_Layer_Source_UP(i,k) = rtv%Planck_Atmosphere(k) * (one - rtv%s_Layer_Trans(i,i,k) )
179  rtv%s_Layer_Source_DOWN(i,k) = rtv%s_Layer_Source_UP(i,k)
180 
181  END DO
182 
183  ! Adding method
184  DO i = 1, rtv%n_Angles
185  rtv%s_Level_Rad_UP(i,k-1 )=rtv%s_Layer_Source_UP(i,k)+ &
186  rtv%s_Layer_Trans(i,i,k)*(sum(rtv%s_Level_Refl_UP(i,1:rtv%n_Angles,k)*rtv%s_Layer_Source_DOWN(1:rtv%n_Angles,k)) &
187  +rtv%s_Level_Rad_UP(i,k ))
188  ENDDO
189  DO i = 1, rtv%n_Angles
190  DO j = 1, rtv%n_Angles
191  rtv%s_Level_Refl_UP(i,j,k-1)=rtv%s_Layer_Trans(i,i,k)*rtv%s_Level_Refl_UP(i,j,k)*rtv%s_Layer_Trans(j,j,k)
192  ENDDO
193  ENDDO
194  ENDIF
195 
196  10 CONTINUE
197 
198  ! Adding reflected cosmic background radiation
199  IF( rtv%mth_Azi == 0 ) THEN
200  DO i = 1, rtv%n_Angles
201  rtv%s_Level_Rad_UP(i,0)=rtv%s_Level_Rad_UP(i,0)+sum(rtv%s_Level_Refl_UP(i,1:rtv%n_Angles,0))*cosmic_background
202  ENDDO
203  END IF
204 
205  RETURN
206 
207  END SUBROUTINE crtm_ada
208 
209  SUBROUTINE crtm_amom_layer( n_streams, & ! Input, number of streams
210  nZ, & ! Input, number of angles
211  KL, & ! Input, KL-th layer
212  single_albedo, & ! Input, single scattering albedo
213  optical_depth, & ! Input, layer optical depth
214  total_opt, & ! Input, accumulated optical depth from the top to current layer top
215  COS_Angle, & ! Input, COSINE of ANGLES
216  COS_Weight, & ! Input, GAUSSIAN Weights
217  ff, & ! Input, Phase matrix (forward part)
218  bb, & ! Input, Phase matrix (backward part)
219  planck_func, & ! Input, Planck for layer temperature
220  rtv) ! Output, layer transmittance, reflectance, and source
221 ! ---------------------------------------------------------------------------------------
222 ! FUNCTION
223 ! Compute layer transmission, reflection matrices and source function
224 ! at the top and bottom of the layer.
225 !
226 ! Method and References
227 ! The transmittance and reflectance matrices is further derived from
228 ! matrix operator method. The matrix operator method is referred to the paper by
229 !
230 ! Weng, F., and Q. Liu, 2003: Satellite Data Assimilation in Numerical Weather Prediction
231 ! Model: Part 1: Forward Radiative Transfer and Jacobian Modeling in Cloudy Atmospheres,
232 ! J. Atmos. Sci., 60, 2633-2646.
233 !
234 ! see also ADA method.
235 ! Quanhua Liu
236 ! Quanhua.Liu@noaa.gov
237 ! ----------------------------------------------------------------------------------------
238  IMPLICIT NONE
239  INTEGER, INTENT(IN) :: n_streams,nZ,KL
240  TYPE(RTV_type), INTENT( INOUT ) :: RTV
241  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
242  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
243  REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt
244 
245  ! internal variables
246  REAL(fp), DIMENSION(nZ,nZ) :: trans, refl, tempo
247  REAL(fp) :: s, c, xx
248  INTEGER :: i,j,N2,N2_1
249  INTEGER :: Error_Status
250  REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ)
251  REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ)
252  CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer'
253  CHARACTER(256) :: Message
254 
255  ! for small layer optical depth, single scattering is applied.
256  IF( optical_depth < delta_optical_depth ) THEN
257  s = optical_depth * single_albedo
258  DO i = 1, nz
259  rtv%Thermal_C(i,kl) = zero
260  c = s/cos_angle(i)
261  DO j = 1, nz
262  rtv%s_Layer_Refl(i,j,kl) = c * bb(i,j) * cos_weight(j)
263  rtv%s_Layer_Trans(i,j,kl) = c * ff(i,j) * cos_weight(j)
264  IF( i == j ) THEN
265  rtv%s_Layer_Trans(i,i,kl) = rtv%s_Layer_Trans(i,i,kl) + &
266  one - optical_depth/cos_angle(i)
267  END IF
268  IF( rtv%mth_Azi == 0 ) THEN
269  rtv%Thermal_C(i,kl) = rtv%Thermal_C(i,kl) + &
270  ( rtv%s_Layer_Refl(i,j,kl) + rtv%s_Layer_Trans(i,j,kl) )
271  END IF
272  ENDDO
273 
274  IF( rtv%mth_Azi == 0 ) THEN
275  rtv%s_Layer_Source_UP(i,kl) = ( one - rtv%Thermal_C(i,kl) ) * planck_func
276  rtv%s_Layer_Source_DOWN(i,kl) = rtv%s_Layer_Source_UP(i,kl)
277  END IF
278  ENDDO
279 
280  RETURN
281 
282  END IF
283  !
284  ! for numerical stability,
285  IF( single_albedo < max_albedo ) THEN
286  s = single_albedo
287  ELSE
288  s = max_albedo
289  END IF
290  !
291  ! building phase matrices
292  DO i = 1, nz
293  c = s/cos_angle(i)
294  DO j = 1, nz
295  rtv%PM(i,j,kl) = c * bb(i,j) * cos_weight(j)
296  rtv%PP(i,j,kl) = c * ff(i,j) * cos_weight(j)
297  ENDDO
298  rtv%PP(i,i,kl) = rtv%PP(i,i,kl) - one/cos_angle(i)
299  ENDDO
300  rtv%PPM(1:nz,1:nz,kl) = rtv%PP(1:nz,1:nz,kl) - rtv%PM(1:nz,1:nz,kl)
301  rtv%i_PPM(1:nz,1:nz,kl) = matinv( rtv%PPM(1:nz,1:nz,kl), error_status )
302  IF( error_status /= success ) THEN
303  WRITE( message,'("Error in matrix inversion matinv( RTV%PPM(1:nZ,1:nZ,KL), Error_Status ) ")' )
304  CALL display_message( routine_name, &
305  trim(message), &
306  error_status )
307  RETURN
308  END IF
309 
310  rtv%PPP(1:nz,1:nz,kl) = rtv%PP(1:nz,1:nz,kl) + rtv%PM(1:nz,1:nz,kl)
311  rtv%HH(1:nz,1:nz,kl) = matmul( rtv%PPM(1:nz,1:nz,kl), rtv%PPP(1:nz,1:nz,kl) )
312  !
313  ! save phase element RTV%HH, call ASYMTX for calculating eigenvalue and vectors.
314  tempo = rtv%HH(1:nz,1:nz,kl)
315  CALL asymtx(tempo,nz,nz,nz,rtv%EigVe(1:nz,1:nz,kl),rtv%EigVa(1:nz,kl),error_status)
316  DO i = 1, nz
317  IF( rtv%EigVa(i,kl) > zero ) THEN
318  rtv%EigValue(i,kl) = sqrt( rtv%EigVa(i,kl) )
319  ELSE
320  rtv%EigValue(i,kl) = zero
321  END IF
322  END DO
323 
324  DO i = 1, nz
325  DO j = 1, nz
326  rtv%EigVeVa(i,j,kl) = rtv%EigVe(i,j,kl) * rtv%EigValue(j,kl)
327  END DO
328  END DO
329  rtv%EigVeF(1:nz,1:nz,kl) = matmul( rtv%i_PPM(1:nz,1:nz,kl), rtv%EigVeVa(1:nz,1:nz,kl) )
330 
331  ! compute layer reflection, transmission and source function
332  rtv%Gp(1:nz,1:nz,kl) = ( rtv%EigVe(1:nz,1:nz,kl) + rtv%EigVeF(1:nz,1:nz,kl) )/2.0_fp
333  rtv%Gm(1:nz,1:nz,kl) = ( rtv%EigVe(1:nz,1:nz,kl) - rtv%EigVeF(1:nz,1:nz,kl) )/2.0_fp
334  rtv%i_Gm(1:nz,1:nz,kl) = matinv( rtv%Gm(1:nz,1:nz,kl), error_status)
335 
336  IF( error_status /= success ) THEN
337  WRITE( message,'("Error in matrix inversion matinv( RTV%Gm(1:nZ,1:nZ,KL), Error_Status) ")' )
338  CALL display_message( routine_name, &
339  trim(message), &
340  error_status )
341  RETURN
342  END IF
343 
344  DO i = 1, nz
345  xx = rtv%EigValue(i,kl)*optical_depth
346  rtv%Exp_x(i,kl) = exp(-xx)
347  END DO
348 
349  DO i = 1, nz
350  DO j = 1, nz
351  rtv%A1(i,j,kl) = rtv%Gp(i,j,kl) * rtv%Exp_x(j,kl)
352  rtv%A4(i,j,kl) = rtv%Gm(i,j,kl) * rtv%Exp_x(j,kl)
353  END DO
354  END DO
355 
356  rtv%A2(1:nz,1:nz,kl) = matmul( rtv%i_Gm(1:nz,1:nz,kl), rtv%A1(1:nz,1:nz,kl) )
357  rtv%A3(1:nz,1:nz,kl) = matmul( rtv%Gp(1:nz,1:nz,kl), rtv%A2(1:nz,1:nz,kl) )
358  rtv%A5(1:nz,1:nz,kl) = matmul( rtv%A1(1:nz,1:nz,kl), rtv%A2(1:nz,1:nz,kl) )
359  rtv%A6(1:nz,1:nz,kl) = matmul( rtv%A4(1:nz,1:nz,kl), rtv%A2(1:nz,1:nz,kl) )
360  rtv%Gm_A5(1:nz,1:nz,kl) = rtv%Gm(1:nz,1:nz,kl) - rtv%A5(1:nz,1:nz,kl)
361  rtv%i_Gm_A5(1:nz,1:nz,kl) = matinv(rtv%Gm_A5(1:nz,1:nz,kl), error_status)
362  IF( error_status /= success ) THEN
363  WRITE( message,'("Error in matrix inversion matinv(RTV%Gm_A5(1:nZ,1:nZ,KL), Error_Status) ")' )
364  CALL display_message( routine_name, &
365  trim(message), &
366  error_status )
367  RETURN
368  END IF
369  trans = matmul( rtv%A4(1:nz,1:nz,kl) - rtv%A3(1:nz,1:nz,kl), rtv%i_Gm_A5(1:nz,1:nz,kl) )
370  refl = matmul( rtv%Gp(1:nz,1:nz,kl) - rtv%A6(1:nz,1:nz,kl), rtv%i_Gm_A5(1:nz,1:nz,kl) )
371 
372  ! post processing
373  rtv%s_Layer_Trans(1:nz,1:nz,kl) = trans(:,:)
374  rtv%s_Layer_Refl(1:nz,1:nz,kl) = refl(:,:)
375  rtv%s_Layer_Source_UP(:,kl) = zero
376  IF( rtv%mth_Azi == 0 ) THEN
377  DO i = 1, nz
378  rtv%Thermal_C(i,kl) = zero
379  DO j = 1, n_streams
380  rtv%Thermal_C(i,kl) = rtv%Thermal_C(i,kl) + (trans(i,j) + refl(i,j) )
381  END DO
382  IF ( i == nz .AND. nz == (n_streams+1) ) THEN
383  rtv%Thermal_C(i,kl) = rtv%Thermal_C(i,kl) + trans(nz,nz)
384  END IF
385  rtv%s_Layer_Source_UP(i,kl) = ( one - rtv%Thermal_C(i,kl) ) * planck_func
386  rtv%s_Layer_Source_DOWN(i,kl) = rtv%s_Layer_Source_UP(i,kl)
387  END DO
388  END IF
389 
390  ! compute visible part for visible channels during daytime
391  IF( rtv%Visible_Flag_true ) THEN
392  n2 = 2 * nz
393  n2_1 = n2 - 1
394  source_up = zero
395  source_down = zero
396  !
397  ! Solar source
398  sfactor = single_albedo*rtv%Solar_irradiance/pi
399  IF( rtv%mth_Azi == 0 ) sfactor = sfactor/two
400  expfactor = exp(-optical_depth/rtv%COS_SUN)
401  s_transmittance = exp(-total_opt/rtv%COS_SUN)
402 
403  DO i = 1, nz
404  solar(i) = -bb(i,nz+1)*sfactor
405  solar(i+nz) = -ff(i,nz+1)*sfactor
406 
407  DO j = 1, nz
408  v0(i,j) = single_albedo * ff(i,j) * cos_weight(j)
409  v0(i+nz,j) = single_albedo * bb(i,j) * cos_weight(j)
410  v0(i,j+nz) = v0(i+nz,j)
411  v0(nz+i,j+nz) = v0(i,j)
412  ENDDO
413  v0(i,i) = v0(i,i) - one - cos_angle(i)/rtv%COS_SUN
414  v0(i+nz,i+nz) = v0(i+nz,i+nz) - one + cos_angle(i)/rtv%COS_SUN
415  ENDDO
416 
417  v1(1:n2_1,1:n2_1) = matinv(v0(1:n2_1,1:n2_1), error_status)
418  IF( error_status /= success ) THEN
419  WRITE( message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' )
420  CALL display_message( routine_name, &
421  trim(message), &
422  error_status )
423  RETURN
424  END IF
425 
426  solar1(1:n2_1) = matmul( v1(1:n2_1,1:n2_1), solar(1:n2_1) )
427  solar1(n2) = zero
428  sfac2 = solar(n2) - sum( v0(n2,1:n2_1)*solar1(1:n2_1) )
429 
430 
431  DO i = 1, nz
432  source_up(i) = solar1(i)
433  source_down(i) = expfactor*solar1(i+nz)
434  DO j = 1, nz
435  source_up(i) =source_up(i)-refl(i,j)*solar1(j+nz)-trans(i,j)*expfactor*solar1(j)
436  source_down(i) =source_down(i) -trans(i,j)*solar1(j+nz) -refl(i,j)*expfactor*solar1(j)
437  END DO
438  END DO
439  ! specific treatment for downeward source function
440  IF( abs( v0(n2,n2) ) > 0.0001_fp ) THEN
441  source_down(nz) =source_down(nz) +(expfactor-trans(nz,nz))*sfac2/v0(n2,n2)
442  ELSE
443  source_down(nz) =source_down(nz) -expfactor*sfac2*optical_depth/cos_angle(nz)
444  END IF
445 
446  source_up(1:nz) = source_up(1:nz)*s_transmittance
447  source_down(1:nz) = source_down(1:nz)*s_transmittance
448 
449  rtv%s_Layer_Source_UP(1:nz,kl) = rtv%s_Layer_Source_UP(1:nz,kl)+source_up(1:nz)
450  rtv%s_Layer_Source_DOWN(1:nz,kl) = rtv%s_Layer_Source_DOWN(1:nz,kl)+source_down(1:nz)
451  END IF
452 
453  RETURN
454 
455  END SUBROUTINE crtm_amom_layer
456 
457  SUBROUTINE crtm_ada_tl(n_Layers, & ! Input number of atmospheric layers
458  w, & ! Input layer scattering albedo
459  T_OD, & ! Input layer optical depth
460  cosmic_background, & ! Input cosmic background radiance
461  emissivity, & ! Input surface emissivity
462  direct_reflectivity, & ! Input direct reflectivity
463  RTV, & ! Input structure containing forward part results
464  Planck_Atmosphere_TL, & ! Input tangent-linear atmospheric layer Planck radiance
465  Planck_Surface_TL, & ! Input TL surface Planck radiance
466  w_TL, & ! Input TL layer scattering albedo
467  T_OD_TL, & ! Input TL layer optical depth
468  emissivity_TL, & ! Input TL surface emissivity
469  reflectivity_TL, & ! Input TL reflectivity
470  direct_reflectivity_TL, & ! Input TL direct reflectivity
471  Pff_TL, & ! Input TL forward phase matrix
472  Pbb_TL, & ! Input TL backward phase matrix
473  s_rad_up_TL) ! Output TL upward radiance
474 ! ------------------------------------------------------------------------- !
475 ! FUNCTION: !
476 ! This subroutine calculates IR/MW tangent-linear radiance at the top of !
477 ! the atmosphere including atmospheric scattering. The structure RTV !
478 ! carried in forward part results. !
479 ! The CRTM_ADA_TL algorithm computes layer tangent-linear reflectance and !
480 ! transmittance as well as source function by the subroutine !
481 ! CRTM_Doubling_layer as source function by the subroutine !
482 ! CRTM_Doubling_layer, then uses !
483 ! an adding method to integrate the layer and surface components. !
484 ! !
485 ! Quanhua Liu Quanhua.Liu@noaa.gov !
486 ! ------------------------------------------------------------------------- !
487  IMPLICIT NONE
488  INTEGER, INTENT(IN) :: n_layers
489  TYPE(rtv_type), INTENT(IN) :: rtv
490  REAL (fp), INTENT(IN), DIMENSION( : ) :: w,t_od
491  REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity,direct_reflectivity
492  REAL (fp), INTENT(IN) :: cosmic_background
493 
494  REAL (fp),INTENT(IN),DIMENSION( :,:,: ) :: pff_tl, pbb_tl
495  REAL (fp),INTENT(IN),DIMENSION( : ) :: w_tl,t_od_tl
496  REAL (fp),INTENT(IN),DIMENSION( 0: ) :: planck_atmosphere_tl
497  REAL (fp),INTENT(IN) :: planck_surface_tl
498  REAL (fp),INTENT(IN),DIMENSION( : ) :: emissivity_tl
499  REAL (fp),INTENT(IN),DIMENSION( :,: ) :: reflectivity_tl
500  REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_tl
501  REAL (fp),INTENT(INOUT),DIMENSION( : ) :: direct_reflectivity_tl
502  ! -------------- internal variables --------------------------------- !
503  ! Abbreviations: !
504  ! s: scattering, rad: radiance, trans: transmission, !
505  ! refl: reflection, up: upward, down: downward !
506  ! --------------------------------------------------------------------!
507  REAL (fp), DIMENSION(RTV%n_Angles,RTV%n_Angles) :: temporal_matrix_tl
508 
509  REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down
510  REAL (fp), DIMENSION( RTV%n_Angles ) :: s_source_up_tl,s_source_down_tl,refl_down_tl
511 
512  REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles ) :: s_trans_tl,s_refl_tl,refl_trans_tl
513  REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles ) :: s_refl_up_tl,inv_gamma_tl,inv_gammat_tl
514  REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_tl
515  INTEGER :: i, j, k, nz
516 !
517  nz = rtv%n_Angles
518 
519  total_opt(0) = zero
520  total_opt_tl(0) = zero
521  DO k = 1, n_layers
522  total_opt(k) = total_opt(k-1) + t_od(k)
523  total_opt_tl(k) = total_opt_tl(k-1) + t_od_tl(k)
524  END DO
525 
526  refl_trans_tl = zero
527  s_rad_up_tl = zero
528  s_refl_up_tl = reflectivity_tl
529  IF( rtv%mth_Azi == 0 ) THEN
530  s_rad_up_tl = emissivity_tl * rtv%Planck_Surface + emissivity * planck_surface_tl
531  END IF
532 
533  IF( rtv%Solar_Flag_true ) THEN
534  s_rad_up_tl = s_rad_up_tl+direct_reflectivity_tl*rtv%COS_SUN*rtv%Solar_irradiance/pi &
535  * exp(-total_opt(n_layers)/rtv%COS_SUN) &
536  - direct_reflectivity * rtv%Solar_irradiance/pi &
537  * total_opt_tl(n_layers) * exp(-total_opt(n_layers)/rtv%COS_SUN)
538  END IF
539 
540  DO 10 k = n_layers, 1, -1
541  s_source_up_tl = zero
542  s_source_down_tl = zero
543  s_trans_tl = zero
544  s_refl_tl = zero
545  inv_gammat_tl = zero
546  inv_gamma_tl = zero
547  refl_down_tl = zero
548 !
549 ! Compute tranmission and reflection matrices for a layer
550  IF(w(k) > scattering_albedo_threshold) THEN
551  ! ----------------------------------------------------------- !
552  ! CALL Doubling algorithm to computing forward and tagent !
553  ! layer transmission, reflection, and source functions. !
554  ! ----------------------------------------------------------- !
555 
556  call crtm_amom_layer_tl(rtv%n_Streams,rtv%n_Angles,k,w(k),t_od(k),total_opt(k-1), & !Input
557  rtv%COS_Angle(1:rtv%n_Angles),rtv%COS_Weight(1:rtv%n_Angles) , & !Input
558  rtv%Pff(:,:,k), rtv%Pbb(:,:,k),rtv%Planck_Atmosphere(k) , & !Input
559  w_tl(k),t_od_tl(k),total_opt_tl(k-1),pff_tl(:,:,k) , & !Input
560  pbb_tl(:,:,k),planck_atmosphere_tl(k),rtv , & !Input
561  s_trans_tl,s_refl_tl,s_source_up_tl,s_source_down_tl) !Output
562 
563 ! Adding method
564  temporal_matrix_tl = -matmul(s_refl_up_tl,rtv%s_Layer_Refl(1:rtv%n_Angles,1:rtv%n_Angles,k)) &
565  - matmul(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k),s_refl_tl)
566 
567  temporal_matrix_tl = matmul(rtv%Inv_Gamma(1:rtv%n_Angles,1:rtv%n_Angles,k),temporal_matrix_tl)
568  inv_gamma_tl = -matmul(temporal_matrix_tl,rtv%Inv_Gamma(1:rtv%n_Angles,1:rtv%n_Angles,k))
569 
570  inv_gammat_tl = matmul(s_trans_tl, rtv%Inv_Gamma(1:rtv%n_Angles,1:rtv%n_Angles,k)) &
571  + matmul(rtv%s_Layer_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k), inv_gamma_tl)
572 
573  refl_down(1:rtv%n_Angles,k) = matmul(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k), &
574  rtv%s_Layer_Source_DOWN(1:rtv%n_Angles,k))
575  refl_down_tl(:) = matmul(s_refl_up_tl,rtv%s_Layer_Source_DOWN(1:rtv%n_Angles,k)) &
576  + matmul(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k),s_source_down_tl(:))
577  s_rad_up_tl(1:rtv%n_Angles)=s_source_up_tl(1:rtv%n_Angles)+ &
578  matmul(inv_gammat_tl,refl_down(:,k)+rtv%s_Level_Rad_UP(1:rtv%n_Angles,k)) &
579  +matmul(rtv%Inv_GammaT(1:rtv%n_Angles,1:rtv%n_Angles,k),refl_down_tl(1:rtv%n_Angles)+s_rad_up_tl(1:rtv%n_Angles))
580 
581  refl_trans_tl = matmul(s_refl_up_tl,rtv%s_Layer_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k)) &
582  + matmul(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k),s_trans_tl)
583 
584  s_refl_up_tl=s_refl_tl+matmul(inv_gammat_tl,rtv%Refl_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k)) &
585  +matmul(rtv%Inv_GammaT(1:rtv%n_Angles,1:rtv%n_Angles,k),refl_trans_tl)
586 
587  refl_trans_tl = zero
588 
589  ELSE
590 
591  DO i = 1, rtv%n_Angles
592  s_trans_tl(i,i) = -t_od_tl(k)/rtv%COS_Angle(i) * rtv%s_Layer_Trans(i,i,k)
593  s_source_up_tl(i) = planck_atmosphere_tl(k) * (one - rtv%s_Layer_Trans(i,i,k) ) &
594  - rtv%Planck_Atmosphere(k) * s_trans_tl(i,i)
595  s_source_down_tl(i) = s_source_up_tl(i)
596  ENDDO
597 
598 ! Adding method
599  DO i = 1, rtv%n_Angles
600  s_rad_up_tl(i)=s_source_up_tl(i) &
601  +s_trans_tl(i,i)*(sum(rtv%s_Level_Refl_UP(i,1:rtv%n_Angles,k) &
602  *rtv%s_Layer_Source_DOWN(1:rtv%n_Angles,k))+rtv%s_Level_Rad_UP(i,k)) &
603  +rtv%s_Layer_Trans(i,i,k) &
604  *(sum(s_refl_up_tl(i,1:rtv%n_Angles)*rtv%s_Layer_Source_DOWN(1:rtv%n_Angles,k) &
605  +rtv%s_Level_Refl_UP(i,1:rtv%n_Angles,k)*s_source_down_tl(1:rtv%n_Angles))+s_rad_up_tl(i))
606 
607  ENDDO
608 
609  DO i = 1, rtv%n_Angles
610  DO j = 1, rtv%n_Angles
611  s_refl_up_tl(i,j)=s_trans_tl(i,i)*rtv%s_Level_Refl_UP(i,j,k) &
612  *rtv%s_Layer_Trans(j,j,k) &
613  +rtv%s_Layer_Trans(i,i,k)*s_refl_up_tl(i,j)*rtv%s_Layer_Trans(j,j,k) &
614  +rtv%s_Layer_Trans(i,i,k)*rtv%s_Level_Refl_UP(i,j,k)*s_trans_tl(j,j)
615  ENDDO
616  ENDDO
617 
618  ENDIF
619  10 CONTINUE
620 !
621 ! Adding reflected cosmic background radiation
622  IF( rtv%mth_Azi == 0 ) THEN
623  DO i = 1, rtv%n_Angles
624  s_rad_up_tl(i)=s_rad_up_tl(i)+sum(s_refl_up_tl(i,:))*cosmic_background
625  ENDDO
626  END IF
627 !
628  RETURN
629  END SUBROUTINE crtm_ada_tl
630 
631  SUBROUTINE crtm_amom_layer_tl( n_streams, & ! Input, number of streams
632  nZ, & ! Input, number of angles
633  KL, & ! Input, KL-th layer
634  single_albedo, & ! Input, single scattering albedo
635  optical_depth, & ! Input, layer optical depth
636  total_opt, & ! Input, accumulated optical depth from the top to current layer top
637  COS_Angle, & ! Input, COSINE of ANGLES
638  COS_Weight, & ! Input, GAUSSIAN Weights
639  ff, & ! Input, Phase matrix (forward part)
640  bb, & ! Input, Phase matrix (backward part)
641  planck_func, & ! Input, Planck for layer temperature
642  single_albedo_tl, & ! Input, tangent-linear single albedo
643  optical_depth_tl, & ! Input, TL layer optical depth
644  total_opt_tl, & ! Input, accumulated TL optical depth from the top to current layer top
645  ff_tl, & ! Input, TL forward Phase matrix
646  bb_tl, & ! Input, TL backward Phase matrix
647  planck_func_tl, & ! Input, TL Planck for layer temperature
648  rtv, & ! Input, structure containing forward results
649  trans_tl, & ! Output, layer tangent-linear trans
650  refl_tl, & ! Output, layer tangent-linear refl
651  source_up_tl, & ! Output, layer tangent-linear source_up
652  source_down_tl) ! Output, layer tangent-linear source_down
653 
654 ! ---------------------------------------------------------------------------------------
655 ! FUNCTION
656 ! Compute TL layer transmission, reflection matrices and source function
657 ! at the top and bottom of the layer.
658 !
659 ! see also ADA method.
660 ! Quanhua Liu
661 ! Quanhua.Liu@noaa.gov
662 ! ----------------------------------------------------------------------------------------
663  IMPLICIT NONE
664  INTEGER, INTENT(IN) :: n_streams,nZ,KL
665  TYPE(RTV_type), INTENT( IN ) :: RTV
666  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
667  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
668  REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt
669  !
670  ! internal variables
671  REAL(fp) :: s, c, s_TL,c_TL,xx_TL
672  INTEGER :: i,j
673  INTEGER :: Error_Status
674  !
675  ! Tangent-Linear Part
676  REAL(fp), INTENT(OUT), DIMENSION( :,: ) :: trans_TL,refl_TL
677  REAL(fp), INTENT(OUT), DIMENSION( : ) :: source_up_TL,source_down_TL
678 
679  REAL(fp), INTENT(IN) :: single_albedo_TL
680  REAL(fp), INTENT(IN) :: optical_depth_TL,Planck_Func_TL,total_opt_TL
681  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff_TL,bb_TL
682 
683  REAL(fp), DIMENSION(nZ) :: Exp_x_TL,EigVa_TL,EigValue_TL
684  REAL(fp), DIMENSION(nZ,nZ) :: i_Gm_TL, Gm_TL, Gp_TL, EigVe_TL, EigVeF_TL, EigVeVa_TL
685  REAL(fp), DIMENSION(nZ,nZ) :: A1_TL,A2_TL,A3_TL,A4_TL,A5_TL,A6_TL,Gm_A5_TL,i_Gm_A5_TL
686  REAL(fp), DIMENSION(nZ,nZ) :: HH_TL,PM_TL,PP_TL,PPM_TL,i_PPM_TL,PPP_TL
687 
688  REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ)
689  REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ)
690  REAL(fp) :: EXPfactor_TL,Sfactor_TL,s_transmittance_TL,Solar_TL(2*nZ),V0_TL(2*nZ,2*nZ),Solar1_TL(2*nZ)
691  REAL(fp) :: Sfac2_TL
692  REAL(fp), DIMENSION( nZ ) :: thermal_up_TL,thermal_down_TL
693  REAL(fp), DIMENSION(nZ,nZ) :: trans, refl
694  INTEGER :: N2, N2_1
695  REAL(fp) :: Thermal_C_TL
696  CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer_TL'
697  CHARACTER(256) :: Message
698  !
699  ! for small layer optical depth, single scattering is applied.
700  IF( optical_depth < delta_optical_depth ) THEN
701  s = optical_depth * single_albedo
702  s_tl = optical_depth_tl * single_albedo + optical_depth * single_albedo_tl
703  DO i = 1, nz
704  thermal_c_tl = zero
705  c = s/cos_angle(i)
706  c_tl = s_tl/cos_angle(i)
707  DO j = 1, nz
708  refl_tl(i,j) = (c_tl * bb(i,j) + c*bb_tl(i,j) )* cos_weight(j)
709  trans_tl(i,j) = (c_tl * ff(i,j) + c*ff_tl(i,j) )* cos_weight(j)
710  IF( i == j ) THEN
711  trans_tl(i,j) = trans_tl(i,j) - optical_depth_tl/cos_angle(i)
712  END IF
713 
714  IF( rtv%mth_Azi == 0 ) THEN
715  thermal_c_tl = thermal_c_tl + refl_tl(i,j) + trans_tl(i,j)
716  END IF
717  ENDDO
718  IF( rtv%mth_Azi == 0 ) THEN
719  source_up_tl(i) = -thermal_c_tl * planck_func + &
720  ( one - rtv%Thermal_C(i,kl) ) * planck_func_tl
721  source_down_tl(i) = source_up_tl(i)
722  END IF
723  ENDDO
724  RETURN
725  END IF
726 
727  !
728  ! for numerical stability,
729  IF( single_albedo < max_albedo ) THEN
730  s = single_albedo
731  s_tl = single_albedo_tl
732  ELSE
733  s = max_albedo
734  s_tl = 0.0_fp
735  END IF
736  !
737  ! building TL phase matrices
738  DO i = 1, nz
739  c = s/cos_angle(i)
740  c_tl = s_tl/cos_angle(i)
741  DO j = 1, nz
742  pm_tl(i,j) = (c_tl * bb(i,j) + c*bb_tl(i,j) ) * cos_weight(j)
743  pp_tl(i,j) = (c_tl * ff(i,j) + c*ff_tl(i,j) ) * cos_weight(j)
744  END DO
745  ENDDO
746 
747  ppm_tl(:,:) = pp_tl(:,:) - pm_tl(:,:)
748  i_ppm_tl(:,:) = - matmul( rtv%i_PPM(1:nz,1:nz,kl), matmul(ppm_tl(:,:),rtv%i_PPM(1:nz,1:nz,kl)) )
749  ppp_tl(:,:) = pp_tl(:,:) + pm_tl(:,:)
750  hh_tl(:,:) = matmul( ppm_tl(:,:), rtv%PPP(1:nz,1:nz,kl) )+matmul( rtv%PPM(1:nz,1:nz,kl), ppp_tl(:,:) )
751  !
752  ! compute TL eigenvectors EigVe, and eigenvalues EigVa
753  CALL asymtx_tl(nz,rtv%EigVe(1:nz,1:nz,kl),rtv%EigVa(1:nz,kl),hh_tl, &
754  eigve_tl,eigva_tl,error_status)
755 
756  DO i = 1, nz
757  IF( rtv%EigVa(i,kl) > zero ) THEN
758  eigvalue_tl(i) = 0.5_fp*eigva_tl(i)/rtv%EigValue(i,kl)
759  ELSE
760  eigvalue_tl(i) = zero
761  END IF
762  END DO
763  eigveva_tl = zero
764 
765  DO i = 1, nz
766  DO j = 1, nz
767  eigveva_tl(i,j) = eigve_tl(i,j) * rtv%EigValue(j,kl)+rtv%EigVe(i,j,kl) * eigvalue_tl(j)
768  END DO
769  END DO
770  eigvef_tl(:,:) = matmul( i_ppm_tl(:,:), rtv%EigVeVa(1:nz,1:nz,kl) ) &
771  + matmul( rtv%i_PPM(1:nz,1:nz,kl), eigveva_tl(:,:) )
772  !
773  ! compute TL reflection and transmission matrices, TL source function
774  gp_tl(:,:) = ( eigve_tl(:,:) + eigvef_tl(:,:) )/2.0_fp
775  gm_tl(:,:) = ( eigve_tl(:,:) - eigvef_tl(:,:) )/2.0_fp
776  i_gm_tl = -matmul( rtv%i_Gm(1:nz,1:nz,kl), matmul(gm_tl,rtv%i_Gm(1:nz,1:nz,kl)) )
777  DO i = 1, nz
778  xx_tl = eigvalue_tl(i)*optical_depth+rtv%EigValue(i,kl)*optical_depth_tl
779  exp_x_tl(i) = -xx_tl*rtv%Exp_x(i,kl)
780  END DO
781 
782  DO i = 1, nz
783  DO j = 1, nz
784  a1_tl(i,j) = gp_tl(i,j)* rtv%Exp_x(j,kl)+ rtv%Gp(i,j,kl)* exp_x_tl(j)
785  a4_tl(i,j) = gm_tl(i,j)* rtv%Exp_x(j,kl)+ rtv%Gm(i,j,kl)* exp_x_tl(j)
786  END DO
787  END DO
788  a2_tl(:,:) = matmul(i_gm_tl(:,:),rtv%A1(1:nz,1:nz,kl))+matmul(rtv%i_Gm(1:nz,1:nz,kl),a1_tl(:,:))
789  a3_tl(:,:) = matmul(gp_tl(:,:),rtv%A2(1:nz,1:nz,kl))+matmul(rtv%Gp(1:nz,1:nz,kl),a2_tl(:,:))
790  a5_tl(:,:) = matmul(a1_tl(:,:),rtv%A2(1:nz,1:nz,kl))+matmul(rtv%A1(1:nz,1:nz,kl),a2_tl(:,:))
791  a6_tl(:,:) = matmul(a4_tl(:,:),rtv%A2(1:nz,1:nz,kl))+matmul(rtv%A4(1:nz,1:nz,kl),a2_tl(:,:))
792 
793  gm_a5_tl(:,:) = gm_tl(:,:) - a5_tl(:,:)
794  i_gm_a5_tl(:,:) = -matmul( rtv%i_Gm_A5(1:nz,1:nz,kl),matmul(gm_a5_tl,rtv%i_Gm_A5(1:nz,1:nz,kl)))
795  !
796  ! T = matmul( RTV%A4(:,:,KL) - RTV%A3(:,:,KL), RTV%i_Gm_A5(:,:,KL) )
797  trans_tl = matmul( a4_tl(:,:) - a3_tl(:,:), rtv%i_Gm_A5(1:nz,1:nz,kl) ) &
798  + matmul( rtv%A4(1:nz,1:nz,kl) - rtv%A3(1:nz,1:nz,kl), i_gm_a5_tl(:,:) )
799  refl_tl = matmul( gp_tl(:,:) - a6_tl(:,:), rtv%i_Gm_A5(1:nz,1:nz,kl) ) &
800  + matmul( rtv%Gp(1:nz,1:nz,kl) - rtv%A6(1:nz,1:nz,kl), i_gm_a5_tl(:,:) )
801 
802  trans(1:nz,1:nz) = rtv%s_Layer_Trans(1:nz,1:nz,kl)
803  refl(1:nz,1:nz) = rtv%s_Layer_Refl(1:nz,1:nz,kl)
804  source_up_tl = zero
805  source_down_tl = zero
806  !
807  ! Thermal part
808  IF( rtv%mth_Azi == 0 ) THEN
809  DO i = 1, nz
810  thermal_c_tl = zero
811  DO j = 1, n_streams
812  thermal_c_tl = thermal_c_tl + (trans_tl(i,j) + refl_tl(i,j))
813  ENDDO
814  IF(i == nz .AND. nz == (n_streams+1)) THEN
815  thermal_c_tl = thermal_c_tl + trans_tl(nz,nz)
816  ENDIF
817  thermal_up_tl(i) = -thermal_c_tl * planck_func &
818  + ( one - rtv%Thermal_C(i,kl) ) * planck_func_tl
819  thermal_down_tl(i) = thermal_up_tl(i)
820  ENDDO
821  END IF
822 
823  !
824  ! for visible channels at daytime
825  IF( rtv%Visible_Flag_true ) THEN
826  n2 = 2 * nz
827  n2_1 = n2 - 1
828  v0 = zero
829  v1 = zero
830  solar = zero
831  solar1 = zero
832  sfac2 = zero
833  v0_tl = zero
834  solar_tl = zero
835  solar1_tl = zero
836  sfac2_tl = zero
837  !
838  ! Solar source
839  sfactor = single_albedo*rtv%Solar_irradiance/pi
840  sfactor_tl = single_albedo_tl*rtv%Solar_irradiance/pi
841 
842  IF( rtv%mth_Azi == 0 ) THEN
843  sfactor = sfactor/two
844  sfactor_tl = sfactor_tl/two
845  END IF
846  expfactor = exp(-optical_depth/rtv%COS_SUN)
847  expfactor_tl = -optical_depth_tl/rtv%COS_SUN*expfactor
848 
849  s_transmittance = exp(-total_opt/rtv%COS_SUN)
850  s_transmittance_tl = -total_opt_tl/rtv%COS_SUN*s_transmittance
851 
852  DO i = 1, nz
853  solar(i) = -bb(i,nz+1)*sfactor
854  solar_tl(i) = -bb_tl(i,nz+1)*sfactor-bb(i,nz+1)*sfactor_tl
855  solar(i+nz) = -ff(i,nz+1)*sfactor
856  solar_tl(i+nz) = -ff_tl(i,nz+1)*sfactor-ff(i,nz+1)*sfactor_tl
857  DO j = 1, nz
858  v0(i,j) = single_albedo * ff(i,j) * cos_weight(j)
859  v0_tl(i,j) = single_albedo_tl*ff(i,j)*cos_weight(j)+single_albedo*ff_tl(i,j)*cos_weight(j)
860  v0(i+nz,j) = single_albedo * bb(i,j) * cos_weight(j)
861  v0_tl(i+nz,j) = single_albedo_tl*bb(i,j)*cos_weight(j)+single_albedo*bb_tl(i,j)*cos_weight(j)
862  v0(i,j+nz) = v0(i+nz,j)
863  v0_tl(i,j+nz) = v0_tl(i+nz,j)
864  v0(nz+i,j+nz) = v0(i,j)
865  v0_tl(nz+i,j+nz) = v0_tl(i,j)
866  ENDDO
867  v0(i,i) = v0(i,i) - one - cos_angle(i)/rtv%COS_SUN
868  v0(i+nz,i+nz) = v0(i+nz,i+nz) - one + cos_angle(i)/rtv%COS_SUN
869  ENDDO
870 
871  v1(1:n2_1,1:n2_1) = matinv(v0(1:n2_1,1:n2_1), error_status)
872  IF( error_status /= success ) THEN
873  WRITE( message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' )
874  CALL display_message( routine_name, &
875  trim(message), &
876  error_status )
877  RETURN
878  END IF
879 
880  solar1(1:n2_1) = matmul( v1(1:n2_1,1:n2_1), solar(1:n2_1) )
881 
882  solar(1:n2_1) = matmul( v1(1:n2_1,1:n2_1),solar(1:n2_1) )
883  solar1_tl(1:n2_1) = matmul( v0_tl(1:n2_1,1:n2_1),solar(1:n2_1) )
884  solar1_tl(1:n2_1) = -matmul( v1(1:n2_1,1:n2_1),solar1_tl(1:n2_1) ) &
885  + matmul( v1(1:n2_1,1:n2_1), solar_tl(1:n2_1) )
886 
887  solar1(n2) = zero
888  solar1_tl(n2) = zero
889  sfac2 = solar(n2) - sum( v0(n2,1:n2_1)*solar1(1:n2_1) )
890  sfac2_tl = solar_tl(n2) - sum( v0_tl(n2,1:n2_1)*solar1(1:n2_1) ) &
891  - sum( v0(n2,1:n2_1)*solar1_tl(1:n2_1) )
892 
893  DO i = 1, nz
894  source_up(i) = solar1(i)
895  source_up_tl(i) = solar1_tl(i)
896  source_down(i) = expfactor*solar1(i+nz)
897  source_down_tl(i) = expfactor_tl*solar1(i+nz)+expfactor*solar1_tl(i+nz)
898  DO j = 1, nz
899  source_up(i) =source_up(i)-refl(i,j)*solar1(j+nz)-trans(i,j)*expfactor*solar1(j)
900  source_up_tl(i) =source_up_tl(i)-refl_tl(i,j)*solar1(j+nz) -refl(i,j)*solar1_tl(j+nz) &
901  - trans_tl(i,j)*expfactor*solar1(j) - trans(i,j)*expfactor_tl*solar1(j) -trans(i,j)*expfactor*solar1_tl(j)
902  source_down(i) =source_down(i) -trans(i,j)*solar1(j+nz) -refl(i,j)*expfactor*solar1(j)
903  source_down_tl(i) =source_down_tl(i) -trans_tl(i,j)*solar1(j+nz) -trans(i,j)*solar1_tl(j+nz) &
904  -refl_tl(i,j)*expfactor*solar1(j) -refl(i,j)*expfactor_tl*solar1(j) -refl(i,j)*expfactor*solar1_tl(j)
905  END DO
906  END DO
907  !
908  ! specific treatment for downeward source function
909  IF( abs( v0(n2,n2) ) > 0.0001_fp ) THEN
910  source_down(nz) =source_down(nz) +(expfactor-trans(nz,nz))*sfac2/v0(n2,n2)
911  source_down_tl(nz) =source_down_tl(nz) +(expfactor_tl-trans_tl(nz,nz))*sfac2/v0(n2,n2) &
912  +(expfactor-trans(nz,nz))*sfac2_tl/v0(n2,n2)-(expfactor-trans(nz,nz))*sfac2*v0_tl(n2,n2)/v0(n2,n2)/v0(n2,n2)
913  ELSE
914  source_down(nz) =source_down(nz) -expfactor*sfac2*optical_depth/cos_angle(nz)
915  source_down_tl(nz) =source_down_tl(nz) -expfactor_tl*sfac2*optical_depth/cos_angle(nz) &
916  -expfactor*sfac2_tl*optical_depth/cos_angle(nz)-expfactor*sfac2*optical_depth_tl/cos_angle(nz)
917  END IF
918 
919  ! source_up(1:nZ) = source_up(1:nZ)*s_transmittance
920  source_up_tl(1:nz) = source_up_tl(1:nz)*s_transmittance+source_up(1:nz)*s_transmittance_tl
921  ! source_down(1:nZ) = source_down(1:nZ)*s_transmittance
922  source_down_tl(1:nz) = source_down_tl(1:nz)*s_transmittance + source_down(1:nz)*s_transmittance_tl
923  END IF
924 
925  source_up_tl(:) = source_up_tl(:) + thermal_up_tl(:)
926  source_down_tl(:) = source_down_tl(:) + thermal_down_tl(:)
927 
928  RETURN
929 
930  END SUBROUTINE crtm_amom_layer_tl
931 
932  SUBROUTINE crtm_ada_ad(n_Layers, & ! Input number of atmospheric layers
933  w, & ! Input layer scattering albedo
934  T_OD, & ! Input layer optical depth
935  cosmic_background, & ! Input cosmic background radiance
936  emissivity, & ! Input surface emissivity
937  direct_reflectivity, & ! surface direct reflectivity
938  RTV, & ! Input structure containing forward results
939  s_rad_up_AD, & ! Input adjoint upward radiance
940  Planck_Atmosphere_AD, & ! Output AD atmospheric layer Planck radiance
941  Planck_Surface_AD, & ! Output AD surface Planck radiance
942  w_AD, & ! Output AD layer scattering albedo
943  T_OD_AD, & ! Output AD layer optical depth
944  emissivity_AD, & ! Output AD surface emissivity
945  reflectivity_AD, & ! Output AD surface reflectivity
946  direct_reflectivity_AD, & ! Output AD surface direct reflectivity
947  Pff_AD, & ! Output AD forward phase matrix
948  Pbb_AD) ! Output AD backward phase matrix
949 ! ------------------------------------------------------------------------- !
950 ! FUNCTION: !
951 ! This subroutine calculates IR/MW adjoint radiance at the top of !
952 ! the atmosphere including atmospheric scattering. The structure RTV !
953 ! carried in forward part results. !
954 ! The CRTM_ADA_AD algorithm computes layer tangent-linear reflectance and !
955 ! transmittance as well as source function by the subroutine !
956 ! CRTM_Doubling_layer as source function by the subroutine !
957 ! CRTM_Doubling_layer, then uses !
958 ! an adding method to integrate the layer and surface components. !
959 ! !
960 ! Quanhua Liu Quanhua.Liu@noaa.gov !
961 ! ------------------------------------------------------------------------- !
962  IMPLICIT NONE
963  INTEGER, INTENT(IN) :: n_layers
964  TYPE(rtv_type), INTENT(IN) :: rtv
965  REAL (fp), INTENT(IN), DIMENSION( : ) :: w,t_od
966  REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity,direct_reflectivity
967  REAL (fp), INTENT(IN) :: cosmic_background
968 
969  REAL (fp),INTENT(INOUT),DIMENSION( :,:,: ) :: pff_ad, pbb_ad
970  REAL (fp),INTENT(INOUT),DIMENSION( : ) :: w_ad,t_od_ad
971  REAL (fp),INTENT(INOUT),DIMENSION( 0: ) :: planck_atmosphere_ad
972  REAL (fp),INTENT(INOUT) :: planck_surface_ad
973  REAL (fp),INTENT(INOUT),DIMENSION( : ) :: emissivity_ad,direct_reflectivity_ad
974  REAL (fp),INTENT(INOUT),DIMENSION( :,: ) :: reflectivity_ad
975  REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_ad
976  ! -------------- internal variables --------------------------------- !
977  ! Abbreviations: !
978  ! s: scattering, rad: radiance, trans: transmission, !
979  ! refl: reflection, up: upward, down: downward !
980  ! --------------------------------------------------------------------!
981  REAL (fp), DIMENSION(RTV%n_Angles,RTV%n_Angles) :: temporal_matrix_ad
982 
983  REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down
984  REAL (fp), DIMENSION( RTV%n_Angles ) :: s_source_up_ad,s_source_down_ad,refl_down_ad
985 
986  REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles) :: s_trans_ad,s_refl_ad,refl_trans_ad
987  REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles) :: s_refl_up_ad,inv_gamma_ad,inv_gammat_ad
988  REAL (fp) :: sum_s_ad, sums_ad, xx
989  REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_ad
990  INTEGER :: i, j, k,nz
991 !
992  nz = rtv%n_Angles
993 
994  total_opt_ad = zero
995  total_opt(0) = zero
996  DO k = 1, n_layers
997  total_opt(k) = total_opt(k-1) + t_od(k)
998  END DO
999 !
1000  s_trans_ad = zero
1001  planck_atmosphere_ad = zero
1002  planck_surface_ad = zero
1003  s_refl_up_ad = zero
1004 
1005  pff_ad = zero
1006  pbb_ad = zero
1007 ! T_OD_AD = ZERO
1008 ! Adding reflected cosmic background radiation
1009 
1010  DO i = 1, rtv%n_Angles
1011  sum_s_ad = s_rad_up_ad(i)*cosmic_background
1012  DO j = 1, rtv%n_Angles
1013  s_refl_up_ad(i,j) = sum_s_ad
1014  ENDDO
1015  ENDDO
1016 
1017 !
1018  DO 10 k = 1, n_layers
1019  s_source_up_ad = zero
1020  s_source_down_ad = zero
1021  s_trans_ad = zero
1022 !
1023 ! Compute tranmission and reflection matrices for a layer
1024  IF(w(k) > scattering_albedo_threshold) THEN
1025 
1026  refl_down(1:rtv%n_Angles,k) = matmul(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k), &
1027  rtv%s_Layer_Source_DOWN(1:rtv%n_Angles,k))
1028 
1029  s_refl_ad = s_refl_up_ad
1030  inv_gammat_ad = matmul(s_refl_up_ad,transpose(rtv%Refl_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k)))
1031  refl_trans_ad = matmul(transpose(rtv%Inv_GammaT(1:rtv%n_Angles,1:rtv%n_Angles,k)),s_refl_up_ad)
1032 
1033  s_refl_up_ad=matmul(refl_trans_ad,transpose(rtv%s_Layer_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k)))
1034  s_trans_ad=matmul(transpose(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k)),refl_trans_ad)
1035 
1036  s_source_up_ad(1:rtv%n_Angles) = s_rad_up_ad(1:rtv%n_Angles)
1037 
1038  DO i = 1, rtv%n_Angles
1039  sums_ad = s_rad_up_ad(i)
1040  DO j = 1, rtv%n_Angles
1041  inv_gammat_ad(i,j)=inv_gammat_ad(i,j)+sums_ad*(refl_down(j,k)+rtv%s_Level_Rad_UP(j,k))
1042  ENDDO
1043  ENDDO
1044 
1045  refl_down_ad(1:rtv%n_Angles)=matmul(transpose(rtv%Inv_GammaT(1:rtv%n_Angles,1:rtv%n_Angles,k)),s_rad_up_ad(1:rtv%n_Angles))
1046  s_rad_up_ad(1:rtv%n_Angles)=matmul(transpose(rtv%Inv_GammaT(1:rtv%n_Angles,1:rtv%n_Angles,k)),s_rad_up_ad(1:rtv%n_Angles))
1047 
1048  DO i = 1, rtv%n_Angles
1049  sums_ad = refl_down_ad(i)
1050  DO j = 1, rtv%n_Angles
1051  s_refl_up_ad(i,j)=s_refl_up_ad(i,j)+sums_ad*rtv%s_Layer_Source_DOWN(j,k)
1052  ENDDO
1053  ENDDO
1054 
1055  s_source_down_ad=matmul(transpose(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k)),refl_down_ad(:))
1056 
1057  s_trans_ad=s_trans_ad+matmul(inv_gammat_ad,transpose(rtv%Inv_Gamma(1:rtv%n_Angles,1:rtv%n_Angles,k)))
1058  inv_gamma_ad= matmul(transpose(rtv%s_Layer_Trans(1:rtv%n_Angles,1:rtv%n_Angles,k)),inv_gammat_ad)
1059 
1060  temporal_matrix_ad= -matmul(inv_gamma_ad,transpose(rtv%Inv_Gamma(1:rtv%n_Angles,1:rtv%n_Angles,k)))
1061  temporal_matrix_ad=matmul(transpose(rtv%Inv_Gamma(1:rtv%n_Angles,1:rtv%n_Angles,k)),temporal_matrix_ad)
1062 
1063  s_refl_up_ad=s_refl_up_ad-matmul(temporal_matrix_ad,transpose(rtv%s_Layer_Refl(1:rtv%n_Angles,1:rtv%n_Angles,k)))
1064  s_refl_ad=s_refl_ad-matmul(transpose(rtv%s_Level_Refl_UP(1:rtv%n_Angles,1:rtv%n_Angles,k)),temporal_matrix_ad)
1065  ! ----------------------------------------------------------- !
1066  ! CALL Doubling algorithm to computing forward and tagent !
1067  ! layer transmission, reflection, and source functions. !
1068  ! ----------------------------------------------------------- !
1069 
1070  call crtm_amom_layer_ad(rtv%n_Streams,rtv%n_Angles,k,w(k),t_od(k),total_opt(k-1) , & !Input
1071  rtv%COS_Angle,rtv%COS_Weight,rtv%Pff(:,:,k),rtv%Pbb(:,:,k),rtv%Planck_Atmosphere(k), & !Input
1072  s_trans_ad,s_refl_ad,s_source_up_ad,s_source_down_ad,rtv,w_ad(k),t_od_ad(k) , &
1073  total_opt_ad(k-1),pff_ad(:,:,k),pbb_ad(:,:,k),planck_atmosphere_ad(k)) !Output
1074 
1075  ELSE
1076  DO i = 1, rtv%n_Angles
1077  DO j = 1, rtv%n_Angles
1078  s_trans_ad(j,j)=s_trans_ad(j,j)+rtv%s_Layer_Trans(i,i,k)*rtv%s_Level_Refl_UP(i,j,k)*s_refl_up_ad(i,j)
1079  s_trans_ad(i,i)=s_trans_ad(i,i)+s_refl_up_ad(i,j)*rtv%s_Level_Refl_UP(i,j,k)*rtv%s_Layer_Trans(j,j,k)
1080  s_refl_up_ad(i,j)=rtv%s_Layer_Trans(i,i,k)*s_refl_up_ad(i,j)*rtv%s_Layer_Trans(j,j,k)
1081  ENDDO
1082  ENDDO
1083 ! Adding method
1084  DO i = 1, rtv%n_Angles
1085 
1086  s_source_up_ad(i)=s_rad_up_ad(i)
1087  s_trans_ad(i,i)=s_trans_ad(i,i)+s_rad_up_ad(i)*(sum(rtv%s_Level_Refl_UP(i,1:rtv%n_Angles,k) &
1088  *rtv%s_Layer_Source_DOWN(1:rtv%n_Angles,k))+rtv%s_Level_Rad_UP(i,k))
1089 
1090  sum_s_ad=rtv%s_Layer_Trans(i,i,k)*s_rad_up_ad(i)
1091  DO j = 1, rtv%n_Angles
1092  s_refl_up_ad(i,j)=s_refl_up_ad(i,j)+sum_s_ad*rtv%s_Layer_Source_DOWN(j,k)
1093  ENDDO
1094 
1095  sum_s_ad=rtv%s_Layer_Trans(i,i,k)*s_rad_up_ad(i)
1096  DO j = 1, rtv%n_Angles
1097  s_source_down_ad(j)=s_source_down_ad(j)+sum_s_ad*rtv%s_Level_Refl_UP(i,j,k)
1098  ENDDO
1099  s_rad_up_ad(i)=rtv%s_Layer_Trans(i,i,k)*s_rad_up_ad(i)
1100 
1101  ENDDO
1102 
1103  DO i = rtv%n_Angles, 1, -1
1104  s_source_up_ad(i) = s_source_up_ad(i) + s_source_down_ad(i)
1105  s_trans_ad(i,i) = s_trans_ad(i,i) - rtv%Planck_Atmosphere(k) * s_source_up_ad(i)
1106  planck_atmosphere_ad(k) = planck_atmosphere_ad(k) + s_source_up_ad(i) * (one - rtv%s_Layer_Trans(i,i,k) )
1107  s_source_up_ad(i) = zero
1108 
1109  t_od_ad(k)=t_od_ad(k)-s_trans_ad(i,i)/rtv%COS_Angle(i)*rtv%s_Layer_Trans(i,i,k)
1110  ENDDO
1111 
1112  ENDIF
1113  10 CONTINUE
1114 
1115 !
1116 
1117  IF( rtv%Solar_Flag_true ) THEN
1118  xx = rtv%Solar_irradiance/pi * exp(-total_opt(n_layers)/rtv%COS_SUN)
1119  total_opt_ad(n_layers) = total_opt_ad(n_layers) &
1120  - xx*sum(direct_reflectivity(1:rtv%n_Angles)*s_rad_up_ad(1:rtv%n_Angles))
1121  direct_reflectivity_ad = direct_reflectivity_ad + s_rad_up_ad * rtv%COS_SUN * xx
1122  END IF
1123 
1124  emissivity_ad = s_rad_up_ad * rtv%Planck_Surface
1125  planck_surface_ad = sum(emissivity(:) * s_rad_up_ad(:) )
1126 
1127 
1128  reflectivity_ad = s_refl_up_ad
1129 !
1130  s_rad_up_ad = zero
1131  s_refl_up_ad = zero
1132 
1133 
1134  DO k = n_layers, 1, -1
1135  t_od_ad(k) = t_od_ad(k) + total_opt_ad(k)
1136  total_opt_ad(k-1) = total_opt_ad(k-1) + total_opt_ad(k)
1137  END DO
1138 
1139  RETURN
1140  END SUBROUTINE crtm_ada_ad
1141 
1142 !
1143 !
1144  SUBROUTINE crtm_amom_layer_ad( n_streams, & ! Input, number of streams
1145  nZ, & ! Input, number of angles
1146  KL, & ! Input, KL-th layer
1147  single_albedo, & ! Input, single scattering albedo
1148  optical_depth, & ! Input, layer optical depth
1149  total_opt, & ! Input,
1150  COS_Angle, & ! Input, COSINE of ANGLES
1151  COS_Weight, & ! Input, GAUSSIAN Weights
1152  ff, & ! Input, Phase matrix (forward part)
1153  bb, & ! Input, Phase matrix (backward part)
1154  planck_func, & ! Input, Planck for layer temperature
1155  trans_ad, & ! Input, layer tangent-linear trans
1156  refl_ad, & ! Input, layer tangent-linear refl
1157  source_up_ad, & ! Input, layer tangent-linear source_up
1158  source_down_ad, & ! Input, layer tangent-linear source_down
1159  rtv, & ! Input, structure containing forward results
1160  single_albedo_ad, & ! Output adjoint single scattering albedo
1161  optical_depth_ad, & ! Output AD layer optical depth
1162  total_opt_ad, & ! Output AD accumulated optical depth ftom TOA to current layer top
1163  ff_ad, & ! Output AD forward Phase matrix
1164  bb_ad, & ! Output AD backward Phase matrix
1165  planck_func_ad) ! Output AD Planck for layer temperature
1166 ! ---------------------------------------------------------------------------------------
1167 ! FUNCTION
1168 ! Compute AD layer transmission, reflection matrices and source function
1169 ! at the top and bottom of the layer.
1170 !
1171 ! see also ADA method.
1172 ! Quanhua Liu
1173 ! Quanhua.Liu@noaa.gov
1174 ! ----------------------------------------------------------------------------------------
1175  IMPLICIT NONE
1176  INTEGER, INTENT(IN) :: n_streams,nZ,KL
1177  TYPE(RTV_type), INTENT( IN ) :: RTV
1178  REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb
1179  REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight
1180  REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt
1181 
1182  ! internal variables
1183  REAL(fp) :: s, c, s_AD,c_AD,xx_AD
1184  INTEGER :: i,j
1185  INTEGER :: Error_Status
1186 
1187  ! Tangent-Linear Part
1188  REAL(fp), INTENT( INOUT ), DIMENSION( :,: ) :: trans_AD,refl_AD
1189  REAL(fp), INTENT( INOUT ), DIMENSION( : ) :: source_up_AD,source_down_AD
1190  REAL(fp), INTENT( INOUT ) :: single_albedo_AD
1191  REAL(fp), INTENT( INOUT ) :: optical_depth_AD,Planck_Func_AD,total_opt_AD
1192  REAL(fp), INTENT(INOUT), DIMENSION(:,:) :: ff_AD,bb_AD
1193 
1194  REAL(fp), DIMENSION(nZ) :: Exp_x_AD,EigVa_AD,EigValue_AD
1195  REAL(fp), DIMENSION(nZ,nZ) :: i_Gm_AD, Gm_AD, Gp_AD, EigVe_AD, EigVeF_AD, EigVeVa_AD
1196  REAL(fp), DIMENSION(nZ,nZ) :: A1_AD,A2_AD,A3_AD,A4_AD,A5_AD,A6_AD,Gm_A5_AD,i_Gm_A5_AD
1197  REAL(fp), DIMENSION(nZ,nZ) :: HH_AD,PM_AD,PP_AD,PPM_AD,i_PPM_AD,PPP_AD
1198 
1199  REAL(fp), DIMENSION(nZ) :: thermal_up_AD, thermal_down_AD
1200  REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ)
1201  REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ)
1202  REAL(fp) :: EXPfactor_AD,Sfactor_AD,s_transmittance_AD,Solar_AD(2*nZ),V0_AD(2*nZ,2*nZ),Solar1_AD(2*nZ)
1203  REAL(fp) :: V1_AD(2*nZ,2*nZ),Sfac2_AD
1204  REAL(fp), DIMENSION(nZ,nZ) :: trans, refl
1205  INTEGER :: N2, N2_1
1206  REAL(fp) :: Thermal_C_AD
1207  CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer_AD'
1208  CHARACTER(256) :: Message
1209 
1210  s_ad = zero
1211  c_ad = zero
1212  !
1213  ! for small layer optical depth, single scattering is applied.
1214  IF( optical_depth < delta_optical_depth ) THEN
1215  s = optical_depth * single_albedo
1216  DO i = 1, nz
1217  c = s/cos_angle(i)
1218  IF( rtv%mth_Azi == 0 ) THEN
1219  source_up_ad(i) = source_up_ad(i) + source_down_ad(i)
1220  source_down_ad(i) = zero
1221  planck_func_ad = planck_func_ad + (one - rtv%Thermal_C(i,kl))*source_up_ad(i)
1222  thermal_c_ad = -source_up_ad(i) * planck_func
1223  END IF
1224 
1225  DO j = 1, nz
1226  IF( rtv%mth_Azi == 0 ) THEN
1227  refl_ad(i,j) = refl_ad(i,j) + thermal_c_ad
1228  trans_ad(i,j) = trans_ad(i,j) + thermal_c_ad
1229  END IF
1230  IF( i == j ) THEN
1231  optical_depth_ad = optical_depth_ad - trans_ad(i,j)/cos_angle(i)
1232  END IF
1233 
1234  c_ad = c_ad + trans_ad(i,j) * ff(i,j) * cos_weight(j)
1235  ff_ad(i,j) = ff_ad(i,j) + c * trans_ad(i,j) * cos_weight(j)
1236  c_ad = c_ad + refl_ad(i,j) * bb(i,j) * cos_weight(j)
1237  bb_ad(i,j) = bb_ad(i,j) + c * refl_ad(i,j) * cos_weight(j)
1238  ENDDO
1239 
1240  source_up_ad(i) = zero
1241  s_ad = s_ad + c_ad/cos_angle(i)
1242  c_ad = zero
1243  ENDDO
1244  optical_depth_ad = optical_depth_ad + s_ad * single_albedo
1245  single_albedo_ad = single_albedo_ad + optical_depth * s_ad
1246  RETURN
1247  END IF
1248 
1249  trans(1:nz,1:nz) = rtv%s_Layer_Trans(1:nz,1:nz,kl)
1250  refl(1:nz,1:nz) = rtv%s_Layer_Refl(1:nz,1:nz,kl)
1251 
1252  thermal_up_ad(:) = source_up_ad(:)
1253  thermal_down_ad(:) = source_down_ad(:)
1254 
1255  IF( rtv%Visible_Flag_true ) THEN
1256  n2 = 2 * nz
1257  n2_1 = n2 - 1
1258 
1259  ! forward part start ********
1260  source_up = zero
1261  source_down = zero
1262  solar_ad = zero
1263  solar1_ad = zero
1264  sfactor_ad = zero
1265  sfac2_ad = zero
1266  expfactor_ad = zero
1267  s_transmittance_ad = zero
1268  v0_ad = zero
1269  v1_ad = zero
1270  !
1271  ! Solar source
1272  sfactor = single_albedo*rtv%Solar_irradiance/pi
1273  IF( rtv%mth_Azi == 0 ) sfactor = sfactor/two
1274  expfactor = exp(-optical_depth/rtv%COS_SUN)
1275  s_transmittance = exp(-total_opt/rtv%COS_SUN)
1276 
1277  DO i = 1, nz
1278  solar(i) = -bb(i,nz+1)*sfactor
1279  solar(i+nz) = -ff(i,nz+1)*sfactor
1280 
1281  DO j = 1, nz
1282  v0(i,j) = single_albedo * ff(i,j) * cos_weight(j)
1283  v0(i+nz,j) = single_albedo * bb(i,j) * cos_weight(j)
1284  v0(i,j+nz) = v0(i+nz,j)
1285  v0(nz+i,j+nz) = v0(i,j)
1286  ENDDO
1287  v0(i,i) = v0(i,i) - one - cos_angle(i)/rtv%COS_SUN
1288  v0(i+nz,i+nz) = v0(i+nz,i+nz) - one + cos_angle(i)/rtv%COS_SUN
1289  ENDDO
1290 
1291  v1(1:n2_1,1:n2_1) = matinv(v0(1:n2_1,1:n2_1), error_status)
1292  IF( error_status /= success ) THEN
1293  WRITE( message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' )
1294  CALL display_message( routine_name, &
1295  trim(message), &
1296  error_status )
1297  RETURN
1298  END IF
1299 
1300  solar1(1:n2_1) = matmul( v1(1:n2_1,1:n2_1), solar(1:n2_1) )
1301  solar1(n2) = zero
1302  sfac2 = solar(n2) - sum( v0(n2,1:n2_1)*solar1(1:n2_1) )
1303 
1304  DO i = 1, nz
1305  source_up(i) = solar1(i)
1306  source_down(i) = expfactor*solar1(i+nz)
1307  DO j = 1, nz
1308  source_up(i) =source_up(i)-refl(i,j)*solar1(j+nz)-trans(i,j)*expfactor*solar1(j)
1309  source_down(i) =source_down(i) -trans(i,j)*solar1(j+nz) -refl(i,j)*expfactor*solar1(j)
1310  END DO
1311  END DO
1312  ! specific treatment for downeward source function
1313  IF( abs( v0(n2,n2) ) > 0.0001_fp ) THEN
1314  source_down(nz) =source_down(nz) +(expfactor-trans(nz,nz))*sfac2/v0(n2,n2)
1315  ELSE
1316  source_down(nz) =source_down(nz) -expfactor*sfac2*optical_depth/cos_angle(nz)
1317  END IF
1318 
1319  ! forward part end ********
1320  !
1321  s_transmittance_ad = s_transmittance_ad+sum(source_down(1:nz)*source_down_ad(1:nz) )
1322  source_down_ad(1:nz) = source_down_ad(1:nz)*s_transmittance
1323  s_transmittance_ad = s_transmittance_ad + sum(source_up(1:nz)*source_up_ad(1:nz) )
1324  source_up_ad(1:nz) = source_up_ad(1:nz)*s_transmittance
1325  !
1326  ! specific treatment for downeward source function
1327  IF( abs( v0(n2,n2) ) > 0.0001_fp ) THEN
1328  v0_ad(n2,n2)=v0_ad(n2,n2)-(expfactor-trans(nz,nz))*sfac2*source_down_ad(nz)/v0(n2,n2)/v0(n2,n2)
1329  sfac2_ad = sfac2_ad+(expfactor-trans(nz,nz))*source_down_ad(nz)/v0(n2,n2)
1330  expfactor_ad = expfactor_ad+source_down_ad(nz)*sfac2/v0(n2,n2)
1331  trans_ad(nz,nz) = trans_ad(nz,nz)-source_down_ad(nz)*sfac2/v0(n2,n2)
1332  ELSE
1333  optical_depth_ad = optical_depth_ad -expfactor*sfac2*source_down_ad(nz)/cos_angle(nz)
1334  sfac2_ad = sfac2_ad-expfactor*source_down_ad(nz)*optical_depth/cos_angle(nz)
1335  expfactor_ad = expfactor_ad-source_down_ad(nz)*sfac2*optical_depth/cos_angle(nz)
1336  END IF
1337 
1338  DO i = nz, 1, -1
1339  DO j = nz, 1, -1
1340  solar1_ad(j)=solar1_ad(j)-refl(i,j)*expfactor*source_down_ad(i)
1341  expfactor_ad = expfactor_ad -refl(i,j)*source_down_ad(i)*solar1(j)
1342  refl_ad(i,j) = refl_ad(i,j) -source_down_ad(i)*expfactor*solar1(j)
1343  solar1_ad(j+nz) = solar1_ad(j+nz) -trans(i,j)*source_down_ad(i)
1344  trans_ad(i,j) = trans_ad(i,j) -source_down_ad(i)*solar1(j+nz)
1345 
1346  solar1_ad(j)=solar1_ad(j)-trans(i,j)*expfactor*source_up_ad(i)
1347  expfactor_ad = expfactor_ad - trans(i,j)*source_up_ad(i)*solar1(j)
1348  trans_ad(i,j)=trans_ad(i,j) - source_up_ad(i)*expfactor*solar1(j)
1349  solar1_ad(j+nz) = solar1_ad(j+nz) -refl(i,j)*source_up_ad(i)
1350  refl_ad(i,j) = refl_ad(i,j) -source_up_ad(i)*solar1(j+nz)
1351  END DO
1352 
1353  solar1_ad(i+nz) = solar1_ad(i+nz) + expfactor * source_down_ad(i)
1354  expfactor_ad = expfactor_ad + source_down_ad(i)*solar1(i+nz)
1355  solar1_ad(i) = solar1_ad(i) + source_up_ad(i)
1356  END DO
1357 
1358  solar1_ad(1:n2_1)=solar1_ad(1:n2_1) -sfac2_ad*v0(n2,1:n2_1)
1359  v0_ad(n2,1:n2_1)=v0_ad(n2,1:n2_1) -sfac2_ad*solar1(1:n2_1)
1360  solar_ad(n2) = solar_ad(n2) + sfac2_ad
1361 
1362  solar1_ad(n2) = zero
1363  solar_ad(1:n2_1)=solar_ad(1:n2_1)+matmul( transpose(v1(1:n2_1,1:n2_1)),solar1_ad(1:n2_1) )
1364  solar(1:n2_1) = matmul( v1(1:n2_1,1:n2_1),solar(1:n2_1) )
1365  solar1_ad(1:n2_1) = -matmul( transpose(v1(1:n2_1,1:n2_1)),solar1_ad(1:n2_1) )
1366  DO i = 1, n2_1
1367  DO j = 1, n2_1
1368  v0_ad(i,j)=v0_ad(i,j)+solar1_ad(i)*solar(j)
1369  END DO
1370  END DO
1371 
1372  ! Solar source
1373  DO i = nz, 1, -1
1374  DO j = nz, 1, -1
1375  v0_ad(i,j)=v0_ad(i,j) + v0_ad(nz+i,j+nz)
1376  v0_ad(i+nz,j)=v0_ad(i+nz,j) + v0_ad(i,j+nz)
1377  bb_ad(i,j)=bb_ad(i,j) + single_albedo*v0_ad(i+nz,j)*cos_weight(j)
1378  single_albedo_ad=single_albedo_ad + v0_ad(i+nz,j)*bb(i,j)*cos_weight(j)
1379  ff_ad(i,j)=ff_ad(i,j) + single_albedo*v0_ad(i,j)*cos_weight(j)
1380  single_albedo_ad=single_albedo_ad +v0_ad(i,j)*ff(i,j)*cos_weight(j)
1381  ENDDO
1382 
1383  sfactor_ad = sfactor_ad -ff(i,nz+1)*solar_ad(i+nz)
1384  ff_ad(i,nz+1) = ff_ad(i,nz+1) -solar_ad(i+nz)*sfactor
1385  sfactor_ad = sfactor_ad -bb(i,nz+1)*solar_ad(i)
1386  bb_ad(i,nz+1)=bb_ad(i,nz+1) - solar_ad(i)*sfactor
1387  ENDDO
1388 
1389  total_opt_ad = total_opt_ad -s_transmittance_ad/rtv%COS_SUN*s_transmittance
1390  optical_depth_ad = optical_depth_ad -expfactor_ad/rtv%COS_SUN*expfactor
1391 
1392  IF( rtv%mth_Azi == 0 ) THEN
1393  sfactor_ad = sfactor_ad/two
1394  END IF
1395 
1396  single_albedo_ad = single_albedo_ad + sfactor_ad*rtv%Solar_irradiance/pi
1397 
1398  END IF
1399 
1400  ! Thermal part
1401  IF( rtv%mth_Azi == 0 ) THEN
1402  DO i = nz, 1, -1
1403  thermal_up_ad(i) = thermal_up_ad(i) + thermal_down_ad(i)
1404  thermal_down_ad(i) = zero
1405  planck_func_ad = planck_func_ad + ( one - rtv%Thermal_C(i,kl) ) * thermal_up_ad(i)
1406  thermal_c_ad = -thermal_up_ad(i) * planck_func
1407 
1408  IF ( i == nz .AND. nz == (n_streams+1) ) THEN
1409  trans_ad(nz,nz) = trans_ad(nz,nz) + thermal_c_ad
1410  END IF
1411 
1412  DO j = n_streams, 1, -1
1413  trans_ad(i,j) = trans_ad(i,j) + thermal_c_ad
1414  refl_ad(i,j) = refl_ad(i,j) + thermal_c_ad
1415  ENDDO
1416  thermal_up_ad(i) = zero
1417  ENDDO
1418  END IF
1419 !
1420  i_gm_a5_ad = matmul( transpose(rtv%Gp(1:nz,1:nz,kl)-rtv%A6(1:nz,1:nz,kl)),refl_ad )
1421  gp_ad(:,:) = matmul( refl_ad, transpose(rtv%i_Gm_A5(1:nz,1:nz,kl)) )
1422 
1423  a6_ad = - gp_ad
1424  i_gm_a5_ad = i_gm_a5_ad + matmul( transpose(rtv%A4(1:nz,1:nz,kl)-rtv%A3(1:nz,1:nz,kl)),trans_ad )
1425  a4_ad(:,:) = matmul( trans_ad, transpose(rtv%i_Gm_A5(1:nz,1:nz,kl)) )
1426  a3_ad = - a4_ad
1427  gm_a5_ad = -matmul( transpose(rtv%i_Gm_A5(1:nz,1:nz,kl)) ,matmul( i_gm_a5_ad,transpose(rtv%i_Gm_A5(1:nz,1:nz,kl)) ) )
1428  gm_ad = gm_a5_ad
1429  a5_ad = - gm_a5_ad
1430 
1431  a4_ad = a4_ad + matmul( a6_ad(:,:), transpose(rtv%A2(1:nz,1:nz,kl)) )
1432  a2_ad = matmul( transpose(rtv%A4(1:nz,1:nz,kl)),a6_ad(:,:) )
1433 
1434  a1_ad = matmul( a5_ad(:,:), transpose(rtv%A2(1:nz,1:nz,kl)) )
1435  a2_ad = a2_ad + matmul( transpose(rtv%A1(1:nz,1:nz,kl)), a5_ad(:,:) )
1436 
1437  gp_ad = gp_ad + matmul( a3_ad(:,:), transpose(rtv%A2(1:nz,1:nz,kl)) )
1438  a2_ad = a2_ad + matmul( transpose(rtv%Gp(1:nz,1:nz,kl)),a3_ad(:,:) )
1439 
1440  i_gm_ad = matmul( a2_ad(:,:), transpose(rtv%A1(1:nz,1:nz,kl)) )
1441  a1_ad = a1_ad + matmul( transpose(rtv%i_Gm(1:nz,1:nz,kl)), a2_ad(:,:) )
1442 
1443  exp_x_ad = zero
1444 
1445  DO i = nz, 1, -1
1446  DO j = nz, 1, -1
1447  gm_ad(i,j) = gm_ad(i,j) + a4_ad(i,j)* rtv%Exp_x(j,kl)
1448  exp_x_ad(j) = exp_x_ad(j) + rtv%Gm(i,j,kl)*a4_ad(i,j)
1449  gp_ad(i,j) = gp_ad(i,j) + a1_ad(i,j)* rtv%Exp_x(j,kl)
1450  exp_x_ad(j) = exp_x_ad(j) + rtv%Gp(i,j,kl)*a1_ad(i,j)
1451  END DO
1452  END DO
1453 
1454  DO i = nz, 1, -1
1455  xx_ad = -exp_x_ad(i)*rtv%Exp_x(i,kl)
1456  exp_x_ad(i) = zero
1457  eigvalue_ad(i) = xx_ad*optical_depth
1458  optical_depth_ad = optical_depth_ad + rtv%EigValue(i,kl)*xx_ad
1459  END DO
1460 
1461  gm_ad = gm_ad -matmul( transpose(rtv%i_Gm(1:nz,1:nz,kl)), matmul( i_gm_ad, transpose(rtv%i_Gm(1:nz,1:nz,kl)) ) )
1462 
1463  eigve_ad(:,:) = gm_ad(:,:)/2.0_fp
1464  eigvef_ad(:,:) = - gm_ad(:,:)/2.0_fp
1465 
1466  eigve_ad = eigve_ad + gp_ad(:,:)/2.0_fp
1467  eigvef_ad = eigvef_ad + gp_ad(:,:)/2.0_fp
1468 
1469  i_ppm_ad(:,:) = matmul( eigvef_ad(:,:), transpose(rtv%EigVeVa(1:nz,1:nz,kl)) )
1470  eigveva_ad(:,:) = matmul( transpose(rtv%i_PPM(1:nz,1:nz,kl)), eigvef_ad(:,:) )
1471 
1472  DO i = nz, 1, -1
1473  DO j = nz, 1, -1
1474  eigve_ad(i,j)=eigve_ad(i,j)+eigveva_ad(i,j)* rtv%EigValue(j,kl)
1475  eigvalue_ad(j) = eigvalue_ad(j)+rtv%EigVe(i,j,kl)*eigveva_ad(i,j)
1476  END DO
1477  END DO
1478 
1479  DO i = nz, 1, -1
1480  IF( rtv%EigVa(i,kl) > zero ) THEN
1481  eigva_ad(i) = 0.5_fp*eigvalue_ad(i)/rtv%EigValue(i,kl)
1482  ELSE
1483  eigvalue_ad(i) = zero
1484  eigva_ad(i) = zero
1485  END IF
1486  END DO
1487 
1488  ! compute eigenvectors EigVe, and eigenvalues EigVa
1489  CALL asymtx_ad(nz,rtv%EigVe(1:nz,1:nz,kl),rtv%EigVa(1:nz,kl), &
1490  eigve_ad,eigva_ad,hh_ad,error_status)
1491 
1492  ppm_ad(:,:) = matmul( hh_ad(:,:), transpose(rtv%PPP(1:nz,1:nz,kl)) )
1493  ppp_ad(:,:) = matmul( transpose(rtv%PPM(1:nz,1:nz,kl)), hh_ad(:,:) )
1494 
1495  pp_ad = ppp_ad
1496  pm_ad = ppp_ad
1497 
1498  ppm_ad(:,:) = ppm_ad(:,:)-matmul( transpose(rtv%i_PPM(1:nz,1:nz,kl)),matmul(i_ppm_ad(:,:),transpose(rtv%i_PPM(1:nz,1:nz,kl))) )
1499 
1500  pp_ad = pp_ad + ppm_ad
1501  pm_ad = pm_ad - ppm_ad
1502 
1503  IF( single_albedo < max_albedo ) THEN
1504  s = single_albedo
1505  ELSE
1506  s = max_albedo
1507  END IF
1508 
1509  c_ad = zero
1510  s_ad = zero
1511  DO i = nz, 1, -1
1512  c = s/cos_angle(i)
1513  DO j = nz, 1, -1
1514  c_ad = c_ad + pp_ad(i,j) * ff(i,j) * cos_weight(j)
1515  ff_ad(i,j) = ff_ad(i,j) + c * pp_ad(i,j) * cos_weight(j)
1516  c_ad = c_ad + pm_ad(i,j) * bb(i,j) * cos_weight(j)
1517  bb_ad(i,j) = bb_ad(i,j) + c * pm_ad(i,j) * cos_weight(j)
1518  END DO
1519  s_ad = s_ad + c_ad/cos_angle(i)
1520  c_ad = zero
1521  ENDDO
1522 !
1523  IF( single_albedo < max_albedo ) THEN
1524  s = single_albedo
1525  single_albedo_ad = s_ad + single_albedo_ad
1526  ELSE
1527  s = max_albedo
1528  s_ad = 0.0_fp
1529  END IF
1530 !
1531  RETURN
1532 
1533  END SUBROUTINE crtm_amom_layer_ad
1534 
1535 
1536 END MODULE ada_module
crtm_utility
Definition: CRTM_Utility.f90:17
crtm_parameters::zero
real(fp), parameter, public zero
Definition: CRTM_Parameters.f90:70
ada_module::crtm_ada
subroutine, public crtm_ada(n_Layers, w, T_OD, cosmic_background, emissivity, reflectivity, direct_reflectivity, RTV)
Definition: ADA_Module.f90:64
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
ada_module::crtm_amom_layer_tl
subroutine crtm_amom_layer_tl(n_streams, nZ, KL, single_albedo, optical_depth, total_opt, COS_Angle, COS_Weight, ff, bb, Planck_Func, single_albedo_TL, optical_depth_TL, total_opt_TL, ff_TL, bb_TL, Planck_Func_TL, RTV, trans_TL, refl_TL, source_up_TL, source_down_TL)
Definition: ADA_Module.f90:653
crtm_utility::asymtx_tl
subroutine, public asymtx_tl(nZ, V, VAL, A_TL, V_TL, VAL_TL, Error_Status)
Definition: CRTM_Utility.f90:822
ada_module::crtm_amom_layer
subroutine crtm_amom_layer(n_streams, nZ, KL, single_albedo, optical_depth, total_opt, COS_Angle, COS_Weight, ff, bb, Planck_Func, RTV)
Definition: ADA_Module.f90:221
message_handler
Definition: Message_Handler.f90:4
ada_module
Definition: ADA_Module.f90:14
ada_module::crtm_amom_layer_ad
subroutine crtm_amom_layer_ad(n_streams, nZ, KL, single_albedo, optical_depth, total_opt, 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, total_opt_AD, ff_AD, bb_AD, Planck_Func_AD)
Definition: ADA_Module.f90:1166
ada_module::crtm_ada_tl
subroutine, public crtm_ada_tl(n_Layers, w, T_OD, cosmic_background, emissivity, direct_reflectivity, RTV, Planck_Atmosphere_TL, Planck_Surface_TL, w_TL, T_OD_TL, emissivity_TL, reflectivity_TL, direct_reflectivity_TL, Pff_TL, Pbb_TL, s_rad_up_TL)
Definition: ADA_Module.f90:474
ada_module::module_version_id
character(*), parameter module_version_id
Definition: ADA_Module.f90:42
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
crtm_parameters::two
real(fp), parameter, public two
Definition: CRTM_Parameters.f90:72
rtv_define::rtv_type
Definition: RTV_Define.f90:106
type_kinds
Definition: Type_Kinds.f90:12
rtv_define::max_albedo
real(fp), parameter, public max_albedo
Definition: RTV_Define.f90:79
crtm_utility::asymtx
subroutine, public asymtx(AAD, M, IA, IEVEC, EVECD, EVALD, IER)
Definition: CRTM_Utility.f90:255
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
crtm_utility::asymtx_ad
subroutine, public asymtx_ad(nZ, V, VAL, V_AD, VAL_AD, A_AD, Error_Status)
Definition: CRTM_Utility.f90:883
ada_module::crtm_ada_ad
subroutine, public crtm_ada_ad(n_Layers, w, T_OD, cosmic_background, emissivity, direct_reflectivity, RTV, s_rad_up_AD, Planck_Atmosphere_AD, Planck_Surface_AD, w_AD, T_OD_AD, emissivity_AD, reflectivity_AD, direct_reflectivity_AD, Pff_AD, Pbb_AD)
Definition: ADA_Module.f90:949