! ! CRTM_Atmosphere ! ! Module for adding layers to the CRTM atmosphere structure as required. ! ! ! CREATION HISTORY: ! Written by: Paul van Delst, 29-Oct-2007 ! paul.vandelst@noaa.gov ! MODULE CRTM_Atmosphere ! ----------------- ! Environment setup ! ----------------- ! Module use USE Type_Kinds , ONLY: fp USE Message_Handler , ONLY: SUCCESS, FAILURE, WARNING, Display_Message USE CRTM_Parameters , ONLY: ZERO, ONE, POINT_5, SET, & TOA_PRESSURE , & MINIMUM_ABSORBER_AMOUNT USE CRTM_Atmosphere_Define, ONLY: CRTM_Atmosphere_type , & OPERATOR(==), & OPERATOR(+), & CRTM_Atmosphere_Associated, & CRTM_Atmosphere_Create, & CRTM_Atmosphere_AddLayerCopy, & CRTM_Atmosphere_Zero USE CRTM_Model_Profiles , ONLY: MODEL_LEVEL_PRESSURE, & CRTM_Get_Model_Profile ! ...Internal variable definition module USE iAtm_Define, ONLY: iAtm_type , & iAtm_Associated, & iAtm_Create , & iAtm_Destroy ! Disable implicit typing IMPLICIT NONE ! ------------ ! Visibilities ! ------------ ! Everything private by default PRIVATE ! Module procedures PUBLIC :: CRTM_Atmosphere_AddLayers PUBLIC :: CRTM_Atmosphere_AddLayers_TL PUBLIC :: CRTM_Atmosphere_AddLayers_AD ! iAtm entities ! ...Structure PUBLIC :: iAtm_type ! ...Procedures PUBLIC :: iAtm_Associated PUBLIC :: iAtm_Create PUBLIC :: iAtm_Destroy ! ----------------- ! Module parameters ! ----------------- ! RCS Id for the module CHARACTER(*), PARAMETER :: MODULE_RCS_ID = & '$Id: CRTM_Atmosphere.f90 22707 2012-11-21 21:09:10Z paul.vandelst@noaa.gov $' ! Message string length INTEGER, PARAMETER :: ML = 256 CONTAINS !################################################################################ !################################################################################ !## ## !## ## PUBLIC MODULE ROUTINES ## ## !## ## !################################################################################ !################################################################################ !-------------------------------------------------------------------------------- !:sdoc+: ! ! NAME: ! CRTM_Atmosphere_AddLayers ! ! PURPOSE: ! Function to copy an atmosphere structure and adding extra layers from ! climatology as required to supplement the upper atmosphere profile data. ! ! CALLING SEQUENCE: ! Error_Status = CRTM_Atmosphere_AddLayers( Atm_In , & ! Input ! Atm_Out ) ! Output ! ! INPUTS: ! Atm_In: Atmosphere structure that is to be supplemented ! or copied. ! UNITS: N/A ! TYPE: CRTM_Atmosphere_type ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! OUTPUTS: ! Atm_Out: Copy of the input atmosphere structure with extra upper ! atmosphere layers added as required. ! UNITS: N/A ! 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 operation was successful ! == FAILURE an error occurred ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! !:sdoc-: !-------------------------------------------------------------------------------- FUNCTION CRTM_Atmosphere_AddLayers( & Atm_In , & ! Input Atm_Out) & ! Output RESULT( err_stat ) ! Arguments TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm_In TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_Out ! Function result INTEGER :: err_stat ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_AddLayers' ! Local variables CHARACTER(ML) :: msg INTEGER :: i, j, k, n TYPE(iAtm_type) :: iAtm ! Set up err_stat = SUCCESS ! If extra layers are NOT needed, ! then simply copy the structure IF ( Atm_In%Level_Pressure(0) <= TOA_PRESSURE) THEN Atm_Out = Atm_In IF ( .NOT. CRTM_Atmosphere_Associated( Atm_Out ) ) THEN err_stat = FAILURE msg = 'Error assigning Atmosphere structure with NO extra layers' CALL Display_Message( ROUTINE_NAME, msg, err_stat ) END IF RETURN END IF ! Determine the number of extra layers required n = Extra_Layers( Atm_In ) IF ( n < 1 ) THEN err_stat = FAILURE msg = 'Error determining extra layer count' CALL Display_Message( ROUTINE_NAME, msg, err_stat ) RETURN END IF ! Allocate the internal variable structure CALL iAtm_Create( iAtm, n, Atm_In%n_Absorbers ) IF ( .NOT. iAtm_Associated( iAtm ) ) THEN err_stat = FAILURE msg = 'Error allocating iAtm internal structure' CALL Display_Message( ROUTINE_NAME, msg, err_stat ) RETURN END IF ! Get the extra layer profiles CALL CRTM_Get_Model_Profile( iAtm%pl, iAtm%tl, iAtm%al, Model=Atm_In%Climatology ) ! First interpolate the extra levels to the user top pressure ! replacing the model data at that array index CALL Interp_LPoly( Atm_In%Level_Pressure(0), iAtm%pl(n-1:n), iAtm%ilpoly ) iAtm%plint_save = Atm_In%Level_Pressure(0) iAtm%pln_save = iAtm%pl(n) iAtm%pl(n) = iAtm%plint_save CALL Interp_Linear( iAtm%ilpoly, iAtm%tl(n-1:n), iAtm%tlint_save ) iAtm%tln_save = iAtm%tl(n) iAtm%tl(n) = iAtm%tlint_save DO j = 1, Atm_In%n_Absorbers CALL Interp_Linear( iAtm%ilpoly, iAtm%al(n-1:n,j), iAtm%alint_save(j) ) iAtm%aln_save(j) = iAtm%al(n,j) iAtm%al(n,j) = iAtm%alint_save(j) END DO ! Now compute the model profile layer averages DO k = 1, n CALL Layer_P(iAtm%pl(k-1:k), iAtm%p(k)) CALL Layer_X(iAtm%tl(k-1:k), iAtm%t(k)) END DO DO j = 1, Atm_In%n_Absorbers DO k = 1, n CALL Layer_X(iAtm%al(k-1:k,j), iAtm%a(k,j)) END DO END DO ! Now, extrapolate user layer profile to get the "layer 0" value and ! use it to shift the model profile to avoid a discontinuity at p(n) CALL Interp_LPoly( iAtm%p(n), Atm_In%Pressure(1:2), iAtm%elpoly ) CALL Shift_Profile( iAtm%elpoly, Atm_In%Temperature(1:2), iAtm%t ) DO j = 1, Atm_In%n_Absorbers CALL Shift_Profile( iAtm%elpoly, Atm_In%Absorber(1:2,j), iAtm%a(:,j) ) END DO ! Make sure the absorber amounts are not negative. ! (Is a further, more physical, check needed here?) iAtm%a_save = iAtm%a WHERE (iAtm%a_save < ZERO) iAtm%a = MINIMUM_ABSORBER_AMOUNT ! Copy over the atmosphere structure with extra layers atm_out = CRTM_Atmosphere_AddLayerCopy( Atm_In, n ) IF ( .NOT. CRTM_Atmosphere_Associated( atm_out ) ) THEN err_stat = FAILURE msg = 'Error copying Atmosphere structure with extra layers' CALL Display_Message( ROUTINE_NAME, msg, err_stat ) RETURN END IF ! Slot the added layers into the output atmosphere structure ! Note: Cloud and Aerosol assignments not really needed (the ! zeroing is handled by the structure allocation) since ! at TOA, typically, there are not any clouds and/or ! aerosols. ! ...Profile Atm_Out%Level_Pressure(0:n) = iAtm%pl Atm_Out%Pressure(1:n) = iAtm%p Atm_Out%Temperature(1:n) = iAtm%t DO j = 1, Atm_Out%n_Absorbers Atm_Out%Absorber(1:n,j) = iAtm%a(:,j) END DO ! ...Clouds IF ( Atm_In%n_Clouds > 0 ) THEN DO i = 1, Atm_In%n_Clouds Atm_Out%Cloud(i)%Effective_Radius(1:n) = ZERO Atm_Out%Cloud(i)%Effective_Variance(1:n) = ZERO Atm_Out%Cloud(i)%Water_Content(1:n) = ZERO END DO END IF ! ...Aerosols IF ( Atm_In%n_Aerosols > 0 ) THEN DO i = 1, Atm_In%n_Aerosols Atm_Out%Aerosol(i)%Effective_Radius(1:n) = ZERO Atm_Out%Aerosol(i)%Concentration(1:n) = ZERO END DO END IF ! Clean up CALL iAtm_Destroy( iAtm ) END FUNCTION CRTM_Atmosphere_AddLayers !-------------------------------------------------------------------------------- !:sdoc+: ! ! NAME: ! CRTM_Atmosphere_AddLayers_TL ! ! PURPOSE: ! Function to copy a tangent-linear atmosphere structure and add extra ! layers as required to supplement the upper atmosphere profile data. ! ! CALLING SEQUENCE: ! Error_Status = CRTM_Atmosphere_AddLayers_TL( Atm_In , & ! FWD Input ! Atm_In_TL , & ! TL Input ! Atm_Out_TL ) ! TL Output ! ! INPUTS: ! Atm_In: Forward model atmosphere structure. ! UNITS: N/A ! TYPE: CRTM_Atmosphere_type ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! Atm_In_TL: Tangent-linear model atmosphere structure that is ! to be copied. ! UNITS: N/A ! TYPE: CRTM_Atmosphere_type ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! OUTPUTS: ! Atm_Out_TL: Copy of the input tangent-linear atmosphere structure ! with extra upper atmosphere layers added as required. ! Note that the tangent-linear values of the added layers ! is *always* zero. ! UNITS: N/A ! 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 operation was successful ! == FAILURE an error occurred ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! !:sdoc-: !-------------------------------------------------------------------------------- FUNCTION CRTM_Atmosphere_AddLayers_TL( & Atm_In , & ! FWD Input Atm_In_TL , & ! TL Input Atm_Out_TL) & ! TL Output RESULT( err_stat ) ! Arguments TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm_In TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm_In_TL TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_Out_TL ! Function result INTEGER :: err_stat ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_AddLayers_TL' ! Local variables CHARACTER(ML) :: msg INTEGER :: n ! Set up err_stat = SUCCESS ! If extra layers are NOT needed, ! then simply copy the structure IF ( Atm_In%Level_Pressure(0) <= TOA_PRESSURE) THEN Atm_Out_TL = Atm_In_TL IF ( .NOT. CRTM_Atmosphere_Associated( Atm_Out_TL ) ) THEN err_stat = FAILURE msg = 'Error assigning Atmosphere structure with NO extra layers' CALL Display_Message( ROUTINE_NAME, msg, err_stat ) END IF RETURN END IF ! Determine how many extra layers are needed n = Extra_Layers( Atm_In ) IF ( n < 1 ) THEN err_stat = FAILURE msg = 'Error determining extra layer count' CALL Display_Message( ROUTINE_NAME, msg, err_stat ) RETURN END IF ! Copy over the atmosphere structure with extra layers ! (which will be zero by definition) atm_out_TL = CRTM_Atmosphere_AddLayerCopy( Atm_In_TL, n ) IF ( .NOT. CRTM_Atmosphere_Associated( atm_out_TL ) ) THEN err_stat = FAILURE msg = 'Error copying Atmosphere structure with extra layers' CALL Display_Message( ROUTINE_NAME, msg, err_stat ) RETURN END IF END FUNCTION CRTM_Atmosphere_AddLayers_TL !-------------------------------------------------------------------------------- !:sdoc+: ! ! NAME: ! CRTM_Atmosphere_AddLayers_AD ! ! PURPOSE: ! Function to copy back an adjoint atmosphere structure removing added ! extra layers as were required to supplement the upper atmosphere ! profile data. ! ! CALLING SEQUENCE: ! Error_Status = CRTM_Atmosphere_AddLayers_AD( Atm_In , & ! FWD Input ! Atm_Out_AD, & ! AD Input ! Atm_In_AD ) ! AD Output ! ! INPUTS: ! Atm_In: Forward model atmosphere structure. ! UNITS: N/A ! TYPE: CRTM_Atmosphere_type ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN) ! ! Atm_Out_AD: Adjoint atmosphere structure that contains the added ! extra layers. ! ** SET TO ZERO ON EXIT ** ! UNITS: N/A ! TYPE: CRTM_Atmosphere_type ! DIMENSION: Scalar ! ATTRIBUTES: INTENT(IN OUT) ! ! OUTPUTS: ! Atm_In_AD: Adjoint atmosphere structure at the original, user ! specified layering. ! ** MUST HAVE VALUE ON ENTRY ** ! UNITS: N/A ! 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 operation was successful ! == FAILURE an error occurred ! UNITS: N/A ! TYPE: INTEGER ! DIMENSION: Scalar ! !:sdoc-: !-------------------------------------------------------------------------------- FUNCTION CRTM_Atmosphere_AddLayers_AD( & Atm_In , & ! FWD Input Atm_Out_AD, & ! AD Input Atm_In_AD ) & ! AD Output RESULT( err_stat ) ! Arguments TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm_In TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_Out_AD TYPE(CRTM_Atmosphere_type), INTENT(IN OUT) :: Atm_In_AD ! Function result INTEGER :: err_stat ! Local parameters CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'CRTM_Atmosphere_AddLayers_AD' ! Local variables CHARACTER(ML) :: msg INTEGER :: i, j, n, no, nt ! Set up err_stat = SUCCESS ! If extra layers are NOT needed, then simply perform ! the adjoint sum. Remember the TL form is ! Atm_Out_TL = Atm_In_TL ! so the adjoint form is ! Atm_In_AD = Atm_In_AD + Atm_Out_AD ! Atm_Out_AD = ZERO IF ( Atm_In%Level_Pressure(0) <= TOA_PRESSURE) THEN Atm_In_AD = Atm_In_AD + Atm_Out_AD CALL CRTM_Atmosphere_Zero( Atm_Out_AD ) RETURN END IF ! Determine how many extra layers have been used n = Extra_Layers( Atm_In ) IF ( n < 1 ) THEN err_stat = FAILURE msg = 'Error determining extra layer count' CALL Display_Message( ROUTINE_NAME, msg, err_stat ) RETURN END IF ! Perform the adjoint summations ! This is the adjoint equivalent of the TL Assign_Atmosphere no = Atm_In_AD%n_Layers nt = n + no ! ...Aerosols IF ( Atm_In_AD%n_Aerosols > 0 ) THEN DO i = 1, Atm_In_AD%n_Aerosols Atm_In_AD%Aerosol(i)%Concentration(1:no) = Atm_In_AD%Aerosol(i)%Concentration(1:no) + & Atm_Out_AD%Aerosol(i)%Concentration(n+1:nt) Atm_In_AD%Aerosol(i)%Effective_Radius(1:no) = Atm_In_AD%Aerosol(i)%Effective_Radius(1:no) + & Atm_Out_AD%Aerosol(i)%Effective_Radius(n+1:nt) Atm_In_AD%Aerosol(i)%Type = Atm_Out_AD%Aerosol(i)%Type END DO END IF ! ...Clouds IF ( Atm_In_AD%n_Clouds > 0 ) THEN DO i = 1, Atm_In_AD%n_Clouds Atm_In_AD%Cloud(i)%Water_Content(1:no) = Atm_In_AD%Cloud(i)%Water_Content(1:no) + & Atm_Out_AD%Cloud(i)%Water_Content(n+1:nt) Atm_In_AD%Cloud(i)%Effective_Variance(1:no) = Atm_In_AD%Cloud(i)%Effective_Variance(1:no) + & Atm_Out_AD%Cloud(i)%Effective_Variance(n+1:nt) Atm_In_AD%Cloud(i)%Effective_Radius(1:no) = Atm_In_AD%Cloud(i)%Effective_Radius(1:no) + & Atm_Out_AD%Cloud(i)%Effective_Radius(n+1:nt) Atm_In_AD%Cloud(i)%Type = Atm_Out_AD%Cloud(i)%Type END DO END IF ! ...Absorber data DO j = 1, Atm_In_AD%n_Absorbers Atm_In_AD%Absorber(1:no,j) = Atm_In_AD%Absorber(1:no,j) + Atm_Out_AD%Absorber(n+1:nt,j) END DO ! ...Temperature data Atm_In_AD%Temperature(1:no) = Atm_In_AD%Temperature(1:no) + Atm_Out_AD%Temperature(n+1:nt) ! ...Pressure data Atm_In_AD%Pressure(1:no) = Atm_In_AD%Pressure(1:no) + Atm_Out_AD%Pressure(n+1:nt) Atm_In_AD%Level_Pressure(0:no) = Atm_In_AD%Level_Pressure(0:no) + Atm_Out_AD%Level_Pressure(n:nt) ! Zero the output atmosphere structure CALL CRTM_Atmosphere_Zero( Atm_Out_AD ) END FUNCTION CRTM_Atmosphere_AddLayers_AD !################################################################################## !################################################################################## !## ## !## ## PRIVATE MODULE ROUTINES ## ## !## ## !################################################################################## !################################################################################## ! Subprogram to determine the number of extra layers required FUNCTION Extra_Layers( Atm ) RESULT( n ) TYPE(CRTM_Atmosphere_type), INTENT(IN) :: Atm INTEGER :: n DO n = 1, SIZE(MODEL_LEVEL_PRESSURE) IF ( MODEL_LEVEL_PRESSURE(n) >= Atm%Level_Pressure(0) ) RETURN END DO n = 0 END FUNCTION Extra_Layers ! Subprogram to compute the average layer pressure SUBROUTINE Layer_P( p, p_layer ) REAL(fp), INTENT(IN) :: p(2) ! Input REAL(fp), INTENT(OUT) :: p_layer ! Output p_layer = (p(2)-p(1))/LOG(p(2)/p(1)) END SUBROUTINE Layer_P ! Subprogram to compute the average layer amount of X SUBROUTINE Layer_X( x, x_layer ) REAL(fp), INTENT(IN) :: x(2) REAL(fp), INTENT(OUT) :: x_layer x_layer = POINT_5*(x(1)+x(2)) END SUBROUTINE Layer_X ! Subprogram to compute the interpolating polynomial linear in log(p) SUBROUTINE Interp_LPoly( p_int, p, lpoly ) REAL(fp), INTENT(IN) :: p_int REAL(fp), INTENT(IN) :: p(2) REAL(fp), INTENT(OUT) :: lpoly lpoly = (LOG(p_int)-LOG(p(1))) / (LOG(p(2))-LOG(p(1))) END SUBROUTINE Interp_LPoly ! Subprogram to perform linear interpolation SUBROUTINE Interp_Linear( lpoly, x, x_int ) REAL(fp), INTENT(IN) :: lpoly REAL(fp), INTENT(IN) :: x(2) REAL(fp), INTENT(OUT) :: x_int x_int = (x(2)-x(1))*lpoly + x(1) END SUBROUTINE Interp_Linear ! Subprogram to shifted the added profile layers SUBROUTINE Shift_Profile( lpoly, x_toa, x_shifted ) REAL(fp), INTENT(IN) :: lpoly REAL(fp), INTENT(IN) :: x_toa(2) REAL(fp), INTENT(IN OUT) :: x_shifted(:) INTEGER :: n REAL(fp) :: x_int, dx n = SIZE(x_shifted) CALL Interp_Linear( lpoly, x_toa, x_int ) dx = x_int - x_shifted(n) x_shifted = x_shifted + dx END SUBROUTINE Shift_Profile END MODULE CRTM_Atmosphere