! ! CRTM_CloudScatter ! ! Module to compute the cloud particle absorption and scattering properties ! required for radiative transfer in a cloudy atmosphere. ! ! ! CREATION HISTORY ! Written by: Quanhua Liu, QSS Group, Inc; Quanhua.Liu@noaa.gov ! Yong Han, NOAA/NESDIS; Yong.Han@noaa.gov ! Paul van Delst, CIMSS/SSEC; paul.vandelst@ssec.wisc.edu ! 02-July-2005 ! MODULE CRTM_CloudScatter ! ----------------- ! Environment setup ! ----------------- ! Module use USE Type_Kinds, ONLY: fp USE Message_Handler, ONLY: SUCCESS, FAILURE, Display_Message USE CRTM_Parameters, ONLY: ZERO, ONE, POINT_5, ONEpointFIVE, & MAX_N_LAYERS, & MAX_N_CLOUDS, & WATER_CONTENT_THRESHOLD, & BS_THRESHOLD, & MAX_N_LEGENDRE_TERMS, & MAX_N_PHASE_ELEMENTS, & HGPHASE ! <<< NEED TO REMOVE THIS IN FUTURE USE CRTM_SpcCoeff, ONLY: SC, & INFRARED_SENSOR, MICROWAVE_SENSOR, VISIBLE_SENSOR USE CRTM_CloudCoeff, ONLY: CloudC USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type, & WATER_CLOUD, & ICE_CLOUD, & RAIN_CLOUD, & SNOW_CLOUD, & GRAUPEL_CLOUD, & HAIL_CLOUD USE CRTM_GeometryInfo_Define, ONLY: CRTM_GeometryInfo_type USE CRTM_Interpolation, ONLY: NPTS , & LPoly_type , & find_index , & interp_1D , & interp_2D , & interp_3D , & interp_2D_TL, & interp_3D_TL, & interp_2D_AD, & interp_3D_AD, & Clear_LPoly , & LPoly , & LPoly_TL , & LPoly_AD USE CRTM_AtmScatter_Define, ONLY: CRTM_AtmScatter_type , & CRTM_Associated_AtmScatter, & CRTM_Destroy_AtmScatter , & CRTM_Allocate_AtmScatter , & CRTM_Assign_AtmScatter , & CRTM_Zero_AtmScatter ! Internal variable definition module USE CSvar_Define, ONLY: CSvar_type, & CSinterp_type, & CSvar_Associated, & CSvar_Destroy , & CSvar_Create ! Disable implicit typing IMPLICIT NONE ! ------------ ! Visibilities ! ------------ ! Everything private by default PRIVATE ! Procedures PUBLIC :: CRTM_Compute_CloudScatter PUBLIC :: CRTM_Compute_CloudScatter_TL PUBLIC :: CRTM_Compute_CloudScatter_AD PUBLIC :: CRTM_Associated_AtmScatter PUBLIC :: CRTM_Destroy_AtmScatter PUBLIC :: CRTM_Allocate_AtmScatter PUBLIC :: CRTM_Assign_AtmScatter PUBLIC :: CRTM_Zero_AtmScatter ! ! FOR TESTING ONLY ! PUBLIC :: Get_Cloud_Opt_IR ! PUBLIC :: Get_Cloud_Opt_IR_TL ! PUBLIC :: Get_Cloud_Opt_IR_AD ! PUBLIC :: Get_Cloud_Opt_MW ! PUBLIC :: Get_Cloud_Opt_MW_TL ! PUBLIC :: Get_Cloud_Opt_MW_AD ! ----------------- ! Module parameters ! ----------------- ! RCS Id for the module CHARACTER(*), PARAMETER :: MODULE_RCS_ID = & '$Id: CRTM_CloudScatter.f90 22707 2012-11-21 21:09:10Z paul.vandelst@noaa.gov $' ! Message string length INTEGER, PARAMETER :: ML = 256 ! Number of stream angle definitions INTEGER, PARAMETER :: TWO_STREAMS = 2 INTEGER, PARAMETER :: FOUR_STREAMS = 4 INTEGER, PARAMETER :: SIX_STREAMS = 6 INTEGER, PARAMETER :: EIGHT_STREAMS = 8 INTEGER, PARAMETER :: SIXTEEN_STREAMS = 16 INTEGER, PARAMETER :: THIRTYTWO_STREAMS = 32 CONTAINS !------------------------------------------------------------------------------ !:sdoc+: ! ! NAME: ! CRTM_Compute_CloudScatter ! ! PURPOSE: ! Function to compute the cloud particle absorption and scattering ! properties and populate the output CloudScatter structure for a ! single channel. ! ! CALLING SEQUENCE: ! Error_Status = CRTM_Compute_CloudScatter( Atmosphere , & ! SensorIndex , & ! ChannelIndex, & ! CloudScatter, & ! CSvar ) ! ! INPUT ARGUMENTS: ! Atmosphere: CRTM_Atmosphere structure containing the atmospheric ! profile data. ! UNITS: N/A ! TYPE: TYPE(CRTM_Atmosphere_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! SensorIndex: Sensor index id. This is a unique index associated ! with a (supported) sensor used to access the ! shared coefficient data for a particular sensor. ! See the ChannelIndex argument. ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! ChannelIndex: Channel index id. This is a unique index associated ! with a (supported) sensor channel used to access the ! shared coefficient data for a particular sensor's ! channel. ! See the SensorIndex argument. ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! OUTPUT ARGUMENTS: ! CloudScatter: CRTM_AtmScatter structure containing the cloud particle ! absorption and scattering properties required for ! radiative transfer. ! UNITS: N/A ! TYPE: TYPE(CRTM_AtmScatter_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN OUT) ! ! CSvar: Structure containing internal variables required for ! subsequent tangent-linear or adjoint model calls. ! UNITS: N/A ! TYPE: CSvar_type ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(OUT) ! ! FUNCTION RESULT: ! Error_Status: The return value is an integer defining the error status. ! The error codes are defined in the Message_Handler module. ! If == SUCCESS the computation was sucessful ! == FAILURE an unrecoverable error occurred ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! !:sdoc-: !------------------------------------------------------------------------------ FUNCTION CRTM_Compute_CloudScatter( & Atm , & ! Input SensorIndex , & ! Input ChannelIndex, & ! Input CScat , & ! Output CSV ) & ! Internal variable output RESULT( Error_Status ) ! Arguments TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm INTEGER , INTENT(IN) :: SensorIndex INTEGER , INTENT(IN) :: ChannelIndex TYPE(CRTM_AtmScatter_type), INTENT(IN OUT) :: CScat TYPE(CSvar_type) , INTENT(IN OUT) :: CSV ! Function result INTEGER :: Error_Status ! Function parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Compute_CloudScatter' ! Local variables CHARACTER(ML) :: Message INTEGER :: k, kc, l, m, n INTEGER :: Sensor_Type REAL(fp) :: Frequency_MW, Frequency_IR LOGICAL :: Layer_Mask(Atm%n_Layers) INTEGER :: Layer_Index(Atm%n_Layers) INTEGER :: nCloud_Layers REAL(fp) :: bs ! ------ ! Set up ! ------ Error_Status = SUCCESS IF (Atm%n_Clouds == 0) RETURN CSV%Total_bs = ZERO ! Spectral variables Sensor_Type = SC(SensorIndex)%Sensor_Type Frequency_MW = SC(SensorIndex)%Frequency(ChannelIndex) Frequency_IR = SC(SensorIndex)%Wavenumber(ChannelIndex) ! Determine offset for Legendre coefficients in ! the CloudC lookup table corresponding to the ! number of streams SELECT CASE(CScat%n_Legendre_Terms) CASE (TWO_STREAMS) ; CScat%lOffset = 0 CASE (FOUR_STREAMS) ; CScat%lOffset = 0 CASE (SIX_STREAMS) ; CScat%lOffset = 5 CASE (EIGHT_STREAMS) ; CScat%lOffset = 12 CASE (SIXTEEN_STREAMS); CScat%lOffset = 21 CASE DEFAULT CScat%lOffset = 0 ! Is this correct? ! Use two-stream model or HG and RAYLEIGH Phase function IF( HGPHASE ) THEN CScat%n_Legendre_Terms = 0 ELSE Error_Status = FAILURE WRITE(Message,'("The n_Legendre_Terms in CloudScatter, ",i0,", do not fit model")') & CScat%n_Legendre_Terms CALL Display_Message( ROUTINE_NAME,Message,Error_Status ) RETURN END IF END SELECT ! --------------------------------------------- ! Loop over the different clouds in the profile ! --------------------------------------------- Cloud_loop: DO n = 1, Atm%n_Clouds ! Only process clouds with more ! than the threshold water amount Layer_Mask = Atm%Cloud(n)%Water_Content > WATER_CONTENT_THRESHOLD nCloud_Layers = COUNT(Layer_Mask) IF ( nCloud_Layers == 0 ) CYCLE Cloud_loop ! ------------------------------------ ! Loop over the current cloud's layers ! ------------------------------------ Layer_Index(1:nCloud_Layers) = PACK((/(k, k=1,Atm%Cloud(n)%n_Layers)/), Layer_Mask) Cloud_Layer_loop: DO k = 1, nCloud_Layers kc = Layer_Index(k) ! Call sensor specific routines SELECT CASE (Sensor_Type) CASE (MICROWAVE_SENSOR) CALL Get_Cloud_Opt_MW(CScat , & ! Input Frequency_MW , & ! Input Atm%Cloud(n)%Type , & ! Input Atm%Cloud(n)%Effective_Radius(kc), & ! Input Atm%Temperature(kc) , & ! Input CSV%ke(kc,n) , & ! Output CSV%w(kc,n) , & ! Output CSV%pcoeff(:,:,kc,n) , & ! Output CSV%csi(kc,n) ) ! Interpolation ! IR and visible use the same cloud optical data file, but distingished with Frequency CASE (INFRARED_SENSOR, VISIBLE_SENSOR) CALL Get_Cloud_Opt_IR(CScat , & ! Input Frequency_IR , & ! Input Atm%Cloud(n)%Type , & ! Input Atm%Cloud(n)%Effective_Radius(kc), & ! Input CSV%ke(kc,n) , & ! Output CSV%w(kc,n) , & ! Output CSV%pcoeff(:,:,kc,n) , & ! Output CSV%csi(kc,n) ) ! Interpolation CASE DEFAULT CSV%ke(kc,n) = ZERO CSV%w(kc,n) = ZERO CSV%pcoeff(:,:,kc,n) = ZERO END SELECT ! interpolation quality control IF( CSV%ke(kc,n) <= ZERO ) THEN CSV%ke(kc,n) = ZERO CSV%w(kc,n) = ZERO END IF IF( CSV%w(kc,n) <= ZERO ) THEN CSV%w(kc,n) = ZERO CSV%pcoeff(:,:,kc,n) = ZERO END IF IF( CSV%w(kc,n) >= ONE ) THEN CSV%w(kc,n) = ONE END IF ! Compute the volume scattering coefficient for the current ! cloud layer and accumulate it for the layer total for the ! profile (i.e. all clouds) ! bs = rho.w.ke ! where ! bs = volume scattering coefficient for a layer [dimensionless] ! rho = integrated cloud water density for a layer (g/m^2) [M.L^-2] ! w = single scatter albedo [dimensionless] ! ke = mass extintion coefficient (m^2/g) [L^2.M^-1] bs = Atm%Cloud(n)%Water_Content(kc) * CSV%ke(kc,n) * CSV%w(kc,n) CSV%Total_bs(kc) = CSV%Total_bs(kc) + bs ! Compute the optical depth (absorption + scattering) ! tau = rho.ke ! where ! rho = integrated cloud water density for a layer (g/m^2) [M.L^-2] ! ke = mass extintion coefficient (m^2/g) [L^2.M^-1] ! Note that since all these computations are done for a given ! layer, the optical depth is the same as the volume extinction ! coefficient, be. Usually, ! tau = be.d(z) ! but we are working with height/thickness independent quantities ! so that ! tau = be ! This is why the optical depth is used in the denominator to ! compute the single scatter albedo in the Layer_loop below. CScat%Optical_Depth(kc) = CScat%Optical_Depth(kc) + & (CSV%ke(kc,n)*Atm%Cloud(n)%Water_Content(kc)) CScat%Single_Scatter_Albedo(kc) = CScat%Single_Scatter_Albedo(kc) + bs ! Compute the phase matrix coefficients ! p = p + p(LUT)*bs ! where ! p(LUT) = the phase coefficient from the LUT IF( CScat%n_Phase_Elements > 0 ) THEN DO m = 1, CScat%n_Phase_Elements DO l = 1, CScat%n_Legendre_Terms CScat%Phase_Coefficient(l,m,kc) = CScat%Phase_Coefficient(l,m,kc) + & (CSV%pcoeff(l,m,kc,n) * bs) END DO END DO END IF END DO Cloud_Layer_loop END DO Cloud_loop END FUNCTION CRTM_Compute_CloudScatter !------------------------------------------------------------------------------ !:sdoc+: ! ! NAME: ! CRTM_Compute_CloudScatter_TL ! ! PURPOSE: ! Function to compute the tangent-linear cloud particle absorption and ! scattering properties and populate the output CloudScatter_TL structure ! for a single channel. ! ! CALLING SEQUENCE: ! Error_Status = CRTM_Compute_CloudScatter_TL( Atmosphere , & ! CloudScatter , & ! Atmosphere_TL , & ! SensorIndex , & ! ChannelIndex , & ! CloudScatter_TL, & ! CSvar ) ! ! INPUT ARGUMENTS: ! Atmosphere: CRTM_Atmosphere structure containing the atmospheric ! profile data. ! UNITS: N/A ! TYPE: TYPE(CRTM_Atmosphere_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! CloudScatter: CRTM_AtmScatter structure containing the forward model ! cloud particle absorption and scattering properties ! required for radiative transfer. ! UNITS: N/A ! TYPE: TYPE(CRTM_AtmScatter_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! Atmosphere_TL: CRTM Atmosphere structure containing the tangent-linear ! atmospheric state data. ! UNITS: N/A ! TYPE: TYPE(CRTM_Atmosphere_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! SensorIndex: Sensor index id. This is a unique index associated ! with a (supported) sensor used to access the ! shared coefficient data for a particular sensor. ! See the ChannelIndex argument. ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! ChannelIndex: Channel index id. This is a unique index associated ! with a (supported) sensor channel used to access the ! shared coefficient data for a particular sensor's ! channel. ! See the SensorIndex argument. ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! CSvar: Structure containing internal variables required for ! subsequent tangent-linear or adjoint model calls. ! UNITS: N/A ! TYPE: CSvar_type ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! OUTPUT ARGUMENTS: ! CloudScatter_TL: CRTM_AtmScatter structure containing the tangent-linear ! cloud particle absorption and scattering properties ! required for radiative transfer. ! UNITS: N/A ! TYPE: TYPE(CRTM_AtmScatter_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN OUT) ! ! ! FUNCTION RESULT: ! Error_Status: The return value is an integer defining the error status. ! The error codes are defined in the Message_Handler module. ! If == SUCCESS the computation was sucessful ! == FAILURE an unrecoverable error occurred ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! !:sdoc-: !------------------------------------------------------------------------------ FUNCTION CRTM_Compute_CloudScatter_TL( & Atm , & ! FWD Input CScat , & ! FWD Input Atm_TL , & ! TL Input SensorIndex , & ! Input ChannelIndex, & ! Input CScat_TL , & ! TL Output CSV ) & ! Internal variable input RESULT( Error_Status ) ! Arguments TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm TYPE(CRTM_AtmScatter_type), INTENT(IN) :: CScat TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm_TL INTEGER , INTENT(IN) :: SensorIndex INTEGER , INTENT(IN) :: ChannelIndex TYPE(CRTM_AtmScatter_type), INTENT(IN OUT) :: CScat_TL TYPE(CSvar_type) , INTENT(IN) :: CSV ! Function result INTEGER :: Error_Status ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Compute_CloudScatter' ! Local variables INTEGER :: k, kc, l, m, n INTEGER :: n_Legendre_Terms, n_Phase_Elements INTEGER :: Sensor_Type REAL(fp) :: Frequency_MW, Frequency_IR LOGICAL :: Layer_Mask(Atm%n_Layers) INTEGER :: Layer_Index(Atm%n_Layers) INTEGER :: nCloud_Layers REAL(fp) :: ke_TL, w_TL REAL(fp) :: pcoeff_TL(0:CScat%n_Legendre_Terms, CScat%n_Phase_Elements) REAL(fp) :: bs, bs_TL ! ------ ! Set up ! ------ Error_Status = SUCCESS IF (Atm%n_Clouds == 0) RETURN ! Spectral variables Sensor_Type = SC(SensorIndex)%Sensor_Type Frequency_MW = SC(SensorIndex)%Frequency(ChannelIndex) Frequency_IR = SC(SensorIndex)%Wavenumber(ChannelIndex) ! Phase matrix dimensions n_Legendre_Terms = CScat_TL%n_Legendre_Terms n_Phase_Elements = CScat_TL%n_Phase_Elements CScat_TL%lOffset = CScat%lOffset ! --------------------------------------------- ! Loop over the different clouds in the profile ! --------------------------------------------- Cloud_loop: DO n = 1, Atm%n_Clouds ! Only process clouds with more ! than the threshold water amount Layer_Mask = Atm%Cloud(n)%Water_Content > WATER_CONTENT_THRESHOLD nCloud_Layers = COUNT(Layer_Mask) IF ( nCloud_Layers == 0 ) CYCLE Cloud_loop ! ------------------------------------ ! Loop over the current cloud's layers ! ------------------------------------ Layer_Index(1:nCloud_Layers) = PACK((/(k, k=1,Atm%Cloud(n)%n_Layers)/), Layer_Mask) Cloud_Layer_loop: DO k = 1, nCloud_Layers kc = Layer_Index(k) ! Call sensor specific routines SELECT CASE (Sensor_Type) CASE (MICROWAVE_SENSOR) CALL Get_Cloud_Opt_MW_TL(CScat_TL , & ! Input Atm%Cloud(n)%Type , & ! Input Atm_TL%Cloud(n)%Effective_Radius(kc), & ! TL Input Atm_TL%Temperature(kc) , & ! TL Input ke_TL , & ! TL Output w_TL , & ! TL Output pcoeff_TL , & ! TL Output CSV%csi(kc,n) ) ! Interpolation CASE (INFRARED_SENSOR, VISIBLE_SENSOR) CALL Get_Cloud_Opt_IR_TL(CScat_TL , & ! Input Atm%Cloud(n)%Type , & ! Input Atm_TL%Cloud(n)%Effective_Radius(kc), & ! TL Input ke_TL , & ! TL Output w_TL , & ! TL Output pcoeff_TL , & ! TL Output CSV%csi(kc,n) ) ! Interpolation CASE DEFAULT ke_TL = ZERO w_TL = ZERO pcoeff_TL = ZERO END SELECT ! interpolation quality control IF( CSV%ke(kc,n) <= ZERO ) THEN ke_TL = ZERO w_TL = ZERO END IF IF( CSV%w(kc,n) <= ZERO ) THEN w_TL = ZERO pcoeff_TL = ZERO END IF IF( CSV%w(kc,n) >= ONE ) THEN w_TL = ZERO END IF ! Compute the volume scattering coefficient bs = Atm%Cloud(n)%Water_Content(kc) * CSV%ke(kc,n) * CSV%w(kc,n) bs_TL = (Atm_TL%Cloud(n)%Water_Content(kc) * CSV%ke(kc,n) * CSV%w(kc,n) ) + & (Atm%Cloud(n)%Water_Content(kc) * ke_TL * CSV%w(kc,n) ) + & (Atm%Cloud(n)%Water_Content(kc) * CSV%ke(kc,n) * w_TL ) CScat_TL%Single_Scatter_Albedo(kc) = CScat_TL%Single_Scatter_Albedo(kc) + bs_TL ! Compute the optical depth (absorption + scattering) CScat_TL%Optical_Depth(kc) = CScat_TL%Optical_Depth(kc) + & (ke_TL * Atm%Cloud(n)%Water_Content(kc)) + & (CSV%ke(kc,n) * Atm_TL%Cloud(n)%Water_Content(kc)) ! Compute the phase matrix coefficients IF( n_Phase_Elements > 0 ) THEN DO m = 1, n_Phase_Elements DO l = 0, n_Legendre_Terms CScat_TL%Phase_Coefficient(l,m,kc) = CScat_TL%Phase_Coefficient(l,m,kc) + & (pcoeff_TL(l,m) * bs ) + & (CSV%pcoeff(l,m,kc,n) * bs_TL) END DO END DO END IF END DO Cloud_Layer_loop END DO Cloud_loop END FUNCTION CRTM_Compute_CloudScatter_TL !------------------------------------------------------------------------------ !:sdoc+: ! ! NAME: ! CRTM_Compute_CloudScatter_AD ! ! PURPOSE: ! Function to compute the adjoint of the cloud particle absorption and ! scattering properties for a single channel. ! ! CALLING SEQUENCE: ! Error_Status = CRTM_Compute_CloudScatter_AD( Atmosphere , & ! CloudScatter , & ! CloudScatter_AD, & ! SensorIndex , & ! ChannelIndex , & ! Atmosphere_AD , & ! CSvar ) ! ! INPUT ARGUMENTS: ! Atmosphere: CRTM_Atmosphere structure containing the atmospheric ! profile data. ! UNITS: N/A ! TYPE: TYPE(CRTM_Atmosphere_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! CloudScatter: CRTM_AtmScatter structure containing the forward model ! cloud particle absorption and scattering properties ! required for radiative transfer. ! UNITS: N/A ! TYPE: TYPE(CRTM_AtmScatter_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! CloudScatter_AD: CRTM_AtmScatter structure containing the adjoint ! of the cloud particle absorption and scattering ! properties required for radiative transfer. ! **NOTE: On EXIT from this function, the contents of ! this structure may be modified (e.g. set to ! zero.) ! UNITS: N/A ! TYPE: TYPE(CRTM_AtmScatter_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN OUT) ! ! SensorIndex: Sensor index id. This is a unique index associated ! with a (supported) sensor used to access the ! shared coefficient data for a particular sensor. ! See the ChannelIndex argument. ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! ChannelIndex: Channel index id. This is a unique index associated ! with a (supported) sensor channel used to access the ! shared coefficient data for a particular sensor's ! channel. ! See the SensorIndex argument. ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! CSvar: Structure containing internal variables required for ! subsequent tangent-linear or adjoint model calls. ! UNITS: N/A ! TYPE: CSvar_type ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! OUTPUT ARGUMENTS: ! Atmosphere_AD: CRTM Atmosphere structure containing the adjoint ! atmospheric state data. ! UNITS: N/A ! TYPE: TYPE(CRTM_Atmosphere_type) ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN OUT) ! ! ! FUNCTION RESULT: ! Error_Status: The return value is an integer defining the error status. ! The error codes are defined in the Message_Handler module. ! If == SUCCESS the computation was sucessful ! == FAILURE an unrecoverable error occurred ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! !:sdoc-: !------------------------------------------------------------------------------ FUNCTION CRTM_Compute_CloudScatter_AD( & Atm , & ! FWD Input CScat , & ! FWD Input CScat_AD , & ! AD Input SensorIndex , & ! Input ChannelIndex, & ! Input Atm_AD , & ! AD Output CSV ) & ! Internal variable input RESULT( Error_Status ) ! Arguments TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm TYPE(CRTM_AtmScatter_type), INTENT(IN) :: CScat TYPE(CRTM_AtmScatter_type), INTENT(IN OUT) :: CScat_AD INTEGER , INTENT(IN) :: SensorIndex INTEGER , INTENT(IN) :: ChannelIndex TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_AD TYPE(CSvar_type) , INTENT(IN) :: CSV ! Function result INTEGER :: Error_Status ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Compute_CloudScatter_AD' ! Local variables INTEGER :: k, kc, l, m, n INTEGER :: n_Legendre_Terms, n_Phase_Elements INTEGER :: Sensor_Type REAL(fp) :: Frequency_MW, Frequency_IR LOGICAL :: Layer_Mask(Atm%n_Layers) INTEGER :: Layer_Index(Atm%n_Layers) INTEGER :: nCloud_Layers REAL(fp) :: ke_AD, w_AD REAL(fp) :: pcoeff_AD(0:CScat%n_Legendre_Terms, CScat%n_Phase_Elements) REAL(fp) :: bs, bs_AD ! ------ ! Set up ! ------ Error_Status = SUCCESS IF (Atm%n_Clouds == 0) RETURN ! Initialize local adjoint variables ! Spectral variables Sensor_Type = SC(SensorIndex)%Sensor_Type Frequency_MW = SC(SensorIndex)%Frequency(ChannelIndex) Frequency_IR = SC(SensorIndex)%Wavenumber(ChannelIndex) ! Phase matrix dimensions n_Legendre_Terms = CScat_AD%n_Legendre_Terms n_Phase_Elements = CScat_AD%n_Phase_Elements CScat_AD%lOffset = CScat%lOffset ! --------------------------------------------- ! Loop over the different clouds in the profile ! --------------------------------------------- Cloud_loop: DO n = 1, Atm%n_Clouds ! Only process clouds with more ! than the threshold water amount Layer_Mask = Atm%Cloud(n)%Water_Content > WATER_CONTENT_THRESHOLD nCloud_Layers = COUNT(Layer_Mask) IF ( nCloud_Layers == 0 ) CYCLE Cloud_loop ! ------------------------------------ ! Loop over the current cloud's layers ! ------------------------------------ Layer_Index(1:nCloud_Layers) = PACK((/(k, k=1,Atm%Cloud(n)%n_Layers)/), Layer_Mask) Cloud_Layer_loop: DO k = nCloud_Layers, 1, -1 !1, nCloud_Layers kc = Layer_Index(k) ! Initialise the individual ! cloud adjoint variables bs_AD = ZERO pcoeff_AD = ZERO ke_AD = ZERO w_AD = ZERO ! Recompute the forward model volume scattering ! coefficient for the current cloud ONLY bs = Atm%Cloud(n)%Water_Content(kc) * CSV%ke(kc,n) * CSV%w(kc,n) ! Compute the adjoint of the ! phase matrix coefficients IF( n_Phase_Elements > 0 ) THEN DO m = 1, n_Phase_Elements DO l = 0, n_Legendre_Terms bs_AD = bs_AD + (CSV%pcoeff(l,m,kc,n) * CScat_AD%Phase_Coefficient(l,m,kc)) pcoeff_AD(l,m) = pcoeff_AD(l,m) + (bs * CScat_AD%Phase_Coefficient(l,m,kc)) END DO END DO END IF ! Compute the adjoint of the optical ! depth (absorption + scattering) Atm_AD%Cloud(n)%Water_Content(kc) = Atm_AD%Cloud(n)%Water_Content(kc) + & (CSV%ke(kc,n) * CScat_AD%Optical_Depth(kc)) ke_AD = ke_AD + (Atm%Cloud(n)%Water_Content(kc) * CScat_AD%Optical_Depth(kc)) ! Compute the adjoint of the volume ! scattering coefficient. ! NOTE: bs_AD is not reinitialised after this ! point since it is reinitialised at the ! start of the Cloud_Layer_loop. bs_AD = bs_AD + CScat_AD%Single_Scatter_Albedo(kc) w_AD = w_AD + (Atm%Cloud(n)%Water_Content(kc) * CSV%ke(kc,n)* bs_AD ) ke_AD = ke_AD + (Atm%Cloud(n)%Water_Content(kc) * bs_AD * CSV%w(kc,n) ) Atm_AD%Cloud(n)%Water_Content(kc) = Atm_AD%Cloud(n)%Water_Content(kc) + & ( bs_AD * CSV%ke(kc,n) * CSV%w(kc,n) ) ! interpolation quality control IF( CSV%w(kc,n) >= ONE ) THEN w_AD = ZERO END IF IF( CSV%ke(kc,n) <= ZERO ) THEN ke_AD = ZERO w_AD = ZERO END IF IF( CSV%w(kc,n) <= ZERO ) THEN w_AD = ZERO pcoeff_AD = ZERO END IF ! Call sensor specific routines SELECT CASE (Sensor_Type) CASE (MICROWAVE_SENSOR) CALL Get_Cloud_Opt_MW_AD(CScat_AD , & ! Input Atm%Cloud(n)%Type , & ! Input ke_AD , & ! AD Input w_AD , & ! AD Input pcoeff_AD , & ! AD Input Atm_AD%Cloud(n)%Effective_Radius(kc), & ! AD Output Atm_AD%Temperature(kc) , & ! AD Output CSV%csi(kc,n) ) ! Interpolation CASE (INFRARED_SENSOR, VISIBLE_SENSOR) CALL Get_Cloud_Opt_IR_AD(CScat_AD , & ! Input Atm%Cloud(n)%Type , & ! Input ke_AD , & ! AD Input w_AD , & ! AD Input pcoeff_AD , & ! AD Input Atm_AD%Cloud(n)%Effective_Radius(kc), & ! AD Output CSV%csi(kc,n) ) ! Interpolation CASE DEFAULT ke_AD = ZERO w_AD = ZERO pcoeff_AD = ZERO END SELECT END DO Cloud_Layer_loop END DO Cloud_loop END FUNCTION CRTM_Compute_CloudScatter_AD !################################################################################ !################################################################################ !## ## !## ## PRIVATE MODULE ROUTINES ## ## !## ## !################################################################################ !################################################################################ ! ------------------------------------------ ! Subroutine to obtain the IR bulk ! optical properties of a cloud: ! extinction coefficient (ke), ! scattering coefficient (w) ! asymmetry factor (g), and ! spherical Legendre coefficients (pcoeff) ! ------------------------------------------ SUBROUTINE Get_Cloud_Opt_IR( CloudScatter, & ! Input CloudScatter structure Frequency , & ! Input Frequency (cm^-1) cloud_type , & ! Input see CRTM_Cloud_Define.f90 Reff , & ! Input effective radius (mm) ke , & ! Output optical depth for 1 mm water content w , & ! Output single scattering albedo pcoeff , & ! Output spherical Legendre coefficients csi ) ! Output interpolation data ! Arguments TYPE(CRTM_AtmScatter_type), INTENT(IN) :: CloudScatter REAL(fp) , INTENT(IN) :: Frequency INTEGER , INTENT(IN) :: Cloud_Type REAL(fp) , INTENT(IN) :: Reff REAL(fp) , INTENT(OUT) :: ke REAL(fp) , INTENT(OUT) :: w REAL(fp) , INTENT(IN OUT) :: pcoeff(0:,:) TYPE(CSinterp_type) , INTENT(IN OUT) :: csi ! Local variables INTEGER :: k, l ! Find the frequency and effective ! radius indices for interpolation ! -------------------------------- csi%f_int = MAX(MIN(CloudC%Frequency_IR(CloudC%n_IR_Frequencies),Frequency),CloudC%Frequency_IR(1)) CALL find_index(CloudC%Frequency_IR, csi%f_int, csi%i1,csi%i2, csi%f_outbound) csi%f = CloudC%Frequency_IR(csi%i1:csi%i2) csi%r_int = MAX(MIN(CloudC%Reff_IR(CloudC%n_IR_Radii),Reff),CloudC%Reff_IR(1)) CALL find_index(CloudC%Reff_IR, csi%r_int, csi%j1,csi%j2, csi%r_outbound) csi%r = CloudC%Reff_IR(csi%j1:csi%j2) ! Calculate the interpolating polynomials ! --------------------------------------- ! Frequency term CALL LPoly( csi%f, csi%f_int, & ! Input csi%wlp ) ! Output ! Effective radius term CALL LPoly( csi%r, csi%r_int, & ! Input csi%xlp ) ! Output ! Determine the density index, k, for the clouds ! based on CloudC LUT organisation ! ---------------------------------------------- SELECT CASE (Cloud_Type) CASE(WATER_CLOUD) ; k=0 ! Liquid CASE(ICE_CLOUD) ; k=3 ! Solid CASE(RAIN_CLOUD) ; k=0 ! Liquid CASE(SNOW_CLOUD) ; k=1 ! Solid CASE(GRAUPEL_CLOUD); k=2 ! Solid CASE(HAIL_CLOUD) ; k=3 ! Solid END SELECT ! Perform interpolation ! --------------------- CALL interp_2D( CloudC%ke_IR(csi%i1:csi%i2,csi%j1:csi%j2,k), csi%wlp, csi%xlp, ke ) CALL interp_2D( CloudC%w_IR(csi%i1:csi%i2,csi%j1:csi%j2,k) , csi%wlp, csi%xlp, w ) IF (CloudScatter%n_Phase_Elements > 0 ) THEN pcoeff(0,1) = POINT_5 DO l = 1, CloudScatter%n_Legendre_Terms CALL interp_2D( CloudC%pcoeff_IR(csi%i1:csi%i2,csi%j1:csi%j2,k,l+CloudScatter%lOffset), & csi%wlp, csi%xlp, pcoeff(l,1) ) END DO END IF END SUBROUTINE Get_Cloud_Opt_IR ! --------------------------------------------- ! Subroutine to obtain the tangent-linear ! IR bulk optical properties of a cloud: ! extinction coefficient (ke_TL), ! scattereing coefficient (w_TL) ! spherical Legendre coefficients (pcoeff_TL) ! --------------------------------------------- SUBROUTINE Get_Cloud_Opt_IR_TL( CloudScatter_TL, & ! Input CloudScatter TL structure cloud_type , & ! Input see CRTM_Cloud_Define.f90 Reff_TL , & ! TL Input effective radius (mm) ke_TL , & ! TL Output extinction coefficient (=~ optical depth for 1 mm water content) w_TL , & ! TL Output single scattering albedo pcoeff_TL , & ! TL Output spherical Legendre coefficients csi ) ! Input interpolation data ! Arguments TYPE(CRTM_AtmScatter_type), INTENT(IN) :: CloudScatter_TL INTEGER , INTENT(IN) :: Cloud_Type REAL(fp), INTENT(IN) :: Reff_TL REAL(fp), INTENT(OUT) :: ke_TL REAL(fp), INTENT(OUT) :: w_TL REAL(fp), INTENT(IN OUT) :: pcoeff_TL(0:,:) TYPE(CSinterp_type), INTENT(IN) :: csi ! Local variables INTEGER :: k, l REAL(fp) :: f_int_TL, r_int_TL REAL(fp) :: f_TL(NPTS), r_TL(NPTS) REAL(fp) :: z_TL(NPTS,NPTS) TYPE(LPoly_type) :: wlp_TL, xlp_TL REAL(fp), POINTER :: z(:,:) => NULL() ! Setup ! ----- ! No TL output when all dimensions ! are outside LUT bounds IF ( csi%f_outbound .AND. csi%r_outbound ) THEN ke_TL = ZERO w_TL = ZERO pcoeff_TL = ZERO RETURN END IF ! The TL inputs f_int_TL = ZERO f_TL = ZERO r_int_TL = Reff_TL r_TL = ZERO z_TL = ZERO ! Calculate the TL interpolating polynomials ! ------------------------------------------ ! Frequency term (always zero. This is a placeholder for testing) CALL LPoly_TL( csi%f, csi%f_int, & ! FWD Input csi%wlp, & ! FWD Input f_TL, f_int_TL, & ! TL Input wlp_TL ) ! TL Output ! Effective radius term CALL LPoly_TL( csi%r, csi%r_int, & ! FWD Input csi%xlp, & ! FWD Input r_TL, r_int_TL, & ! TL Input xlp_TL ) ! TL Output ! Determine the density index, k, for the clouds ! based on CloudC LUT organisation ! ---------------------------------------------- SELECT CASE (Cloud_Type) CASE(WATER_CLOUD) ; k=0 ! Liquid CASE(ICE_CLOUD) ; k=3 ! Solid CASE(RAIN_CLOUD) ; k=0 ! Liquid CASE(SNOW_CLOUD) ; k=1 ! Solid CASE(GRAUPEL_CLOUD); k=2 ! Solid CASE(HAIL_CLOUD) ; k=3 ! Solid END SELECT ! Perform interpolation based on cloud type ! ----------------------------------------- ! Extinction coefficient z => CloudC%ke_IR(csi%i1:csi%i2,csi%j1:csi%j2,k) CALL interp_2D_TL( z , csi%wlp, csi%xlp, & ! FWD Input z_TL, wlp_TL , xlp_TL , & ! TL Input ke_TL ) ! TL Output ! Single scatter albedo z => CloudC%w_IR(csi%i1:csi%i2,csi%j1:csi%j2,k) CALL interp_2D_TL( z , csi%wlp, csi%xlp, & ! FWD Input z_TL, wlp_TL , xlp_TL , & ! TL Input w_TL ) ! TL Output ! Phase matrix coefficients IF ( CloudScatter_TL%n_Phase_Elements > 0 ) THEN pcoeff_TL(0,1) = ZERO DO l = 1, CloudScatter_TL%n_Legendre_Terms z => CloudC%pcoeff_IR(csi%i1:csi%i2,csi%j1:csi%j2,k,l+CloudScatter_TL%lOffset) CALL interp_2D_TL( z , csi%wlp, csi%xlp, & ! FWD Input z_TL, wlp_TL , xlp_TL , & ! TL Input pcoeff_TL(l,1) ) ! TL Output END DO END IF NULLIFY(z) END SUBROUTINE Get_Cloud_Opt_IR_TL ! --------------------------------------------- ! Subroutine to obtain the adjoint of the ! IR bulk optical properties of a cloud: ! effective radius (Reff_AD) ! --------------------------------------------- SUBROUTINE Get_Cloud_Opt_IR_AD( CloudScatter_AD, & ! Input CloudScatter AD structure cloud_type , & ! Input see CRTM_Cloud_Define.f90 ke_AD , & ! AD Input extinction coefficient (=~ optical depth for 1 mm water content) w_AD , & ! AD Input single scattering albedo pcoeff_AD , & ! AD Input spherical Legendre coefficients Reff_AD , & ! AD Output effective radius (mm) csi ) ! Input interpolation data ! Arguments TYPE(CRTM_AtmScatter_type), INTENT(IN) :: CloudScatter_AD INTEGER , INTENT(IN) :: Cloud_Type REAL(fp), INTENT(IN OUT) :: ke_AD ! AD Input REAL(fp), INTENT(IN OUT) :: w_AD ! AD Input REAL(fp), INTENT(IN OUT) :: pcoeff_AD(0:,:) ! AD Input REAL(fp), INTENT(IN OUT) :: Reff_AD ! AD Output TYPE(CSinterp_type), INTENT(IN) :: csi ! Local variables INTEGER :: k, l REAL(fp) :: f_int_AD, r_int_AD REAL(fp) :: f_AD(NPTS), r_AD(NPTS) REAL(fp) :: z_AD(NPTS,NPTS) TYPE(LPoly_type) :: wlp_AD, xlp_AD REAL(fp), POINTER :: z(:,:) => NULL() ! Setup ! ----- ! No AD output when all dimensions ! are outside LUT bounds IF ( csi%f_outbound .AND. csi%r_outbound ) THEN Reff_AD = ZERO ke_AD = ZERO w_AD = ZERO pcoeff_AD = ZERO RETURN END IF ! Initialise local adjoint variables f_int_AD = ZERO r_int_AD = ZERO f_AD = ZERO r_AD = ZERO z_AD = ZERO CALL Clear_LPoly(wlp_AD) CALL Clear_LPoly(xlp_AD) ! Determine the density index, k, for the clouds ! based on CloudC LUT organisation ! ---------------------------------------------- SELECT CASE (Cloud_Type) CASE(WATER_CLOUD) ; k=0 ! Liquid CASE(ICE_CLOUD) ; k=3 ! Solid CASE(RAIN_CLOUD) ; k=0 ! Liquid CASE(SNOW_CLOUD) ; k=1 ! Solid CASE(GRAUPEL_CLOUD); k=2 ! Solid CASE(HAIL_CLOUD) ; k=3 ! Solid END SELECT ! Perform interpolation ! --------------------- ! Phase matrix coefficients IF (CloudScatter_AD%n_Phase_Elements > 0 ) THEN DO l = 1, CloudScatter_AD%n_Legendre_Terms z => CloudC%pcoeff_IR(csi%i1:csi%i2,csi%j1:csi%j2,k,l+CloudScatter_AD%lOffset) CALL interp_2D_AD( z , csi%wlp , csi%xlp , & ! FWD Input pcoeff_AD(l,1) , & ! AD Input z_AD, wlp_AD, xlp_AD ) ! AD Output END DO pcoeff_AD(0,1) = ZERO END IF ! Single scatter albedo z => CloudC%w_IR(csi%i1:csi%i2,csi%j1:csi%j2,k) CALL interp_2D_AD( z , csi%wlp, csi%xlp, & ! FWD Input w_AD , & ! AD Input z_AD, wlp_AD, xlp_AD ) ! AD Output ! Extinction coefficient z => CloudC%ke_IR(csi%i1:csi%i2,csi%j1:csi%j2,k) CALL interp_2D_AD( z , csi%wlp, csi%xlp, & ! FWD Input ke_AD , & ! AD Input z_AD, wlp_AD, xlp_AD ) ! AD Output NULLIFY(z) ! Compute the AD of the interpolating polynomials ! ----------------------------------------------- ! Efective radius term CALL LPoly_AD( csi%r, csi%r_int, & ! FWD Input csi%xlp, & ! FWD Input xlp_AD, & ! AD Input r_AD, r_int_AD ) ! AD Output ! Frequency term (always zero. This is a placeholder for testing) CALL LPoly_AD( csi%f, csi%f_int, & ! FWD Input csi%wlp, & ! FWD Input wlp_AD, & ! AD Input f_AD, f_int_AD ) ! AD Output ! The AD outputs ! -------------- Reff_AD = Reff_AD + r_int_AD END SUBROUTINE Get_Cloud_Opt_IR_AD ! ------------------------------------------ ! Subroutine to obtain the MW bulk ! optical properties of a cloud: ! extinction coefficient (ke), ! scattereing coefficient (w) ! asymmetry factor (g), and ! spherical Legendre coefficients (pcoeff) ! ------------------------------------------ SUBROUTINE Get_Cloud_Opt_MW( CloudScatter, & ! Input CloudScatter structure Frequency , & ! Input Frequency (GHz) cloud_type , & ! Input see CRTM_Cloud_Define.f90 Reff , & ! Input effective radius (mm) Temperature , & ! Input cloudy temperature ke , & ! Input optical depth for 1 mm water content w , & ! Input single scattering albedo pcoeff , & ! Output spherical Legendre coefficients csi ) ! Output interpolation data ! Arguments TYPE(CRTM_AtmScatter_type), INTENT(IN) :: CloudScatter REAL(fp) , INTENT(IN) :: Frequency INTEGER , INTENT(IN) :: Cloud_Type REAL(fp) , INTENT(IN) :: Reff REAL(fp) , INTENT(IN) :: Temperature REAL(fp) , INTENT(OUT) :: ke REAL(fp) , INTENT(OUT) :: w REAL(fp) , INTENT(IN OUT) :: pcoeff(0:,:) TYPE(CSinterp_type) , INTENT(IN OUT) :: csi ! Local variables INTEGER :: j, k, l, m ! Initialise results that may ! not be interpolated ! --------------------------- w = ZERO pcoeff = ZERO ! Find the frequency, effective radius ! and temperature indices for interpolation ! ----------------------------------------- csi%f_int = MAX(MIN(CloudC%Frequency_MW(CloudC%n_MW_Frequencies),Frequency),CloudC%Frequency_MW(1)) CALL find_index(CloudC%Frequency_MW, csi%f_int, csi%i1,csi%i2, csi%f_outbound) csi%f = CloudC%Frequency_MW(csi%i1:csi%i2) csi%r_int = MAX(MIN(CloudC%Reff_MW(CloudC%n_MW_Radii),Reff),CloudC%Reff_MW(1)) CALL find_index(CloudC%Reff_MW, csi%r_int, csi%j1,csi%j2, csi%r_outbound) csi%r = CloudC%Reff_MW(csi%j1:csi%j2) csi%t_int = MAX(MIN(CloudC%Temperature(CloudC%n_Temperatures),Temperature),CloudC%Temperature(1)) CALL find_index(CloudC%Temperature, csi%t_int, csi%k1,csi%k2, csi%t_outbound) csi%t = CloudC%Temperature(csi%k1:csi%k2) ! Calculate the interpolating polynomials ! --------------------------------------- ! Frequency term CALL LPoly( csi%f, csi%f_int, & ! Input csi%wlp ) ! Output ! Effective radius term CALL LPoly( csi%r, csi%r_int, & ! Input csi%xlp ) ! Output ! Temperature term CALL LPoly( csi%t, csi%t_int, & ! Input csi%ylp ) ! Output ! Perform interpolation based on cloud type ! ----------------------------------------- SELECT CASE (Cloud_Type) ! Only 2-D interpolation of extinction coefficient as a ! fn. of frequency and temperature for water cloud; i.e. ! we're only considering Rayleigh scattering here. CASE (WATER_CLOUD) j = 1 CALL interp_2D( CloudC%ke_L_MW(csi%i1:csi%i2,j,csi%k1:csi%k2), csi%wlp, csi%ylp, ke ) ! All 3-D interpolations for rain cloud! CASE (RAIN_CLOUD) CALL interp_3D( CloudC%ke_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2), csi%wlp, csi%xlp, csi%ylp, ke ) CALL interp_3D( CloudC%w_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2) , csi%wlp, csi%xlp, csi%ylp, w ) IF ( CloudScatter%n_Phase_Elements > 0 ) THEN pcoeff(0,1) = POINT_5 DO m = 1, CloudScatter%n_Phase_Elements DO l = 1, CloudScatter%n_Legendre_Terms CALL interp_3D( CloudC%pcoeff_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2,l+CloudScatter%lOffset,m), & csi%wlp, csi%xlp, csi%ylp, pcoeff(l,m) ) END DO END DO END IF ! Only 1-D interpolation of extinction coefficient as a ! fn. of frequency for ice cloud CASE (ICE_CLOUD) j = 1; k = 3 CALL interp_1D( CloudC%ke_S_MW(csi%i1:csi%i2,j,k), csi%wlp, ke ) ! The remaining cloud types have 2-D interpolation ! as a fn. of frequency and radius CASE DEFAULT ! Select the LUT density SELECT CASE (Cloud_Type) CASE (GRAUPEL_CLOUD); k = 2 CASE (HAIL_CLOUD) ; k = 3 CASE DEFAULT ; k = 1 END SELECT ! Perform interpolation CALL interp_2D( CloudC%ke_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k), csi%wlp, csi%xlp, ke ) CALL interp_2D( CloudC%w_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k) , csi%wlp, csi%xlp, w ) IF (CloudScatter%n_Phase_Elements > 0 ) THEN pcoeff(0,1) = POINT_5 DO m = 1, CloudScatter%n_Phase_Elements DO l = 1, CloudScatter%n_Legendre_Terms CALL interp_2D( CloudC%pcoeff_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k,l+CloudScatter%lOffset,m), & csi%wlp, csi%xlp, pcoeff(l,m) ) END DO END DO END IF END SELECT END SUBROUTINE Get_Cloud_Opt_MW ! --------------------------------------------- ! Subroutine to obtain the tangent-linear ! MW bulk optical properties of a cloud: ! extinction coefficient (ke_TL), ! scattereing coefficient (w_TL) ! spherical Legendre coefficients (pcoeff_TL) ! --------------------------------------------- SUBROUTINE Get_Cloud_Opt_MW_TL( CloudScatter_TL, & ! Input CloudScatter TL structure cloud_type , & ! Input see CRTM_Cloud_Define.f90 Reff_TL , & ! TL Input effective radius (mm) Temperature_TL , & ! TL Input cloudy temperature ke_TL , & ! TL Output extinction coefficient (=~ optical depth for 1 mm water content) w_TL , & ! TL Output single scattering albedo pcoeff_TL , & ! TL Output spherical Legendre coefficients csi ) ! Input interpolation data ! Arguments TYPE(CRTM_AtmScatter_type), INTENT(IN) :: CloudScatter_TL INTEGER , INTENT(IN) :: Cloud_Type REAL(fp), INTENT(IN) :: Reff_TL REAL(fp), INTENT(IN) :: Temperature_TL REAL(fp), INTENT(OUT) :: ke_TL REAL(fp), INTENT(OUT) :: w_TL REAL(fp), INTENT(IN OUT) :: pcoeff_TL(0:,:) TYPE(CSinterp_type), INTENT(IN) :: csi ! Local variables INTEGER :: j, k, l, m REAL(fp) :: f_int_TL, r_int_TL, t_int_TL REAL(fp) :: f_TL(NPTS), r_TL(NPTS), t_TL(NPTS) REAL(fp) :: z2_TL(NPTS,NPTS) REAL(fp) :: z3_TL(NPTS,NPTS,NPTS) TYPE(LPoly_type) :: wlp_TL, xlp_TL, ylp_TL REAL(fp), POINTER :: z2(:,:) => NULL() REAL(fp), POINTER :: z3(:,:,:) => NULL() ! Setup ! ----- ! Initialise results that may ! not be interpolated w_TL = ZERO pcoeff_TL = ZERO ! The local TL inputs f_int_TL = ZERO f_TL = ZERO r_int_TL = Reff_TL r_TL = ZERO t_int_TL = Temperature_TL t_TL = ZERO z2_TL = ZERO z3_TL = ZERO ! Calculate the TL interpolating polynomials ! ------------------------------------------ ! Frequency term (always zero. This is a placeholder for testing) CALL LPoly_TL( csi%f, csi%f_int, & ! FWD Input csi%wlp, & ! FWD Input f_TL, f_int_TL, & ! TL Input wlp_TL ) ! TL Output ! Effective radius term CALL LPoly_TL( csi%r, csi%r_int, & ! FWD Input csi%xlp, & ! FWD Input r_TL, r_int_TL, & ! TL Input xlp_TL ) ! TL Output ! Temperature term CALL LPoly_TL( csi%t, csi%t_int, & ! FWD Input csi%ylp, & ! FWD Input t_TL, t_int_TL, & ! TL Input ylp_TL ) ! TL Output ! Perform interpolation based on cloud type ! ----------------------------------------- SELECT CASE (Cloud_Type) ! Only 2-D interpolation of extinction coefficient as a ! fn. of frequency and temperature for water cloud; i.e. ! we're only considering Rayleigh scattering here. CASE (WATER_CLOUD) ! No TL output when all dimensions ! are outside LUT bounds IF ( csi%f_outbound .AND. csi%t_outbound ) THEN ke_TL = ZERO RETURN END IF j = 1 z2 => CloudC%ke_L_MW(csi%i1:csi%i2,j,csi%k1:csi%k2) CALL interp_2D_TL( z2 , csi%wlp, csi%ylp, & ! FWD Input z2_TL, wlp_TL , ylp_TL , & ! TL Input ke_TL ) ! TL Output ! All 3-D interpolations for rain cloud! CASE (RAIN_CLOUD) ! No TL output when all dimensions ! are outside LUT bounds IF ( csi%f_outbound .AND. csi%r_outbound .AND. csi%t_outbound ) THEN ke_TL = ZERO RETURN END IF ! Extinction coefficient z3 => CloudC%ke_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2) CALL interp_3D_TL( z3 , csi%wlp, csi%xlp, csi%ylp, & ! FWD Input z3_TL, wlp_TL , xlp_TL , ylp_TL , & ! TL Input ke_TL ) ! TL Output ! Single scatter albedo z3 => CloudC%w_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2) CALL interp_3D_TL( z3 , csi%wlp, csi%xlp, csi%ylp, & ! FWD Input z3_TL, wlp_TL , xlp_TL , ylp_TL , & ! TL Input w_TL ) ! TL Output ! Phase matrix coefficients IF ( CloudScatter_TL%n_Phase_Elements > 0 ) THEN pcoeff_TL(0,1) = ZERO DO m = 1, CloudScatter_TL%n_Phase_Elements DO l = 1, CloudScatter_TL%n_Legendre_Terms z3 => CloudC%pcoeff_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2,l+CloudScatter_TL%lOffset,m) CALL interp_3D_TL( z3 , csi%wlp, csi%xlp, csi%ylp, & ! FWD Input z3_TL, wlp_TL , xlp_TL , ylp_TL , & ! TL Input pcoeff_TL(l,m) ) ! TL Output END DO END DO END IF ! No TL interpolation of extinction coefficient as it ! is only a fn. of frequency for ice cloud CASE (ICE_CLOUD) ke_TL = ZERO ! The remaining cloud types have 2-D interpolation ! as a fn. of frequency and radius CASE DEFAULT ! No TL output when all dimensions ! are outside LUT bounds IF ( csi%f_outbound .AND. csi%r_outbound ) THEN ke_TL = ZERO RETURN END IF ! Select the LUT temperature SELECT CASE (Cloud_Type) CASE (GRAUPEL_CLOUD); k = 2 CASE (HAIL_CLOUD) ; k = 3 CASE DEFAULT ; k = 1 END SELECT ! Extinction coefficient z2 => CloudC%ke_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k) CALL interp_2D_TL( z2 , csi%wlp, csi%xlp, & ! FWD Input z2_TL, wlp_TL , xlp_TL , & ! TL Input ke_TL ) ! TL Output ! Single scatter albedo z2 => CloudC%w_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k) CALL interp_2D_TL( z2 , csi%wlp, csi%xlp, & ! FWD Input z2_TL, wlp_TL , xlp_TL , & ! TL Input w_TL ) ! TL Output ! Phase matrix coefficients IF ( CloudScatter_TL%n_Phase_Elements > 0 ) THEN pcoeff_TL(0,1) = ZERO DO m = 1, CloudScatter_TL%n_Phase_Elements DO l = 1, CloudScatter_TL%n_Legendre_Terms z2 => CloudC%pcoeff_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k,l+CloudScatter_TL%lOffset,m) CALL interp_2D_TL( z2 , csi%wlp, csi%xlp, & ! FWD Input z2_TL, wlp_TL , xlp_TL , & ! TL Input pcoeff_TL(l,m) ) ! TL Output END DO END DO END IF END SELECT NULLIFY(z2, z3) END SUBROUTINE Get_Cloud_Opt_MW_TL ! --------------------------------------------- ! Subroutine to obtain the adjoint of the ! MW bulk optical properties of a cloud: ! effective radius (Reff_AD), ! temperature (temperature_AD) ! --------------------------------------------- SUBROUTINE Get_Cloud_Opt_MW_AD(CloudScatter_AD, & ! Input CloudScatter AD structure cloud_type , & ! Input see CRTM_Cloud_Define.f90 ke_AD , & ! AD Input extinction coefficient (=~ optical depth for 1 mm water content) w_AD , & ! AD Input single scattering albedo pcoeff_AD , & ! AD Input spherical Legendre coefficients Reff_AD , & ! AD Output effective radius (mm) Temperature_AD , & ! AD Output temperature csi ) ! Input interpolation data ! Arguments TYPE(CRTM_AtmScatter_type), INTENT(IN) :: CloudScatter_AD INTEGER , INTENT(IN) :: Cloud_Type REAL(fp), INTENT(IN OUT) :: ke_AD ! AD Input REAL(fp), INTENT(IN OUT) :: w_AD ! AD Input REAL(fp), INTENT(IN OUT) :: pcoeff_AD(0:,:) ! AD Input REAL(fp), INTENT(IN OUT) :: Reff_AD ! AD Output REAL(fp), INTENT(IN OUT) :: Temperature_AD ! AD Output TYPE(CSinterp_type), INTENT(IN) :: csi ! Local variables INTEGER :: j, k, l, m REAL(fp) :: f_int_AD, r_int_AD, t_int_AD REAL(fp) :: f_AD(NPTS), r_AD(NPTS), t_AD(NPTS) REAL(fp) :: z2_AD(NPTS,NPTS) REAL(fp) :: z3_AD(NPTS,NPTS,NPTS) TYPE(LPoly_type) :: wlp_AD, xlp_AD, ylp_AD REAL(fp), POINTER :: z2(:,:) => NULL() REAL(fp), POINTER :: z3(:,:,:) => NULL() ! Setup ! ----- ! Initialise local adjoint variables f_int_AD = ZERO f_AD = ZERO r_int_AD = ZERO r_AD = ZERO t_int_AD = ZERO t_AD = ZERO z2_AD = ZERO z3_AD = ZERO CALL Clear_LPoly(wlp_AD) CALL Clear_LPoly(xlp_AD) CALL Clear_LPoly(ylp_AD) ! Perform interpolation based on cloud type ! ----------------------------------------- SELECT CASE (Cloud_Type) ! Only 2-D interpolation of extinction coefficient as a ! fn. of frequency and temperature for water cloud; i.e. ! we're only considering Rayleigh scattering here. ! ------------------------------------------------------ CASE (WATER_CLOUD) ! No AD output when all dimensions ! are outside LUT bounds IF ( csi%f_outbound .AND. csi%t_outbound ) THEN ke_AD = ZERO w_AD = ZERO pcoeff_AD = ZERO RETURN END IF ! Perform the AD interpolations j = 1 z2 => CloudC%ke_L_MW(csi%i1:csi%i2,j,csi%k1:csi%k2) CALL interp_2D_AD( z2 , csi%wlp, csi%ylp, & ! FWD Input ke_AD , & ! AD Input z2_AD, wlp_AD, ylp_AD ) ! AD Output ! Compute the AD of the interpolating polynomials ! Temperature term CALL LPoly_AD( csi%t, csi%t_int, & ! FWD Input csi%ylp, & ! FWD Input ylp_AD, & ! AD Input t_AD, t_int_AD ) ! AD Output ! Frequency term (always zero. This is a placeholder for testing) CALL LPoly_AD( csi%f, csi%f_int, & ! FWD Input csi%wlp, & ! FWD Input wlp_AD, & ! AD Input f_AD, f_int_AD ) ! AD Output ! The AD outputs Temperature_AD = Temperature_AD + t_int_AD ! All 3-D interpolations for rain cloud! ! -------------------------------------- CASE (RAIN_CLOUD) ! No AD output when all dimensions ! are outside LUT bounds IF ( csi%f_outbound .AND. csi%r_outbound .AND. csi%t_outbound ) THEN ke_AD = ZERO w_AD = ZERO pcoeff_AD = ZERO RETURN END IF ! Perform the AD interpolations ! Phase matrix coefficients IF (CloudScatter_AD%n_Phase_Elements > 0 ) THEN DO m = 1, CloudScatter_AD%n_Phase_Elements DO l = 1, CloudScatter_AD%n_Legendre_Terms z3 => CloudC%pcoeff_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2,l+CloudScatter_AD%lOffset,m) CALL interp_3D_AD( z3 , csi%wlp, csi%xlp, csi%ylp, & ! FWD Input pcoeff_AD(l,m) , & ! AD Input z3_AD, wlp_AD , xlp_AD , ylp_AD ) ! AD Output END DO END DO pcoeff_AD(0,1) = ZERO END IF ! Single scatter albedo z3 => CloudC%w_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2) CALL interp_3D_AD( z3 , csi%wlp, csi%xlp, csi%ylp, & ! FWD Input w_AD , & ! AD Input z3_AD, wlp_AD , xlp_AD , ylp_AD ) ! AD Output ! Extinction coefficient z3 => CloudC%ke_L_MW(csi%i1:csi%i2,csi%j1:csi%j2,csi%k1:csi%k2) CALL interp_3D_AD( z3 , csi%wlp, csi%xlp, csi%ylp, & ! FWD Input ke_AD , & ! AD Input z3_AD, wlp_AD , xlp_AD , ylp_AD ) ! AD Output ! Compute the AD of the interpolating polynomials ! Temperature term CALL LPoly_AD( csi%t, csi%t_int, & ! FWD Input csi%ylp, & ! FWD Input ylp_AD, & ! AD Input t_AD, t_int_AD ) ! AD Output ! Effective radius term CALL LPoly_AD( csi%r, csi%r_int, & ! FWD Input csi%xlp, & ! FWD Input xlp_AD, & ! AD Input r_AD, r_int_AD ) ! AD Output ! Frequency term (always zero. This is a placeholder for testing) CALL LPoly_AD( csi%f, csi%f_int, & ! FWD Input csi%wlp, & ! FWD Input wlp_AD, & ! AD Input f_AD, f_int_AD ) ! AD Output ! The AD outputs Temperature_AD = Temperature_AD + t_int_AD Reff_AD = Reff_AD + r_int_AD ! No AD interpolation as it is only a fn. ! of frequency for ice cloud ! --------------------------------------- CASE (ICE_CLOUD) ke_AD = ZERO w_AD = ZERO pcoeff_AD = ZERO ! The remaining cloud types have 2-D interpolation ! as a fn. of frequency and radius ! ------------------------------------------------ CASE DEFAULT ! No TL output when all dimensions ! are outside LUT bounds IF ( csi%f_outbound .AND. csi%r_outbound ) THEN ke_AD = ZERO w_AD = ZERO pcoeff_AD = ZERO RETURN END IF ! Select the LUT temperature SELECT CASE (Cloud_Type) CASE (GRAUPEL_CLOUD); k = 2 CASE (HAIL_CLOUD) ; k = 3 CASE DEFAULT ; k = 1 END SELECT ! Perform the AD interpolations ! Phase matrix coefficients IF (CloudScatter_AD%n_Phase_Elements > 0 ) THEN DO m = 1, CloudScatter_AD%n_Phase_Elements DO l = 1, CloudScatter_AD%n_Legendre_Terms z2 => CloudC%pcoeff_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k,l+CloudScatter_AD%lOffset,m) CALL interp_2D_AD( z2 , csi%wlp, csi%xlp, & ! FWD Input pcoeff_AD(l,m) , & ! AD Input z2_AD, wlp_AD , xlp_AD ) ! AD Output END DO END DO pcoeff_AD(0,1) = ZERO END IF ! Single scatter albedo z2 => CloudC%w_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k) CALL interp_2D_AD( z2 , csi%wlp, csi%xlp, & ! FWD Input w_AD , & ! AD Input z2_AD, wlp_AD , xlp_AD ) ! AD Output ! Extinction coefficient z2 => CloudC%ke_S_MW(csi%i1:csi%i2,csi%j1:csi%j2,k) CALL interp_2D_AD( z2 , csi%wlp, csi%xlp, & ! FWD Input ke_AD , & ! AD Input z2_AD, wlp_AD , xlp_AD ) ! AD Output ! Compute the AD of the interpolating polynomials ! Effective radius term CALL LPoly_AD( csi%r, csi%r_int, & ! FWD Input csi%xlp, & ! FWD Input xlp_AD, & ! AD Input r_AD, r_int_AD ) ! AD Output ! Frequency term (always zero. This is a placeholder for testing) CALL LPoly_AD( csi%f, csi%f_int, & ! FWD Input csi%wlp, & ! FWD Input wlp_AD, & ! AD Input f_AD, f_int_AD ) ! AD Output ! The AD outputs Reff_AD = Reff_AD + r_int_AD END SELECT NULLIFY(z2, z3) END SUBROUTINE Get_Cloud_Opt_MW_AD END MODULE CRTM_CloudScatter