! ! ADA_Module ! ! Module containing the Adding Doubling Adding (ADA) radiative ! transfer solution procedures used in the CRTM. ! ! ! CREATION HISTORY: ! Written by: Quanhua Liu, QSS at JCSDA; quanhua.liu@noaa.gov ! Yong Han, NOAA/NESDIS; yong.han@noaa.gov ! Paul van Delst; CIMMS/SSEC; paul.vandelst@noaa.gov ! 08-Jun-2004 MODULE ADA_Module ! ------------------ ! Environemnt set up ! ------------------ ! Module use statements USE RTV_Define USE CRTM_Parameters USE Type_Kinds USE Message_Handler USE CRTM_Utility IMPLICIT NONE ! -------------------- ! Default visibilities ! -------------------- ! Everything private by default PRIVATE PUBLIC CRTM_ADA PUBLIC CRTM_ADA_TL PUBLIC CRTM_ADA_AD ! ----------------- ! Module parameters ! ----------------- ! Version Id for the module CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = & '$Id: $' CONTAINS !################################################################################ !################################################################################ !## ## !## ## PUBLIC MODULE ROUTINES ## ## !## ## !################################################################################ !################################################################################ SUBROUTINE CRTM_ADA(n_Layers, & ! Input number of atmospheric layers w, & ! Input layer scattering albedo T_OD, & ! Input layer optical depth cosmic_background, & ! Input cosmic background radiance emissivity, & ! Input surface emissivity reflectivity, & ! Input surface reflectivity matrix direct_reflectivity, & ! Input surface direct reflectivity RTV) ! IN/Output upward radiance and others ! ------------------------------------------------------------------------- ! ! FUNCTION: ! ! This subroutine calculates IR/MW radiance at the top of the atmosphere ! ! including atmospheric scattering. The scheme will include solar part. ! ! The ADA algorithm computes layer reflectance and transmittance as well ! ! as source function by the subroutine CRTM_Doubling_layer, then uses ! ! an adding method to integrate the layer and surface components. ! ! ! ! Quanhua Liu Quanhua.Liu@noaa.gov ! ! ------------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, INTENT(IN) :: n_Layers INTEGER nZ TYPE(RTV_type), INTENT( INOUT ) :: RTV REAL (fp), INTENT(IN), DIMENSION( : ) :: w,T_OD REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity, direct_reflectivity REAL (fp), INTENT(IN), DIMENSION( :,: ) :: reflectivity REAL (fp), INTENT(IN) :: cosmic_background ! -------------- internal variables --------------------------------- ! ! Abbreviations: ! ! s: scattering, rad: radiance, trans: transmission, ! ! refl: reflection, up: upward, down: downward ! ! --------------------------------------------------------------------! REAL (fp), DIMENSION(RTV%n_Angles, RTV%n_Angles) :: temporal_matrix REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down REAL (fp), DIMENSION(0:n_Layers) :: total_opt INTEGER :: i, j, k, Error_Status CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_ADA' CHARACTER(256) :: Message total_opt(0) = ZERO DO k = 1, n_Layers total_opt(k) = total_opt(k-1) + T_OD(k) END DO nZ = RTV%n_Angles RTV%s_Layer_Trans = ZERO RTV%s_Layer_Refl = ZERO RTV%s_Level_Refl_UP = ZERO RTV%s_Level_Rad_UP = ZERO RTV%s_Layer_Source_UP = ZERO RTV%s_Layer_Source_DOWN = ZERO RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,n_Layers)=reflectivity(1:RTV%n_Angles,1:RTV%n_Angles) IF( RTV%mth_Azi == 0 ) THEN RTV%s_Level_Rad_UP(1:RTV%n_Angles,n_Layers ) = emissivity(1:RTV%n_Angles)*RTV%Planck_Surface END IF IF( RTV%Solar_Flag_true ) THEN RTV%s_Level_Rad_UP(1:nZ,n_Layers ) = RTV%s_Level_Rad_UP(1:nZ,n_Layers )+direct_reflectivity(1:nZ)* & RTV%COS_SUN*RTV%Solar_irradiance/PI*exp(-total_opt(n_Layers)/RTV%COS_SUN) END IF ! UPWARD ADDING LOOP STARTS FROM BOTTOM LAYER TO ATMOSPHERIC TOP LAYER. DO 10 k = n_Layers, 1, -1 ! Compute tranmission and reflection matrices for a layer IF(w(k) > SCATTERING_ALBEDO_THRESHOLD) THEN ! ----------------------------------------------------------- ! ! CALL multiple-stream algorithm for computing layer ! ! transmission, reflection, and source functions. ! ! ----------------------------------------------------------- ! CALL CRTM_AMOM_layer( & RTV%n_Streams, & RTV%n_Angles,k,w(k), & T_OD(k), & total_opt(k-1), & RTV%COS_Angle, & ! Input RTV%COS_Weight, & RTV%Pff(:,:,k), & RTV%Pbb(:,:,k), & ! Input RTV%Planck_Atmosphere(k), & ! Input RTV ) ! Internal variable ! ----------------------------------------------------------- ! ! Adding method to add the layer to the present level ! ! to compute upward radiances and reflection matrix ! ! at new level. ! ! ----------------------------------------------------------- ! temporal_matrix = -matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k), & RTV%s_Layer_Refl(1:RTV%n_Angles,1:RTV%n_Angles,k)) DO i = 1, RTV%n_Angles temporal_matrix(i,i) = ONE + temporal_matrix(i,i) END DO RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k) = matinv(temporal_matrix, Error_Status) IF( Error_Status /= SUCCESS ) THEN WRITE( Message,'("Error in matrix inversion matinv(temporal_matrix, Error_Status) ")' ) CALL Display_Message( ROUTINE_NAME, & TRIM(Message), & Error_Status ) RETURN END IF RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k) = & 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)) refl_down(1:RTV%n_Angles,k) = matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k), & RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k)) RTV%s_Level_Rad_UP(1:RTV%n_Angles,k-1 )=RTV%s_Layer_Source_UP(1:RTV%n_Angles,k)+ & matmul(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k),refl_down(1:RTV%n_Angles,k) & +RTV%s_Level_Rad_UP(1:RTV%n_Angles,k )) 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), & RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k)) 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) + & 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)) ELSE DO i = 1, RTV%n_Angles RTV%s_Layer_Trans(i,i,k) = exp(-T_OD(k)/RTV%COS_Angle(i)) RTV%s_Layer_Source_UP(i,k) = RTV%Planck_Atmosphere(k) * (ONE - RTV%s_Layer_Trans(i,i,k) ) RTV%s_Layer_Source_DOWN(i,k) = RTV%s_Layer_Source_UP(i,k) END DO ! Adding method DO i = 1, RTV%n_Angles RTV%s_Level_Rad_UP(i,k-1 )=RTV%s_Layer_Source_UP(i,k)+ & 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)) & +RTV%s_Level_Rad_UP(i,k )) ENDDO DO i = 1, RTV%n_Angles DO j = 1, RTV%n_Angles 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) ENDDO ENDDO ENDIF 10 CONTINUE ! Adding reflected cosmic background radiation IF( RTV%mth_Azi == 0 ) THEN DO i = 1, RTV%n_Angles 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 ENDDO END IF RETURN END SUBROUTINE CRTM_ADA SUBROUTINE CRTM_AMOM_layer( n_streams, & ! Input, number of streams nZ, & ! Input, number of angles KL, & ! Input, KL-th layer single_albedo, & ! Input, single scattering albedo optical_depth, & ! Input, layer optical depth total_opt, & ! Input, accumulated optical depth from the top to current layer top COS_Angle, & ! Input, COSINE of ANGLES COS_Weight, & ! Input, GAUSSIAN Weights ff, & ! Input, Phase matrix (forward part) bb, & ! Input, Phase matrix (backward part) Planck_Func, & ! Input, Planck for layer temperature RTV) ! Output, layer transmittance, reflectance, and source ! --------------------------------------------------------------------------------------- ! FUNCTION ! Compute layer transmission, reflection matrices and source function ! at the top and bottom of the layer. ! ! Method and References ! The transmittance and reflectance matrices is further derived from ! matrix operator method. The matrix operator method is referred to the paper by ! ! Weng, F., and Q. Liu, 2003: Satellite Data Assimilation in Numerical Weather Prediction ! Model: Part 1: Forward Radiative Transfer and Jacobian Modeling in Cloudy Atmospheres, ! J. Atmos. Sci., 60, 2633-2646. ! ! see also ADA method. ! Quanhua Liu ! Quanhua.Liu@noaa.gov ! ---------------------------------------------------------------------------------------- IMPLICIT NONE INTEGER, INTENT(IN) :: n_streams,nZ,KL TYPE(RTV_type), INTENT( INOUT ) :: RTV REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt ! internal variables REAL(fp), DIMENSION(nZ,nZ) :: trans, refl, tempo REAL(fp) :: s, c, xx INTEGER :: i,j,N2,N2_1 INTEGER :: Error_Status REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ) REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ) CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer' CHARACTER(256) :: Message ! for small layer optical depth, single scattering is applied. IF( optical_depth < DELTA_OPTICAL_DEPTH ) THEN s = optical_depth * single_albedo DO i = 1, nZ RTV%Thermal_C(i,KL) = ZERO c = s/COS_Angle(i) DO j = 1, nZ RTV%s_Layer_Refl(i,j,KL) = c * bb(i,j) * COS_Weight(j) RTV%s_Layer_Trans(i,j,KL) = c * ff(i,j) * COS_Weight(j) IF( i == j ) THEN RTV%s_Layer_Trans(i,i,KL) = RTV%s_Layer_Trans(i,i,KL) + & ONE - optical_depth/COS_Angle(i) END IF IF( RTV%mth_Azi == 0 ) THEN RTV%Thermal_C(i,KL) = RTV%Thermal_C(i,KL) + & ( RTV%s_Layer_Refl(i,j,KL) + RTV%s_Layer_Trans(i,j,KL) ) END IF ENDDO IF( RTV%mth_Azi == 0 ) THEN RTV%s_Layer_Source_UP(i,KL) = ( ONE - RTV%Thermal_C(i,KL) ) * Planck_Func RTV%s_Layer_Source_DOWN(i,KL) = RTV%s_Layer_Source_UP(i,KL) END IF ENDDO RETURN END IF ! ! for numerical stability, IF( single_albedo < max_albedo ) THEN s = single_albedo ELSE s = max_albedo END IF ! ! building phase matrices DO i = 1, nZ c = s/COS_Angle(i) DO j = 1, nZ RTV%PM(i,j,KL) = c * bb(i,j) * COS_Weight(j) RTV%PP(i,j,KL) = c * ff(i,j) * COS_Weight(j) ENDDO RTV%PP(i,i,KL) = RTV%PP(i,i,KL) - ONE/COS_Angle(i) ENDDO RTV%PPM(1:nZ,1:nZ,KL) = RTV%PP(1:nZ,1:nZ,KL) - RTV%PM(1:nZ,1:nZ,KL) RTV%i_PPM(1:nZ,1:nZ,KL) = matinv( RTV%PPM(1:nZ,1:nZ,KL), Error_Status ) IF( Error_Status /= SUCCESS ) THEN WRITE( Message,'("Error in matrix inversion matinv( RTV%PPM(1:nZ,1:nZ,KL), Error_Status ) ")' ) CALL Display_Message( ROUTINE_NAME, & TRIM(Message), & Error_Status ) RETURN END IF RTV%PPP(1:nZ,1:nZ,KL) = RTV%PP(1:nZ,1:nZ,KL) + RTV%PM(1:nZ,1:nZ,KL) RTV%HH(1:nZ,1:nZ,KL) = matmul( RTV%PPM(1:nZ,1:nZ,KL), RTV%PPP(1:nZ,1:nZ,KL) ) ! ! save phase element RTV%HH, call ASYMTX for calculating eigenvalue and vectors. tempo = RTV%HH(1:nZ,1:nZ,KL) CALL ASYMTX(tempo,nZ,nZ,nZ,RTV%EigVe(1:nZ,1:nZ,KL),RTV%EigVa(1:nZ,KL),Error_Status) DO i = 1, nZ IF( RTV%EigVa(i,KL) > ZERO ) THEN RTV%EigValue(i,KL) = sqrt( RTV%EigVa(i,KL) ) ELSE RTV%EigValue(i,KL) = ZERO END IF END DO DO i = 1, nZ DO j = 1, nZ RTV%EigVeVa(i,j,KL) = RTV%EigVe(i,j,KL) * RTV%EigValue(j,KL) END DO END DO RTV%EigVeF(1:nZ,1:nZ,KL) = matmul( RTV%i_PPM(1:nZ,1:nZ,KL), RTV%EigVeVa(1:nZ,1:nZ,KL) ) ! compute layer reflection, transmission and source function RTV%Gp(1:nZ,1:nZ,KL) = ( RTV%EigVe(1:nZ,1:nZ,KL) + RTV%EigVeF(1:nZ,1:nZ,KL) )/2.0_fp RTV%Gm(1:nZ,1:nZ,KL) = ( RTV%EigVe(1:nZ,1:nZ,KL) - RTV%EigVeF(1:nZ,1:nZ,KL) )/2.0_fp RTV%i_Gm(1:nZ,1:nZ,KL) = matinv( RTV%Gm(1:nZ,1:nZ,KL), Error_Status) IF( Error_Status /= SUCCESS ) THEN WRITE( Message,'("Error in matrix inversion matinv( RTV%Gm(1:nZ,1:nZ,KL), Error_Status) ")' ) CALL Display_Message( ROUTINE_NAME, & TRIM(Message), & Error_Status ) RETURN END IF DO i = 1, nZ xx = RTV%EigValue(i,KL)*optical_depth RTV%Exp_x(i,KL) = exp(-xx) END DO DO i = 1, nZ DO j = 1, nZ RTV%A1(i,j,KL) = RTV%Gp(i,j,KL) * RTV%Exp_x(j,KL) RTV%A4(i,j,KL) = RTV%Gm(i,j,KL) * RTV%Exp_x(j,KL) END DO END DO RTV%A2(1:nZ,1:nZ,KL) = matmul( RTV%i_Gm(1:nZ,1:nZ,KL), RTV%A1(1:nZ,1:nZ,KL) ) RTV%A3(1:nZ,1:nZ,KL) = matmul( RTV%Gp(1:nZ,1:nZ,KL), RTV%A2(1:nZ,1:nZ,KL) ) RTV%A5(1:nZ,1:nZ,KL) = matmul( RTV%A1(1:nZ,1:nZ,KL), RTV%A2(1:nZ,1:nZ,KL) ) RTV%A6(1:nZ,1:nZ,KL) = matmul( RTV%A4(1:nZ,1:nZ,KL), RTV%A2(1:nZ,1:nZ,KL) ) RTV%Gm_A5(1:nZ,1:nZ,KL) = RTV%Gm(1:nZ,1:nZ,KL) - RTV%A5(1:nZ,1:nZ,KL) RTV%i_Gm_A5(1:nZ,1:nZ,KL) = matinv(RTV%Gm_A5(1:nZ,1:nZ,KL), Error_Status) IF( Error_Status /= SUCCESS ) THEN WRITE( Message,'("Error in matrix inversion matinv(RTV%Gm_A5(1:nZ,1:nZ,KL), Error_Status) ")' ) CALL Display_Message( ROUTINE_NAME, & TRIM(Message), & Error_Status ) RETURN END IF 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) ) 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) ) ! post processing RTV%s_Layer_Trans(1:nZ,1:nZ,KL) = trans(:,:) RTV%s_Layer_Refl(1:nZ,1:nZ,KL) = refl(:,:) RTV%s_Layer_Source_UP(:,KL) = ZERO IF( RTV%mth_Azi == 0 ) THEN DO i = 1, nZ RTV%Thermal_C(i,KL) = ZERO DO j = 1, n_Streams RTV%Thermal_C(i,KL) = RTV%Thermal_C(i,KL) + (trans(i,j) + refl(i,j) ) END DO IF ( i == nZ .AND. nZ == (n_Streams+1) ) THEN RTV%Thermal_C(i,KL) = RTV%Thermal_C(i,KL) + trans(nZ,nZ) END IF RTV%s_Layer_Source_UP(i,KL) = ( ONE - RTV%Thermal_C(i,KL) ) * Planck_Func RTV%s_Layer_Source_DOWN(i,KL) = RTV%s_Layer_Source_UP(i,KL) END DO END IF ! compute visible part for visible channels during daytime IF( RTV%Solar_Flag_true ) THEN N2 = 2 * nZ N2_1 = N2 - 1 source_up = ZERO source_down = ZERO ! ! Solar source Sfactor = single_albedo*RTV%Solar_irradiance/PI IF( RTV%mth_Azi == 0 ) Sfactor = Sfactor/TWO EXPfactor = exp(-optical_depth/RTV%COS_SUN) s_transmittance = exp(-total_opt/RTV%COS_SUN) DO i = 1, nZ Solar(i) = -bb(i,nZ+1)*Sfactor Solar(i+nZ) = -ff(i,nZ+1)*Sfactor DO j = 1, nZ V0(i,j) = single_albedo * ff(i,j) * COS_Weight(j) V0(i+nZ,j) = single_albedo * bb(i,j) * COS_Weight(j) V0(i,j+nZ) = V0(i+nZ,j) V0(nZ+i,j+nZ) = V0(i,j) ENDDO V0(i,i) = V0(i,i) - ONE - COS_Angle(i)/RTV%COS_SUN V0(i+nZ,i+nZ) = V0(i+nZ,i+nZ) - ONE + COS_Angle(i)/RTV%COS_SUN ENDDO V1(1:N2_1,1:N2_1) = matinv(V0(1:N2_1,1:N2_1), Error_Status) IF( Error_Status /= SUCCESS ) THEN WRITE( Message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' ) CALL Display_Message( ROUTINE_NAME, & TRIM(Message), & Error_Status ) RETURN END IF Solar1(1:N2_1) = matmul( V1(1:N2_1,1:N2_1), Solar(1:N2_1) ) Solar1(N2) = ZERO Sfac2 = Solar(N2) - sum( V0(N2,1:N2_1)*Solar1(1:N2_1) ) DO i = 1, nZ source_up(i) = Solar1(i) source_down(i) = EXPfactor*Solar1(i+nZ) DO j = 1, nZ source_up(i) =source_up(i)-refl(i,j)*Solar1(j+nZ)-trans(i,j)*EXPfactor*Solar1(j) source_down(i) =source_down(i) -trans(i,j)*Solar1(j+nZ) -refl(i,j)*EXPfactor*Solar1(j) END DO END DO ! specific treatment for downeward source function IF( abs( V0(N2,N2) ) > 0.0001_fp ) THEN source_down(nZ) =source_down(nZ) +(EXPfactor-trans(nZ,nZ))*Sfac2/V0(N2,N2) ELSE source_down(nZ) =source_down(nZ) -EXPfactor*Sfac2*optical_depth/COS_Angle(nZ) END IF source_up(1:nZ) = source_up(1:nZ)*s_transmittance source_down(1:nZ) = source_down(1:nZ)*s_transmittance RTV%s_Layer_Source_UP(1:nZ,KL) = RTV%s_Layer_Source_UP(1:nZ,KL)+source_up(1:nZ) RTV%s_Layer_Source_DOWN(1:nZ,KL) = RTV%s_Layer_Source_DOWN(1:nZ,KL)+source_down(1:nZ) END IF RETURN END SUBROUTINE CRTM_AMOM_layer SUBROUTINE CRTM_ADA_TL(n_Layers, & ! Input number of atmospheric layers w, & ! Input layer scattering albedo T_OD, & ! Input layer optical depth cosmic_background, & ! Input cosmic background radiance emissivity, & ! Input surface emissivity direct_reflectivity, & ! Input direct reflectivity RTV, & ! Input structure containing forward part results Planck_Atmosphere_TL, & ! Input tangent-linear atmospheric layer Planck radiance Planck_Surface_TL, & ! Input TL surface Planck radiance w_TL, & ! Input TL layer scattering albedo T_OD_TL, & ! Input TL layer optical depth emissivity_TL, & ! Input TL surface emissivity reflectivity_TL, & ! Input TL reflectivity direct_reflectivity_TL, & ! Input TL direct reflectivity Pff_TL, & ! Input TL forward phase matrix Pbb_TL, & ! Input TL backward phase matrix s_rad_up_TL) ! Output TL upward radiance ! ------------------------------------------------------------------------- ! ! FUNCTION: ! ! This subroutine calculates IR/MW tangent-linear radiance at the top of ! ! the atmosphere including atmospheric scattering. The structure RTV ! ! carried in forward part results. ! ! The CRTM_ADA_TL algorithm computes layer tangent-linear reflectance and ! ! transmittance as well as source function by the subroutine ! ! CRTM_Doubling_layer as source function by the subroutine ! ! CRTM_Doubling_layer, then uses ! ! an adding method to integrate the layer and surface components. ! ! ! ! Quanhua Liu Quanhua.Liu@noaa.gov ! ! ------------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, INTENT(IN) :: n_Layers TYPE(RTV_type), INTENT(IN) :: RTV REAL (fp), INTENT(IN), DIMENSION( : ) :: w,T_OD REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity,direct_reflectivity REAL (fp), INTENT(IN) :: cosmic_background REAL (fp),INTENT(IN),DIMENSION( :,:,: ) :: Pff_TL, Pbb_TL REAL (fp),INTENT(IN),DIMENSION( : ) :: w_TL,T_OD_TL REAL (fp),INTENT(IN),DIMENSION( 0: ) :: Planck_Atmosphere_TL REAL (fp),INTENT(IN) :: Planck_Surface_TL REAL (fp),INTENT(IN),DIMENSION( : ) :: emissivity_TL REAL (fp),INTENT(IN),DIMENSION( :,: ) :: reflectivity_TL REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_TL REAL (fp),INTENT(INOUT),DIMENSION( : ) :: direct_reflectivity_TL ! -------------- internal variables --------------------------------- ! ! Abbreviations: ! ! s: scattering, rad: radiance, trans: transmission, ! ! refl: reflection, up: upward, down: downward ! ! --------------------------------------------------------------------! REAL (fp), DIMENSION(RTV%n_Angles,RTV%n_Angles) :: temporal_matrix_TL REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down REAL (fp), DIMENSION( RTV%n_Angles ) :: s_source_up_TL,s_source_down_TL,refl_down_TL REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles ) :: s_trans_TL,s_refl_TL,Refl_Trans_TL REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles ) :: s_refl_up_TL,Inv_Gamma_TL,Inv_GammaT_TL REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_TL INTEGER :: i, j, k, nZ ! nZ = RTV%n_Angles total_opt(0) = ZERO total_opt_TL(0) = ZERO DO k = 1, n_Layers total_opt(k) = total_opt(k-1) + T_OD(k) total_opt_TL(k) = total_opt_TL(k-1) + T_OD_TL(k) END DO Refl_Trans_TL = ZERO s_rad_up_TL = ZERO s_refl_up_TL = reflectivity_TL IF( RTV%mth_Azi == 0 ) THEN s_rad_up_TL = emissivity_TL * RTV%Planck_Surface + emissivity * Planck_Surface_TL END IF IF( RTV%Solar_Flag_true ) THEN s_rad_up_TL = s_rad_up_TL+direct_reflectivity_TL*RTV%COS_SUN*RTV%Solar_irradiance/PI & * exp(-total_opt(n_Layers)/RTV%COS_SUN) & - direct_reflectivity * RTV%Solar_irradiance/PI & * total_opt_TL(n_Layers) * exp(-total_opt(n_Layers)/RTV%COS_SUN) END IF DO 10 k = n_Layers, 1, -1 s_source_up_TL = ZERO s_source_down_TL = ZERO s_trans_TL = ZERO s_refl_TL = ZERO Inv_GammaT_TL = ZERO Inv_Gamma_TL = ZERO refl_down_TL = ZERO ! ! Compute tranmission and reflection matrices for a layer IF(w(k) > SCATTERING_ALBEDO_THRESHOLD) THEN ! ----------------------------------------------------------- ! ! CALL Doubling algorithm to computing forward and tagent ! ! layer transmission, reflection, and source functions. ! ! ----------------------------------------------------------- ! call CRTM_AMOM_layer_TL(RTV%n_Streams,RTV%n_Angles,k,w(k),T_OD(k),total_opt(k-1), & !Input RTV%COS_Angle(1:RTV%n_Angles),RTV%COS_Weight(1:RTV%n_Angles) , & !Input RTV%Pff(:,:,k), RTV%Pbb(:,:,k),RTV%Planck_Atmosphere(k) , & !Input w_TL(k),T_OD_TL(k),total_opt_TL(k-1),Pff_TL(:,:,k) , & !Input Pbb_TL(:,:,k),Planck_Atmosphere_TL(k),RTV , & !Input s_trans_TL,s_refl_TL,s_source_up_TL,s_source_down_TL) !Output ! Adding method temporal_matrix_TL = -matmul(s_refl_up_TL,RTV%s_Layer_Refl(1:RTV%n_Angles,1:RTV%n_Angles,k)) & - matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),s_refl_TL) temporal_matrix_TL = matmul(RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k),temporal_matrix_TL) Inv_Gamma_TL = -matmul(temporal_matrix_TL,RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k)) Inv_GammaT_TL = matmul(s_trans_TL, RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k)) & + matmul(RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k), Inv_Gamma_TL) refl_down(1:RTV%n_Angles,k) = matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k), & RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k)) refl_down_TL(:) = matmul(s_refl_up_TL,RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k)) & + matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k),s_source_down_TL(:)) s_rad_up_TL(1:RTV%n_Angles)=s_source_up_TL(1:RTV%n_Angles)+ & matmul(Inv_GammaT_TL,refl_down(:,k)+RTV%s_Level_Rad_UP(1:RTV%n_Angles,k)) & +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)) Refl_Trans_TL = matmul(s_refl_up_TL,RTV%s_Layer_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),s_trans_TL) s_refl_up_TL=s_refl_TL+matmul(Inv_GammaT_TL,RTV%Refl_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k)) & +matmul(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k),Refl_Trans_TL) Refl_Trans_TL = ZERO ELSE DO i = 1, RTV%n_Angles s_trans_TL(i,i) = -T_OD_TL(k)/RTV%COS_Angle(i) * RTV%s_Layer_Trans(i,i,k) s_source_up_TL(i) = Planck_Atmosphere_TL(k) * (ONE - RTV%s_Layer_Trans(i,i,k) ) & - RTV%Planck_Atmosphere(k) * s_trans_TL(i,i) s_source_down_TL(i) = s_source_up_TL(i) ENDDO ! Adding method DO i = 1, RTV%n_Angles s_rad_up_TL(i)=s_source_up_TL(i) & +s_trans_TL(i,i)*(sum(RTV%s_Level_Refl_UP(i,1:RTV%n_Angles,k) & *RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k))+RTV%s_Level_Rad_UP(i,k)) & +RTV%s_Layer_Trans(i,i,k) & *(sum(s_refl_up_TL(i,1:RTV%n_Angles)*RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k) & +RTV%s_Level_Refl_UP(i,1:RTV%n_Angles,k)*s_source_down_TL(1:RTV%n_Angles))+s_rad_up_TL(i)) ENDDO DO i = 1, RTV%n_Angles DO j = 1, RTV%n_Angles s_refl_up_TL(i,j)=s_trans_TL(i,i)*RTV%s_Level_Refl_UP(i,j,k) & *RTV%s_Layer_Trans(j,j,k) & +RTV%s_Layer_Trans(i,i,k)*s_refl_up_TL(i,j)*RTV%s_Layer_Trans(j,j,k) & +RTV%s_Layer_Trans(i,i,k)*RTV%s_Level_Refl_UP(i,j,k)*s_trans_TL(j,j) ENDDO ENDDO ENDIF 10 CONTINUE ! ! Adding reflected cosmic background radiation IF( RTV%mth_Azi == 0 ) THEN DO i = 1, RTV%n_Angles s_rad_up_TL(i)=s_rad_up_TL(i)+sum(s_refl_up_TL(i,:))*cosmic_background ENDDO END IF ! RETURN END SUBROUTINE CRTM_ADA_TL SUBROUTINE CRTM_AMOM_layer_TL( n_streams, & ! Input, number of streams nZ, & ! Input, number of angles KL, & ! Input, KL-th layer single_albedo, & ! Input, single scattering albedo optical_depth, & ! Input, layer optical depth total_opt, & ! Input, accumulated optical depth from the top to current layer top COS_Angle, & ! Input, COSINE of ANGLES COS_Weight, & ! Input, GAUSSIAN Weights ff, & ! Input, Phase matrix (forward part) bb, & ! Input, Phase matrix (backward part) Planck_Func, & ! Input, Planck for layer temperature single_albedo_TL, & ! Input, tangent-linear single albedo optical_depth_TL, & ! Input, TL layer optical depth total_opt_TL, & ! Input, accumulated TL optical depth from the top to current layer top ff_TL, & ! Input, TL forward Phase matrix bb_TL, & ! Input, TL backward Phase matrix Planck_Func_TL, & ! Input, TL Planck for layer temperature RTV, & ! Input, structure containing forward results trans_TL, & ! Output, layer tangent-linear trans refl_TL, & ! Output, layer tangent-linear refl source_up_TL, & ! Output, layer tangent-linear source_up source_down_TL) ! Output, layer tangent-linear source_down ! --------------------------------------------------------------------------------------- ! FUNCTION ! Compute TL layer transmission, reflection matrices and source function ! at the top and bottom of the layer. ! ! see also ADA method. ! Quanhua Liu ! Quanhua.Liu@noaa.gov ! ---------------------------------------------------------------------------------------- IMPLICIT NONE INTEGER, INTENT(IN) :: n_streams,nZ,KL TYPE(RTV_type), INTENT( IN ) :: RTV REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt ! ! internal variables REAL(fp) :: s, c, s_TL,c_TL,xx_TL INTEGER :: i,j INTEGER :: Error_Status ! ! Tangent-Linear Part REAL(fp), INTENT(OUT), DIMENSION( :,: ) :: trans_TL,refl_TL REAL(fp), INTENT(OUT), DIMENSION( : ) :: source_up_TL,source_down_TL REAL(fp), INTENT(IN) :: single_albedo_TL REAL(fp), INTENT(IN) :: optical_depth_TL,Planck_Func_TL,total_opt_TL REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff_TL,bb_TL REAL(fp), DIMENSION(nZ) :: Exp_x_TL,EigVa_TL,EigValue_TL REAL(fp), DIMENSION(nZ,nZ) :: i_Gm_TL, Gm_TL, Gp_TL, EigVe_TL, EigVeF_TL, EigVeVa_TL REAL(fp), DIMENSION(nZ,nZ) :: A1_TL,A2_TL,A3_TL,A4_TL,A5_TL,A6_TL,Gm_A5_TL,i_Gm_A5_TL REAL(fp), DIMENSION(nZ,nZ) :: HH_TL,PM_TL,PP_TL,PPM_TL,i_PPM_TL,PPP_TL REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ) REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ) REAL(fp) :: EXPfactor_TL,Sfactor_TL,s_transmittance_TL,Solar_TL(2*nZ),V0_TL(2*nZ,2*nZ),Solar1_TL(2*nZ) REAL(fp) :: Sfac2_TL REAL(fp), DIMENSION( nZ ) :: thermal_up_TL,thermal_down_TL REAL(fp), DIMENSION(nZ,nZ) :: trans, refl INTEGER :: N2, N2_1 REAL(fp) :: Thermal_C_TL CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer_TL' CHARACTER(256) :: Message ! ! for small layer optical depth, single scattering is applied. IF( optical_depth < DELTA_OPTICAL_DEPTH ) THEN s = optical_depth * single_albedo s_TL = optical_depth_TL * single_albedo + optical_depth * single_albedo_TL DO i = 1, nZ Thermal_C_TL = ZERO c = s/COS_Angle(i) c_TL = s_TL/COS_Angle(i) DO j = 1, nZ refl_TL(i,j) = (c_TL * bb(i,j) + c*bb_TL(i,j) )* COS_Weight(j) trans_TL(i,j) = (c_TL * ff(i,j) + c*ff_TL(i,j) )* COS_Weight(j) IF( i == j ) THEN trans_TL(i,j) = trans_TL(i,j) - optical_depth_TL/COS_Angle(i) END IF IF( RTV%mth_Azi == 0 ) THEN Thermal_C_TL = Thermal_C_TL + refl_TL(i,j) + trans_TL(i,j) END IF ENDDO IF( RTV%mth_Azi == 0 ) THEN source_up_TL(i) = -Thermal_C_TL * Planck_Func + & ( ONE - RTV%Thermal_C(i,KL) ) * Planck_Func_TL source_down_TL(i) = source_up_TL(i) END IF ENDDO RETURN END IF ! ! for numerical stability, IF( single_albedo < max_albedo ) THEN s = single_albedo s_TL = single_albedo_TL ELSE s = max_albedo s_TL = 0.0_fp END IF ! ! building TL phase matrices DO i = 1, nZ c = s/COS_Angle(i) c_TL = s_TL/COS_Angle(i) DO j = 1, nZ PM_TL(i,j) = (c_TL * bb(i,j) + c*bb_TL(i,j) ) * COS_Weight(j) PP_TL(i,j) = (c_TL * ff(i,j) + c*ff_TL(i,j) ) * COS_Weight(j) END DO ENDDO PPM_TL(:,:) = PP_TL(:,:) - PM_TL(:,:) i_PPM_TL(:,:) = - matmul( RTV%i_PPM(1:nZ,1:nZ,KL), matmul(PPM_TL(:,:),RTV%i_PPM(1:nZ,1:nZ,KL)) ) PPP_TL(:,:) = PP_TL(:,:) + PM_TL(:,:) HH_TL(:,:) = matmul( PPM_TL(:,:), RTV%PPP(1:nZ,1:nZ,KL) )+matmul( RTV%PPM(1:nZ,1:nZ,KL), PPP_TL(:,:) ) ! ! compute TL eigenvectors EigVe, and eigenvalues EigVa CALL ASYMTX_TL(nZ,RTV%EigVe(1:nZ,1:nZ,KL),RTV%EigVa(1:nZ,KL),HH_TL, & EigVe_TL,EigVa_TL,Error_Status) DO i = 1, nZ IF( RTV%EigVa(i,KL) > ZERO ) THEN EigValue_TL(i) = 0.5_fp*EigVa_TL(i)/RTV%EigValue(i,KL) ELSE EigValue_TL(i) = ZERO END IF END DO EigVeVa_TL = ZERO DO i = 1, nZ DO j = 1, nZ EigVeVa_TL(i,j) = EigVe_TL(i,j) * RTV%EigValue(j,KL)+RTV%EigVe(i,j,KL) * EigValue_TL(j) END DO END DO EigVeF_TL(:,:) = matmul( i_PPM_TL(:,:), RTV%EigVeVa(1:nZ,1:nZ,KL) ) & + matmul( RTV%i_PPM(1:nZ,1:nZ,KL), EigVeVa_TL(:,:) ) ! ! compute TL reflection and transmission matrices, TL source function Gp_TL(:,:) = ( EigVe_TL(:,:) + EigVeF_TL(:,:) )/2.0_fp Gm_TL(:,:) = ( EigVe_TL(:,:) - EigVeF_TL(:,:) )/2.0_fp i_Gm_TL = -matmul( RTV%i_Gm(1:nZ,1:nZ,KL), matmul(Gm_TL,RTV%i_Gm(1:nZ,1:nZ,KL)) ) DO i = 1, nZ xx_TL = EigValue_TL(i)*optical_depth+RTV%EigValue(i,KL)*optical_depth_TL Exp_x_TL(i) = -xx_TL*RTV%Exp_x(i,KL) END DO DO i = 1, nZ DO j = 1, nZ A1_TL(i,j) = Gp_TL(i,j)* RTV%Exp_x(j,KL)+ RTV%Gp(i,j,KL)* Exp_x_TL(j) A4_TL(i,j) = Gm_TL(i,j)* RTV%Exp_x(j,KL)+ RTV%Gm(i,j,KL)* Exp_x_TL(j) END DO END DO A2_TL(:,:) = matmul(i_Gm_TL(:,:),RTV%A1(1:nZ,1:nZ,KL))+matmul(RTV%i_Gm(1:nZ,1:nZ,KL),A1_TL(:,:)) A3_TL(:,:) = matmul(Gp_TL(:,:),RTV%A2(1:nZ,1:nZ,KL))+matmul(RTV%Gp(1:nZ,1:nZ,KL),A2_TL(:,:)) A5_TL(:,:) = matmul(A1_TL(:,:),RTV%A2(1:nZ,1:nZ,KL))+matmul(RTV%A1(1:nZ,1:nZ,KL),A2_TL(:,:)) A6_TL(:,:) = matmul(A4_TL(:,:),RTV%A2(1:nZ,1:nZ,KL))+matmul(RTV%A4(1:nZ,1:nZ,KL),A2_TL(:,:)) Gm_A5_TL(:,:) = Gm_TL(:,:) - A5_TL(:,:) 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))) ! ! T = matmul( RTV%A4(:,:,KL) - RTV%A3(:,:,KL), RTV%i_Gm_A5(:,:,KL) ) trans_TL = matmul( A4_TL(:,:) - A3_TL(:,:), RTV%i_Gm_A5(1:nZ,1:nZ,KL) ) & + matmul( RTV%A4(1:nZ,1:nZ,KL) - RTV%A3(1:nZ,1:nZ,KL), i_Gm_A5_TL(:,:) ) refl_TL = matmul( Gp_TL(:,:) - A6_TL(:,:), RTV%i_Gm_A5(1:nZ,1:nZ,KL) ) & + matmul( RTV%Gp(1:nZ,1:nZ,KL) - RTV%A6(1:nZ,1:nZ,KL), i_Gm_A5_TL(:,:) ) trans(1:nZ,1:nZ) = RTV%s_Layer_Trans(1:nZ,1:nZ,KL) refl(1:nZ,1:nZ) = RTV%s_Layer_Refl(1:nZ,1:nZ,KL) Source_UP_TL = ZERO Source_DOWN_TL = ZERO ! ! Thermal part IF( RTV%mth_Azi == 0 ) THEN DO i = 1, nZ Thermal_C_TL = ZERO DO j = 1, n_Streams Thermal_C_TL = Thermal_C_TL + (trans_TL(i,j) + refl_TL(i,j)) ENDDO IF(i == nZ .AND. nZ == (n_Streams+1)) THEN Thermal_C_TL = Thermal_C_TL + trans_TL(nZ,nZ) ENDIF thermal_up_TL(i) = -Thermal_C_TL * Planck_Func & + ( ONE - RTV%Thermal_C(i,KL) ) * Planck_Func_TL thermal_down_TL(i) = thermal_up_TL(i) ENDDO END IF ! ! for visible channels at daytime IF( RTV%Solar_Flag_true ) THEN N2 = 2 * nZ N2_1 = N2 - 1 V0 = ZERO V1 = ZERO Solar = ZERO Solar1 = ZERO Sfac2 = ZERO V0_TL = ZERO Solar_TL = ZERO Solar1_TL = ZERO Sfac2_TL = ZERO ! ! Solar source Sfactor = single_albedo*RTV%Solar_irradiance/PI Sfactor_TL = single_albedo_TL*RTV%Solar_irradiance/PI IF( RTV%mth_Azi == 0 ) THEN Sfactor = Sfactor/TWO Sfactor_TL = Sfactor_TL/TWO END IF EXPfactor = exp(-optical_depth/RTV%COS_SUN) EXPfactor_TL = -optical_depth_TL/RTV%COS_SUN*EXPfactor s_transmittance = exp(-total_opt/RTV%COS_SUN) s_transmittance_TL = -total_opt_TL/RTV%COS_SUN*s_transmittance DO i = 1, nZ Solar(i) = -bb(i,nZ+1)*Sfactor Solar_TL(i) = -bb_TL(i,nZ+1)*Sfactor-bb(i,nZ+1)*Sfactor_TL Solar(i+nZ) = -ff(i,nZ+1)*Sfactor Solar_TL(i+nZ) = -ff_TL(i,nZ+1)*Sfactor-ff(i,nZ+1)*Sfactor_TL DO j = 1, nZ V0(i,j) = single_albedo * ff(i,j) * COS_Weight(j) V0_TL(i,j) = single_albedo_TL*ff(i,j)*COS_Weight(j)+single_albedo*ff_TL(i,j)*COS_Weight(j) V0(i+nZ,j) = single_albedo * bb(i,j) * COS_Weight(j) V0_TL(i+nZ,j) = single_albedo_TL*bb(i,j)*COS_Weight(j)+single_albedo*bb_TL(i,j)*COS_Weight(j) V0(i,j+nZ) = V0(i+nZ,j) V0_TL(i,j+nZ) = V0_TL(i+nZ,j) V0(nZ+i,j+nZ) = V0(i,j) V0_TL(nZ+i,j+nZ) = V0_TL(i,j) ENDDO V0(i,i) = V0(i,i) - ONE - COS_Angle(i)/RTV%COS_SUN V0(i+nZ,i+nZ) = V0(i+nZ,i+nZ) - ONE + COS_Angle(i)/RTV%COS_SUN ENDDO V1(1:N2_1,1:N2_1) = matinv(V0(1:N2_1,1:N2_1), Error_Status) IF( Error_Status /= SUCCESS ) THEN WRITE( Message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' ) CALL Display_Message( ROUTINE_NAME, & TRIM(Message), & Error_Status ) RETURN END IF Solar1(1:N2_1) = matmul( V1(1:N2_1,1:N2_1), Solar(1:N2_1) ) Solar(1:N2_1) = matmul( V1(1:N2_1,1:N2_1),Solar(1:N2_1) ) Solar1_TL(1:N2_1) = matmul( V0_TL(1:N2_1,1:N2_1),Solar(1:N2_1) ) Solar1_TL(1:N2_1) = -matmul( V1(1:N2_1,1:N2_1),Solar1_TL(1:N2_1) ) & + matmul( V1(1:N2_1,1:N2_1), Solar_TL(1:N2_1) ) Solar1(N2) = ZERO Solar1_TL(N2) = ZERO Sfac2 = Solar(N2) - sum( V0(N2,1:N2_1)*Solar1(1:N2_1) ) Sfac2_TL = Solar_TL(N2) - sum( V0_TL(N2,1:N2_1)*Solar1(1:N2_1) ) & - sum( V0(N2,1:N2_1)*Solar1_TL(1:N2_1) ) DO i = 1, nZ source_up(i) = Solar1(i) source_up_TL(i) = Solar1_TL(i) source_down(i) = EXPfactor*Solar1(i+nZ) source_down_TL(i) = EXPfactor_TL*Solar1(i+nZ)+EXPfactor*Solar1_TL(i+nZ) DO j = 1, nZ source_up(i) =source_up(i)-refl(i,j)*Solar1(j+nZ)-trans(i,j)*EXPfactor*Solar1(j) source_up_TL(i) =source_up_TL(i)-refl_TL(i,j)*Solar1(j+nZ) -refl(i,j)*Solar1_TL(j+nZ) & - trans_TL(i,j)*EXPfactor*Solar1(j) - trans(i,j)*EXPfactor_TL*Solar1(j) -trans(i,j)*EXPfactor*Solar1_TL(j) source_down(i) =source_down(i) -trans(i,j)*Solar1(j+nZ) -refl(i,j)*EXPfactor*Solar1(j) source_down_TL(i) =source_down_TL(i) -trans_TL(i,j)*Solar1(j+nZ) -trans(i,j)*Solar1_TL(j+nZ) & -refl_TL(i,j)*EXPfactor*Solar1(j) -refl(i,j)*EXPfactor_TL*Solar1(j) -refl(i,j)*EXPfactor*Solar1_TL(j) END DO END DO ! ! specific treatment for downeward source function IF( abs( V0(N2,N2) ) > 0.0001_fp ) THEN source_down(nZ) =source_down(nZ) +(EXPfactor-trans(nZ,nZ))*Sfac2/V0(N2,N2) source_down_TL(nZ) =source_down_TL(nZ) +(EXPfactor_TL-trans_TL(nZ,nZ))*Sfac2/V0(N2,N2) & +(EXPfactor-trans(nZ,nZ))*Sfac2_TL/V0(N2,N2)-(EXPfactor-trans(nZ,nZ))*Sfac2*V0_TL(N2,N2)/V0(N2,N2)/V0(N2,N2) ELSE source_down(nZ) =source_down(nZ) -EXPfactor*Sfac2*optical_depth/COS_Angle(nZ) source_down_TL(nZ) =source_down_TL(nZ) -EXPfactor_TL*Sfac2*optical_depth/COS_Angle(nZ) & -EXPfactor*Sfac2_TL*optical_depth/COS_Angle(nZ)-EXPfactor*Sfac2*optical_depth_TL/COS_Angle(nZ) END IF ! source_up(1:nZ) = source_up(1:nZ)*s_transmittance source_up_TL(1:nZ) = source_up_TL(1:nZ)*s_transmittance+source_up(1:nZ)*s_transmittance_TL ! source_down(1:nZ) = source_down(1:nZ)*s_transmittance source_down_TL(1:nZ) = source_down_TL(1:nZ)*s_transmittance + source_down(1:nZ)*s_transmittance_TL END IF source_up_TL(:) = source_up_TL(:) + thermal_up_TL(:) source_down_TL(:) = source_down_TL(:) + thermal_down_TL(:) RETURN END SUBROUTINE CRTM_AMOM_layer_TL SUBROUTINE CRTM_ADA_AD(n_Layers, & ! Input number of atmospheric layers w, & ! Input layer scattering albedo T_OD, & ! Input layer optical depth cosmic_background, & ! Input cosmic background radiance emissivity, & ! Input surface emissivity direct_reflectivity, & ! surface direct reflectivity RTV, & ! Input structure containing forward results s_rad_up_AD, & ! Input adjoint upward radiance Planck_Atmosphere_AD, & ! Output AD atmospheric layer Planck radiance Planck_Surface_AD, & ! Output AD surface Planck radiance w_AD, & ! Output AD layer scattering albedo T_OD_AD, & ! Output AD layer optical depth emissivity_AD, & ! Output AD surface emissivity reflectivity_AD, & ! Output AD surface reflectivity direct_reflectivity_AD, & ! Output AD surface direct reflectivity Pff_AD, & ! Output AD forward phase matrix Pbb_AD) ! Output AD backward phase matrix ! ------------------------------------------------------------------------- ! ! FUNCTION: ! ! This subroutine calculates IR/MW adjoint radiance at the top of ! ! the atmosphere including atmospheric scattering. The structure RTV ! ! carried in forward part results. ! ! The CRTM_ADA_AD algorithm computes layer tangent-linear reflectance and ! ! transmittance as well as source function by the subroutine ! ! CRTM_Doubling_layer as source function by the subroutine ! ! CRTM_Doubling_layer, then uses ! ! an adding method to integrate the layer and surface components. ! ! ! ! Quanhua Liu Quanhua.Liu@noaa.gov ! ! ------------------------------------------------------------------------- ! IMPLICIT NONE INTEGER, INTENT(IN) :: n_Layers TYPE(RTV_type), INTENT(IN) :: RTV REAL (fp), INTENT(IN), DIMENSION( : ) :: w,T_OD REAL (fp), INTENT(IN), DIMENSION( : ) :: emissivity,direct_reflectivity REAL (fp), INTENT(IN) :: cosmic_background REAL (fp),INTENT(INOUT),DIMENSION( :,:,: ) :: Pff_AD, Pbb_AD REAL (fp),INTENT(INOUT),DIMENSION( : ) :: w_AD,T_OD_AD REAL (fp),INTENT(INOUT),DIMENSION( 0: ) :: Planck_Atmosphere_AD REAL (fp),INTENT(INOUT) :: Planck_Surface_AD REAL (fp),INTENT(INOUT),DIMENSION( : ) :: emissivity_AD,direct_reflectivity_AD REAL (fp),INTENT(INOUT),DIMENSION( :,: ) :: reflectivity_AD REAL (fp),INTENT(INOUT),DIMENSION( : ) :: s_rad_up_AD ! -------------- internal variables --------------------------------- ! ! Abbreviations: ! ! s: scattering, rad: radiance, trans: transmission, ! ! refl: reflection, up: upward, down: downward ! ! --------------------------------------------------------------------! REAL (fp), DIMENSION(RTV%n_Angles,RTV%n_Angles) :: temporal_matrix_AD REAL (fp), DIMENSION( RTV%n_Angles, n_Layers) :: refl_down REAL (fp), DIMENSION( RTV%n_Angles ) :: s_source_up_AD,s_source_down_AD,refl_down_AD REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles) :: s_trans_AD,s_refl_AD,Refl_Trans_AD REAL (fp), DIMENSION( RTV%n_Angles, RTV%n_Angles) :: s_refl_up_AD,Inv_Gamma_AD,Inv_GammaT_AD REAL (fp) :: sum_s_AD, sums_AD, xx REAL (fp), DIMENSION(0:n_Layers) :: total_opt, total_opt_AD INTEGER :: i, j, k,nZ ! nZ = RTV%n_Angles total_opt_AD = ZERO total_opt(0) = ZERO DO k = 1, n_Layers total_opt(k) = total_opt(k-1) + T_OD(k) END DO ! s_trans_AD = ZERO Planck_Atmosphere_AD = ZERO Planck_Surface_AD = ZERO s_refl_up_AD = ZERO Pff_AD = ZERO Pbb_AD = ZERO ! T_OD_AD = ZERO ! Adding reflected cosmic background radiation DO i = 1, RTV%n_Angles sum_s_AD = s_rad_up_AD(i)*cosmic_background DO j = 1, RTV%n_Angles s_refl_up_AD(i,j) = sum_s_AD ENDDO ENDDO ! DO 10 k = 1, n_Layers s_source_up_AD = ZERO s_source_down_AD = ZERO s_trans_AD = ZERO ! ! Compute tranmission and reflection matrices for a layer IF(w(k) > SCATTERING_ALBEDO_THRESHOLD) THEN refl_down(1:RTV%n_Angles,k) = matmul(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k), & RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k)) s_refl_AD = s_refl_up_AD Inv_GammaT_AD = matmul(s_refl_up_AD,transpose(RTV%Refl_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k))) Refl_Trans_AD = matmul(transpose(RTV%Inv_GammaT(1:RTV%n_Angles,1:RTV%n_Angles,k)),s_refl_up_AD) s_refl_up_AD=matmul(Refl_Trans_AD,transpose(RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k))) s_trans_AD=matmul(transpose(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k)),Refl_Trans_AD) s_source_up_AD(1:RTV%n_Angles) = s_rad_up_AD(1:RTV%n_Angles) DO i = 1, RTV%n_Angles sums_AD = s_rad_up_AD(i) DO j = 1, RTV%n_Angles Inv_GammaT_AD(i,j)=Inv_GammaT_AD(i,j)+sums_AD*(refl_down(j,k)+RTV%s_Level_Rad_UP(j,k)) ENDDO ENDDO 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)) 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)) DO i = 1, RTV%n_Angles sums_AD = refl_down_AD(i) DO j = 1, RTV%n_Angles s_refl_up_AD(i,j)=s_refl_up_AD(i,j)+sums_AD*RTV%s_Layer_Source_DOWN(j,k) ENDDO ENDDO s_source_down_AD=matmul(transpose(RTV%s_Level_Refl_UP(1:RTV%n_Angles,1:RTV%n_Angles,k)),refl_down_AD(:)) s_trans_AD=s_trans_AD+matmul(Inv_GammaT_AD,transpose(RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k))) Inv_Gamma_AD= matmul(transpose(RTV%s_Layer_Trans(1:RTV%n_Angles,1:RTV%n_Angles,k)),Inv_GammaT_AD) temporal_matrix_AD= -matmul(Inv_Gamma_AD,transpose(RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k))) temporal_matrix_AD=matmul(transpose(RTV%Inv_Gamma(1:RTV%n_Angles,1:RTV%n_Angles,k)),temporal_matrix_AD) 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))) 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) ! ----------------------------------------------------------- ! ! CALL Doubling algorithm to computing forward and tagent ! ! layer transmission, reflection, and source functions. ! ! ----------------------------------------------------------- ! call CRTM_AMOM_layer_AD(RTV%n_Streams,RTV%n_Angles,k,w(k),T_OD(k),total_opt(k-1) , & !Input RTV%COS_Angle,RTV%COS_Weight,RTV%Pff(:,:,k),RTV%Pbb(:,:,k),RTV%Planck_Atmosphere(k), & !Input s_trans_AD,s_refl_AD,s_source_up_AD,s_source_down_AD,RTV,w_AD(k),T_OD_AD(k) , & total_opt_AD(k-1),Pff_AD(:,:,k),Pbb_AD(:,:,k),Planck_Atmosphere_AD(k)) !Output ELSE DO i = 1, RTV%n_Angles DO j = 1, RTV%n_Angles 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) 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) 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) ENDDO ENDDO ! Adding method DO i = 1, RTV%n_Angles s_source_up_AD(i)=s_rad_up_AD(i) 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) & *RTV%s_Layer_Source_DOWN(1:RTV%n_Angles,k))+RTV%s_Level_Rad_UP(i,k)) sum_s_AD=RTV%s_Layer_Trans(i,i,k)*s_rad_up_AD(i) DO j = 1, RTV%n_Angles s_refl_up_AD(i,j)=s_refl_up_AD(i,j)+sum_s_AD*RTV%s_Layer_Source_DOWN(j,k) ENDDO sum_s_AD=RTV%s_Layer_Trans(i,i,k)*s_rad_up_AD(i) DO j = 1, RTV%n_Angles s_source_down_AD(j)=s_source_down_AD(j)+sum_s_AD*RTV%s_Level_Refl_UP(i,j,k) ENDDO s_rad_up_AD(i)=RTV%s_Layer_Trans(i,i,k)*s_rad_up_AD(i) ENDDO DO i = RTV%n_Angles, 1, -1 s_source_up_AD(i) = s_source_up_AD(i) + s_source_down_AD(i) s_trans_AD(i,i) = s_trans_AD(i,i) - RTV%Planck_Atmosphere(k) * s_source_up_AD(i) Planck_Atmosphere_AD(k) = Planck_Atmosphere_AD(k) + s_source_up_AD(i) * (ONE - RTV%s_Layer_Trans(i,i,k) ) s_source_up_AD(i) = ZERO T_OD_AD(k)=T_OD_AD(k)-s_trans_AD(i,i)/RTV%COS_Angle(i)*RTV%s_Layer_Trans(i,i,k) ENDDO ENDIF 10 CONTINUE ! IF( RTV%Solar_Flag_true ) THEN xx = RTV%Solar_irradiance/PI * exp(-total_opt(n_Layers)/RTV%COS_SUN) total_opt_AD(n_Layers) = total_opt_AD(n_Layers) & - xx*sum(direct_reflectivity(1:RTV%n_Angles)*s_rad_up_AD(1:RTV%n_Angles)) direct_reflectivity_AD = direct_reflectivity_AD + s_rad_up_AD * RTV%COS_SUN * xx END IF emissivity_AD = s_rad_up_AD * RTV%Planck_Surface Planck_Surface_AD = sum(emissivity(:) * s_rad_up_AD(:) ) reflectivity_AD = s_refl_up_AD ! s_rad_up_AD = ZERO s_refl_up_AD = ZERO DO k = n_Layers, 1, -1 T_OD_AD(k) = T_OD_AD(k) + total_opt_AD(k) total_opt_AD(k-1) = total_opt_AD(k-1) + total_opt_AD(k) END DO RETURN END SUBROUTINE CRTM_ADA_AD ! ! SUBROUTINE CRTM_AMOM_layer_AD( n_streams, & ! Input, number of streams nZ, & ! Input, number of angles KL, & ! Input, KL-th layer single_albedo, & ! Input, single scattering albedo optical_depth, & ! Input, layer optical depth total_opt, & ! Input, COS_Angle, & ! Input, COSINE of ANGLES COS_Weight, & ! Input, GAUSSIAN Weights ff, & ! Input, Phase matrix (forward part) bb, & ! Input, Phase matrix (backward part) Planck_Func, & ! Input, Planck for layer temperature trans_AD, & ! Input, layer tangent-linear trans refl_AD, & ! Input, layer tangent-linear refl source_up_AD, & ! Input, layer tangent-linear source_up source_down_AD, & ! Input, layer tangent-linear source_down RTV, & ! Input, structure containing forward results single_albedo_AD, & ! Output adjoint single scattering albedo optical_depth_AD, & ! Output AD layer optical depth total_opt_AD, & ! Output AD accumulated optical depth ftom TOA to current layer top ff_AD, & ! Output AD forward Phase matrix bb_AD, & ! Output AD backward Phase matrix Planck_Func_AD) ! Output AD Planck for layer temperature ! --------------------------------------------------------------------------------------- ! FUNCTION ! Compute AD layer transmission, reflection matrices and source function ! at the top and bottom of the layer. ! ! see also ADA method. ! Quanhua Liu ! Quanhua.Liu@noaa.gov ! ---------------------------------------------------------------------------------------- IMPLICIT NONE INTEGER, INTENT(IN) :: n_streams,nZ,KL TYPE(RTV_type), INTENT( IN ) :: RTV REAL(fp), INTENT(IN), DIMENSION(:,:) :: ff,bb REAL(fp), INTENT(IN), DIMENSION(:) :: COS_Angle, COS_Weight REAL(fp) :: single_albedo,optical_depth,Planck_Func,total_opt ! internal variables REAL(fp) :: s, c, s_AD,c_AD,xx_AD INTEGER :: i,j INTEGER :: Error_Status ! Tangent-Linear Part REAL(fp), INTENT( INOUT ), DIMENSION( :,: ) :: trans_AD,refl_AD REAL(fp), INTENT( INOUT ), DIMENSION( : ) :: source_up_AD,source_down_AD REAL(fp), INTENT( INOUT ) :: single_albedo_AD REAL(fp), INTENT( INOUT ) :: optical_depth_AD,Planck_Func_AD,total_opt_AD REAL(fp), INTENT(INOUT), DIMENSION(:,:) :: ff_AD,bb_AD REAL(fp), DIMENSION(nZ) :: Exp_x_AD,EigVa_AD,EigValue_AD REAL(fp), DIMENSION(nZ,nZ) :: i_Gm_AD, Gm_AD, Gp_AD, EigVe_AD, EigVeF_AD, EigVeVa_AD REAL(fp), DIMENSION(nZ,nZ) :: A1_AD,A2_AD,A3_AD,A4_AD,A5_AD,A6_AD,Gm_A5_AD,i_Gm_A5_AD REAL(fp), DIMENSION(nZ,nZ) :: HH_AD,PM_AD,PP_AD,PPM_AD,i_PPM_AD,PPP_AD REAL(fp), DIMENSION(nZ) :: thermal_up_AD, thermal_down_AD REAL(fp) :: EXPfactor,Sfactor,s_transmittance,Solar(2*nZ),V0(2*nZ,2*nZ),Solar1(2*nZ) REAL(fp) :: V1(2*nZ,2*nZ),Sfac2,source_up(nZ),source_down(nZ) REAL(fp) :: EXPfactor_AD,Sfactor_AD,s_transmittance_AD,Solar_AD(2*nZ),V0_AD(2*nZ,2*nZ),Solar1_AD(2*nZ) REAL(fp) :: V1_AD(2*nZ,2*nZ),Sfac2_AD REAL(fp), DIMENSION(nZ,nZ) :: trans, refl INTEGER :: N2, N2_1 REAL(fp) :: Thermal_C_AD CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_AMOM_layer_AD' CHARACTER(256) :: Message s_AD = ZERO c_AD = ZERO ! ! for small layer optical depth, single scattering is applied. IF( optical_depth < DELTA_OPTICAL_DEPTH ) THEN s = optical_depth * single_albedo DO i = 1, nZ c = s/COS_Angle(i) IF( RTV%mth_Azi == 0 ) THEN source_up_AD(i) = source_up_AD(i) + source_down_AD(i) source_down_AD(i) = ZERO Planck_Func_AD = Planck_Func_AD + (ONE - RTV%Thermal_C(i,KL))*source_up_AD(i) Thermal_C_AD = -source_up_AD(i) * Planck_Func END IF DO j = 1, nZ IF( RTV%mth_Azi == 0 ) THEN refl_AD(i,j) = refl_AD(i,j) + Thermal_C_AD trans_AD(i,j) = trans_AD(i,j) + Thermal_C_AD END IF IF( i == j ) THEN optical_depth_AD = optical_depth_AD - trans_AD(i,j)/COS_Angle(i) END IF c_AD = c_AD + trans_AD(i,j) * ff(i,j) * COS_Weight(j) ff_AD(i,j) = ff_AD(i,j) + c * trans_AD(i,j) * COS_Weight(j) c_AD = c_AD + refl_AD(i,j) * bb(i,j) * COS_Weight(j) bb_AD(i,j) = bb_AD(i,j) + c * refl_AD(i,j) * COS_Weight(j) ENDDO source_up_AD(i) = ZERO s_AD = s_AD + c_AD/COS_Angle(i) c_AD = ZERO ENDDO optical_depth_AD = optical_depth_AD + s_AD * single_albedo single_albedo_AD = single_albedo_AD + optical_depth * s_AD RETURN END IF trans(1:nZ,1:nZ) = RTV%s_Layer_Trans(1:nZ,1:nZ,KL) refl(1:nZ,1:nZ) = RTV%s_Layer_Refl(1:nZ,1:nZ,KL) thermal_up_AD(:) = source_up_AD(:) thermal_down_AD(:) = source_down_AD(:) IF( RTV%Solar_Flag_true ) THEN N2 = 2 * nZ N2_1 = N2 - 1 ! forward part start ******** source_up = ZERO source_down = ZERO Solar_AD = ZERO Solar1_AD = ZERO Sfactor_AD = ZERO Sfac2_AD = ZERO EXPfactor_AD = ZERO s_transmittance_AD = ZERO V0_AD = ZERO V1_AD = ZERO ! ! Solar source Sfactor = single_albedo*RTV%Solar_irradiance/PI IF( RTV%mth_Azi == 0 ) Sfactor = Sfactor/TWO EXPfactor = exp(-optical_depth/RTV%COS_SUN) s_transmittance = exp(-total_opt/RTV%COS_SUN) DO i = 1, nZ Solar(i) = -bb(i,nZ+1)*Sfactor Solar(i+nZ) = -ff(i,nZ+1)*Sfactor DO j = 1, nZ V0(i,j) = single_albedo * ff(i,j) * COS_Weight(j) V0(i+nZ,j) = single_albedo * bb(i,j) * COS_Weight(j) V0(i,j+nZ) = V0(i+nZ,j) V0(nZ+i,j+nZ) = V0(i,j) ENDDO V0(i,i) = V0(i,i) - ONE - COS_Angle(i)/RTV%COS_SUN V0(i+nZ,i+nZ) = V0(i+nZ,i+nZ) - ONE + COS_Angle(i)/RTV%COS_SUN ENDDO V1(1:N2_1,1:N2_1) = matinv(V0(1:N2_1,1:N2_1), Error_Status) IF( Error_Status /= SUCCESS ) THEN WRITE( Message,'("Error in matrix inversion matinv(V0(1:N2_1,1:N2_1), Error_Status) ")' ) CALL Display_Message( ROUTINE_NAME, & TRIM(Message), & Error_Status ) RETURN END IF Solar1(1:N2_1) = matmul( V1(1:N2_1,1:N2_1), Solar(1:N2_1) ) Solar1(N2) = ZERO Sfac2 = Solar(N2) - sum( V0(N2,1:N2_1)*Solar1(1:N2_1) ) DO i = 1, nZ source_up(i) = Solar1(i) source_down(i) = EXPfactor*Solar1(i+nZ) DO j = 1, nZ source_up(i) =source_up(i)-refl(i,j)*Solar1(j+nZ)-trans(i,j)*EXPfactor*Solar1(j) source_down(i) =source_down(i) -trans(i,j)*Solar1(j+nZ) -refl(i,j)*EXPfactor*Solar1(j) END DO END DO ! specific treatment for downeward source function IF( abs( V0(N2,N2) ) > 0.0001_fp ) THEN source_down(nZ) =source_down(nZ) +(EXPfactor-trans(nZ,nZ))*Sfac2/V0(N2,N2) ELSE source_down(nZ) =source_down(nZ) -EXPfactor*Sfac2*optical_depth/COS_Angle(nZ) END IF ! forward part end ******** ! s_transmittance_AD = s_transmittance_AD+sum (source_down(1:nZ)*source_down_AD(1:nZ) ) source_down_AD(1:nZ) = source_down_AD(1:nZ)*s_transmittance s_transmittance_AD = s_transmittance_AD + sum (source_up(1:nZ)*source_up_AD(1:nZ) ) source_up_AD(1:nZ) = source_up_AD(1:nZ)*s_transmittance ! ! specific treatment for downeward source function IF( abs( V0(N2,N2) ) > 0.0001_fp ) THEN V0_AD(N2,N2)=V0_AD(N2,N2)-(EXPfactor-trans(nZ,nZ))*Sfac2*source_down_AD(nZ)/V0(N2,N2)/V0(N2,N2) Sfac2_AD = Sfac2_AD+(EXPfactor-trans(nZ,nZ))*source_down_AD(nZ)/V0(N2,N2) EXPfactor_AD = EXPfactor_AD+source_down_AD(nZ)*Sfac2/V0(N2,N2) trans_AD(nZ,nZ) = trans_AD(nZ,nZ)-source_down_AD(nZ)*Sfac2/V0(N2,N2) ELSE optical_depth_AD = optical_depth_AD -EXPfactor*Sfac2*source_down_AD(nZ)/COS_Angle(nZ) Sfac2_AD = Sfac2_AD-EXPfactor*source_down_AD(nZ)*optical_depth/COS_Angle(nZ) EXPfactor_AD = EXPfactor_AD-source_down_AD(nZ)*Sfac2*optical_depth/COS_Angle(nZ) END IF DO i = nZ, 1, -1 DO j = nZ, 1, -1 Solar1_AD(j)=Solar1_AD(j)-refl(i,j)*EXPfactor*source_down_AD(i) EXPfactor_AD = EXPfactor_AD -refl(i,j)*source_down_AD(i)*Solar1(j) refl_AD(i,j) = refl_AD(i,j) -source_down_AD(i)*EXPfactor*Solar1(j) Solar1_AD(j+nZ) = Solar1_AD(j+nZ) -trans(i,j)*source_down_AD(i) trans_AD(i,j) = trans_AD(i,j) -source_down_AD(i)*Solar1(j+nZ) Solar1_AD(j)=Solar1_AD(j)-trans(i,j)*EXPfactor*source_up_AD(i) EXPfactor_AD = EXPfactor_AD - trans(i,j)*source_up_AD(i)*Solar1(j) trans_AD(i,j)=trans_AD(i,j) - source_up_AD(i)*EXPfactor*Solar1(j) Solar1_AD(j+nZ) = Solar1_AD(j+nZ) -refl(i,j)*source_up_AD(i) refl_AD(i,j) = refl_AD(i,j) -source_up_AD(i)*Solar1(j+nZ) END DO Solar1_AD(i+nZ) = Solar1_AD(i+nZ) + EXPfactor * source_down_AD(i) EXPfactor_AD = EXPfactor_AD + source_down_AD(i)*Solar1(i+nZ) Solar1_AD(i) = Solar1_AD(i) + source_up_AD(i) END DO Solar1_AD(1:N2_1)=Solar1_AD(1:N2_1) -Sfac2_AD*V0(N2,1:N2_1) V0_AD(N2,1:N2_1)=V0_AD(N2,1:N2_1) -Sfac2_AD*Solar1(1:N2_1) Solar_AD(N2) = Solar_AD(N2) + Sfac2_AD Solar1_AD(N2) = ZERO Solar_AD(1:N2_1)=Solar_AD(1:N2_1)+matmul( transpose(V1(1:N2_1,1:N2_1)),Solar1_AD(1:N2_1) ) Solar(1:N2_1) = matmul( V1(1:N2_1,1:N2_1),Solar(1:N2_1) ) Solar1_AD(1:N2_1) = -matmul( transpose(V1(1:N2_1,1:N2_1)),Solar1_AD(1:N2_1) ) DO i = 1, N2_1 DO j = 1, N2_1 V0_AD(i,j)=V0_AD(i,j)+Solar1_AD(i)*Solar(j) END DO END DO ! Solar source DO i = nZ, 1, -1 DO j = nZ, 1, -1 V0_AD(i,j)=V0_AD(i,j) + V0_AD(nZ+i,j+nZ) V0_AD(i+nZ,j)=V0_AD(i+nZ,j) + V0_AD(i,j+nZ) bb_AD(i,j)=bb_AD(i,j) + single_albedo*V0_AD(i+nZ,j)*COS_Weight(j) single_albedo_AD=single_albedo_AD + V0_AD(i+nZ,j)*bb(i,j)*COS_Weight(j) ff_AD(i,j)=ff_AD(i,j) + single_albedo*V0_AD(i,j)*COS_Weight(j) single_albedo_AD=single_albedo_AD +V0_AD(i,j)*ff(i,j)*COS_Weight(j) ENDDO Sfactor_AD = Sfactor_AD -ff(i,nZ+1)*Solar_AD(i+nZ) ff_AD(i,nZ+1) = ff_AD(i,nZ+1) -Solar_AD(i+nZ)*Sfactor Sfactor_AD = Sfactor_AD -bb(i,nZ+1)*Solar_AD(i) bb_AD(i,nZ+1)=bb_AD(i,nZ+1) - Solar_AD(i)*Sfactor ENDDO total_opt_AD = total_opt_AD -s_transmittance_AD/RTV%COS_SUN*s_transmittance optical_depth_AD = optical_depth_AD -EXPfactor_AD/RTV%COS_SUN*EXPfactor IF( RTV%mth_Azi == 0 ) THEN Sfactor_AD = Sfactor_AD/TWO END IF single_albedo_AD = single_albedo_AD + Sfactor_AD*RTV%Solar_irradiance/PI END IF ! Thermal part IF( RTV%mth_Azi == 0 ) THEN DO i = nZ, 1, -1 thermal_up_AD(i) = thermal_up_AD(i) + thermal_down_AD(i) thermal_down_AD(i) = ZERO Planck_Func_AD = Planck_Func_AD + ( ONE - RTV%Thermal_C(i,KL) ) * thermal_up_AD(i) Thermal_C_AD = -thermal_up_AD(i) * Planck_Func IF ( i == nZ .AND. nZ == (n_Streams+1) ) THEN trans_AD(nZ,nZ) = trans_AD(nZ,nZ) + Thermal_C_AD END IF DO j = n_Streams, 1, -1 trans_AD(i,j) = trans_AD(i,j) + Thermal_C_AD refl_AD(i,j) = refl_AD(i,j) + Thermal_C_AD ENDDO thermal_up_AD(i) = ZERO ENDDO END IF ! i_Gm_A5_AD = matmul( transpose(RTV%Gp(1:nZ,1:nZ,KL)-RTV%A6(1:nZ,1:nZ,KL)),refl_AD ) Gp_AD(:,:) = matmul( refl_AD, transpose(RTV%i_Gm_A5(1:nZ,1:nZ,KL)) ) A6_AD = - GP_AD 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 ) A4_AD(:,:) = matmul( trans_AD, transpose(RTV%i_Gm_A5(1:nZ,1:nZ,KL)) ) A3_AD = - A4_AD 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)) ) ) Gm_AD = Gm_A5_AD A5_AD = - Gm_A5_AD A4_AD = A4_AD + matmul( A6_AD(:,:), transpose(RTV%A2(1:nZ,1:nZ,KL)) ) A2_AD = matmul( transpose(RTV%A4(1:nZ,1:nZ,KL)),A6_AD(:,:) ) A1_AD = matmul( A5_AD(:,:), transpose(RTV%A2(1:nZ,1:nZ,KL)) ) A2_AD = A2_AD + matmul( transpose(RTV%A1(1:nZ,1:nZ,KL)), A5_AD(:,:) ) Gp_AD = Gp_AD + matmul( A3_AD(:,:), transpose(RTV%A2(1:nZ,1:nZ,KL)) ) A2_AD = A2_AD + matmul( transpose(RTV%Gp(1:nZ,1:nZ,KL)),A3_AD(:,:) ) i_Gm_AD = matmul( A2_AD(:,:), transpose(RTV%A1(1:nZ,1:nZ,KL)) ) A1_AD = A1_AD + matmul( transpose(RTV%i_Gm(1:nZ,1:nZ,KL)), A2_AD(:,:) ) Exp_x_AD = ZERO DO i = nZ, 1, -1 DO j = nZ, 1, -1 Gm_AD(i,j) = Gm_AD(i,j) + A4_AD(i,j)* RTV%Exp_x(j,KL) Exp_x_AD(j) = Exp_x_AD(j) + RTV%Gm(i,j,KL)*A4_AD(i,j) Gp_AD(i,j) = Gp_AD(i,j) + A1_AD(i,j)* RTV%Exp_x(j,KL) Exp_x_AD(j) = Exp_x_AD(j) + RTV%Gp(i,j,KL)*A1_AD(i,j) END DO END DO DO i = nZ, 1, -1 xx_AD = -Exp_x_AD(i)*RTV%Exp_x(i,KL) Exp_x_AD(i) = ZERO EigValue_AD(i) = xx_AD*optical_depth optical_depth_AD = optical_depth_AD + RTV%EigValue(i,KL)*xx_AD END DO 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)) ) ) EigVe_AD(:,:) = Gm_AD(:,:)/2.0_fp EigVeF_AD(:,:) = - Gm_AD(:,:)/2.0_fp EigVe_AD = EigVe_AD + Gp_AD(:,:)/2.0_fp EigVeF_AD = EigVeF_AD + Gp_AD(:,:)/2.0_fp i_PPM_AD(:,:) = matmul( EigVeF_AD(:,:), transpose(RTV%EigVeVa(1:nZ,1:nZ,KL)) ) EigVeVa_AD(:,:) = matmul( transpose(RTV%i_PPM(1:nZ,1:nZ,KL)), EigVeF_AD(:,:) ) DO i = nZ, 1, -1 DO j = nZ, 1, -1 EigVe_AD(i,j)=EigVe_AD(i,j)+EigVeVa_AD(i,j)* RTV%EigValue(j,KL) EigValue_AD(j) = EigValue_AD(j)+RTV%EigVe(i,j,KL)*EigVeVa_AD(i,j) END DO END DO DO i = nZ, 1, -1 IF( RTV%EigVa(i,KL) > ZERO ) THEN EigVa_AD(i) = 0.5_fp*EigValue_AD(i)/RTV%EigValue(i,KL) ELSE EigValue_AD(i) = ZERO EigVa_AD(i) = ZERO END IF END DO ! compute eigenvectors EigVe, and eigenvalues EigVa CALL ASYMTX_AD(nZ,RTV%EigVe(1:nZ,1:nZ,KL),RTV%EigVa(1:nZ,KL), & EigVe_AD,EigVa_AD,HH_AD,Error_Status) PPM_AD(:,:) = matmul( HH_AD(:,:), transpose(RTV%PPP(1:nZ,1:nZ,KL)) ) PPP_AD(:,:) = matmul( transpose(RTV%PPM(1:nZ,1:nZ,KL)), HH_AD(:,:) ) PP_AD = PPP_AD PM_AD = PPP_AD 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))) ) PP_AD = PP_AD + PPM_AD PM_AD = PM_AD - PPM_AD IF( single_albedo < max_albedo ) THEN s = single_albedo ELSE s = max_albedo END IF c_AD = ZERO s_AD = ZERO DO i = nZ, 1, -1 c = s/COS_Angle(i) DO j = nZ, 1, -1 c_AD = c_AD + PP_AD(i,j) * ff(i,j) * COS_Weight(j) ff_AD(i,j) = ff_AD(i,j) + c * PP_AD(i,j) * COS_Weight(j) c_AD = c_AD + PM_AD(i,j) * bb(i,j) * COS_Weight(j) bb_AD(i,j) = bb_AD(i,j) + c * PM_AD(i,j) * COS_Weight(j) END DO s_AD = s_AD + c_AD/COS_Angle(i) c_AD = ZERO ENDDO ! IF( single_albedo < max_albedo ) THEN s = single_albedo single_albedo_AD = s_AD + single_albedo_AD ELSE s = max_albedo s_AD = 0.0_fp END IF ! RETURN END SUBROUTINE CRTM_AMOM_layer_AD END MODULE ADA_Module