! $Id: rt_utilities_mod.f90 3599 2019-11-20 18:49:25Z yli $ !-------------------------------------------------------------------------------------- ! Clouds from AVHRR Extended (CLAVR-x) 1b PROCESSING SOFTWARE Version 6.0 ! ! NAME: rt_utilities_mod.f90 (src) ! RT_UTILITIES_MOD (program) ! ! PURPOSE: Perform needed Radiative Transfer Functions ! ! DESCRIPTION: CLAVR-x uses much radiative transfer. This module stores and serves ! the RT variables to other modules. It serves as the interface ! to external RT models (PFAAST) and holds convenient RT-specific ! routines ! ! AUTHORS: ! Andrew Heidinger, Andrew.Heidinger@noaa.gov ! Andi Walther, CIMSS, andi.walther@ssec.wisc.edu ! Denis Botambekov, CIMSS, denis.botambekov@ssec.wisc.edu ! William Straka, CIMSS, wstraka@ssec.wisc.edu ! ! COPYRIGHT ! THIS SOFTWARE AND ITS DOCUMENTATION ARE CONSIDERED TO BE IN THE PUBLIC ! DOMAIN AND THUS ARE AVAILABLE FOR UNRESTRICTED PUBLIC USE. THEY ARE ! FURNISHED "AS IS." THE AUTHORS, THE UNITED STATES GOVERNMENT, ITS ! INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND AGENTS MAKE NO WARRANTY, ! EXPRESS OR IMPLIED, AS TO THE USEFULNESS OF THE SOFTWARE AND ! DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME NO RESPONSIBILITY (1) FOR ! THE USE OF THE SOFTWARE AND DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL ! SUPPORT TO USERS. ! ! Description of RTM Structure Members given in rtm_common.f90 ! ! CLAVR-x has 44 channels. ! Channels 1-36 are MODIS ! Channels 37-38 are ABI channels not on MODIS ! Channels 39-43 are the VIIRS I-bands ! Channel 44 is the VIIRS DNB ! ! Not all members of the RTM structure are populated for all channels. ! However, all members are allocated for any active cell ! ! Here is the current implementation ! ! There are 6 types of channels (FIX THIS) ! 1. MODIS IR-only channels = 21-36 excluding 26. ! 2. Non-MODIS IR-only channels 37-38 ! 3. MODIS (Solar + IR) channels = Channel 20 ! 4. Supported Solar Channels = Channels 1,2,5,6,7 and 44 ! 5. Unsupported Solar Channels = Channels 3,4,8-19,26,39-41 ! 6. Unsupported IR Channels = Channels 42 and 43 ! ! For each type described above, the following profiles are made: ! 1. Rad_Atm_Prof, Trans_Atm_Profile, Rad_BB_Cloud_Profile ! 2. All Profiles ! 3. Trans_Atm_Profile, Trans_Atm_Solar_Profile, Trans_Atm_Total_Profile ! 4. No Profiles ! 5. No Profiles ! ! Viirs I-band support is limited (channels 39-43) and is being developed. No ! pixel-level toa clear-sky fields are generated for the I-band variables yet. ! !-------------------------------------------------------------------------------------- module RT_UTILITIES_MOD use CONSTANTS_MOD, only: & Real4 & , Int4 & , Int1 & , Sym & , Dtor & , Missing_Value_Real4 & , Exe_Prompt & , G & , PI & , SOLAR_OBS_TYPE & , LUNAR_OBS_TYPE & , MIXED_OBS_TYPE & , THERMAL_OBS_TYPE & , Nchan_Clavrx use NWP_COMMON_MOD use SURFACE_PROPERTIES_MOD use PIXEL_COMMON_MOD, only: & NWP_PIX & , Sensor & , Ch & , Geo & , Sfc & , Image & , Bad_Pixel_Mask & , Zen_Idx_Rtm & , Sfc_Level_Rtm_Pixel & , Solar_Rtm & , Use_Sst_Anal & , Sst_Anal & , Trans_Atm_Ch20_Solar_Rtm & , Trans_Atm_Ch20_Solar_Total_Rtm & , Beta_11um_12um_Tropo_Rtm & , Beta_11um_104um_Tropo_Rtm & , Beta_11um_85um_Tropo_Rtm & , Beta_11um_67um_Tropo_Rtm & , Beta_11um_133um_Tropo_Rtm & , Pixel_Local_Time_Hours & , Ancil_Data_Dir & , Pc_Opaque_Cloud & , Zc_Opaque_Cloud & , Beta_104um_12um_Tropo_Rtm & , Beta_104um_85um_Tropo_Rtm & , Beta_104um_67um_Tropo_Rtm & , Beta_104um_133um_Tropo_Rtm & , Beta_104um_11um_Tropo_Rtm use NUMERICAL_ROUTINES_MOD , only: & LOCATE use CX_DATE_TIME_TOOLS_MOD , only: & COMPUTE_TIME_HOURS use CX_SCIENCE_TOOLS_MOD , only: & VAPOR use PLANCK_MOD, only: & PLANCK_RAD_FAST & , PLANCK_TEMP_FAST use CALIBRATION_CONSTANTS_MOD, only: & SOLAR_CH20_NU use RTM_COMMON_MOD, only: & Nlevels_Rtm & , Rtm & , P_Std_Rtm & , T_Std_Rtm & , Wvmr_Std_Rtm & , Ozmr_Std_Rtm use CX_PFAAST_MOD, only: & COMPUTE_TRANSMISSION_PFAAST use CLAVRX_MESSAGE_MOD implicit none private private:: EMISSIVITY, & NADIR_EMISSIVITY, & BETA_RATIO, & INVERSION_PROFILE public:: & SETUP_SOLAR_RTM, & MAP_NWP_RTM, & CREATE_TEMP_NWP_VECTORS, & DESTROY_TEMP_NWP_VECTORS, & GET_PIXEL_NWP_RTM, & ALLOCATE_RTM, & DEALLOCATE_RTM, & DEALLOCATE_RTM_VARS, & DEALLOCATE_RTM_CELL, & COMPUTE_CLEAR_SKY_SCATTER, & ATMOS_CORR integer, parameter, public:: Rtm_Nvzen = 50 real, parameter :: Co2_Ratio = 380.0 !in ppmv integer, parameter :: Chan_Idx_Min = 1 integer, parameter :: Chan_Idx_Max = 44 real(kind=real4), save, dimension(Chan_Idx_Min:Chan_Idx_Max):: Gamma_Trans_Factor real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Trans_Atm_Prof real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Trans_Atm_Solar_Prof real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Trans_Atm_Total_Prof real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Trans_Atm_Prof_Prev real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Trans_Atm_Solar_Prof_Prev real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Trans_Atm_Total_Prof_Prev real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Rad_Atm_Prof real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Rad_Atm_Dwn_Prof real, dimension(NLevels_Rtm,Chan_Idx_Min:Chan_Idx_Max), save:: Rad_BB_Cloud_Prof integer, parameter:: Ilon_Stride = 0 integer, parameter:: Ilat_Stride = 0 integer, parameter:: Ivza_Stride = 0 integer, dimension(NLevels_Rtm) :: k_Rtm_Nwp real, dimension(NLevels_Rtm) :: x_Rtm_Nwp real, dimension(NLevels_Rtm), save:: & T_Prof_Rtm, & Z_Prof_Rtm, & Wvmr_Prof_Rtm, & Ozmr_Prof_Rtm, & Tpw_Prof_Rtm, & Trans_Prof_Rtm integer, dimension(:), allocatable, save:: k_Nwp_Rtm real, dimension(:), allocatable, save:: x_Nwp_Rtm real, dimension(:), allocatable, save:: Wvmr_Nwp !----- PFAAST Specific Settings (all private to this module) integer, dimension(Chan_Idx_Min:Chan_Idx_Max), save:: Rtm_Chan_Idx character(len=20), dimension(Chan_Idx_Min:Chan_Idx_Max), save:: Pfaast_Name integer, save:: Sc_Id_Rtm character(len=20), save:: Sc_Name_Rtm real, parameter:: Rtm_Vza_Binsize = 0.02 contains !==================================================================== ! subroutine Name: CREATE_TEMP_NWP_VECTORS ! ! Function: ! create and initialize NWP vectors used for RTM calcs ! !==================================================================== subroutine CREATE_TEMP_NWP_VECTORS() allocate( Wvmr_Nwp(NWP%Nlevels), & k_Nwp_Rtm(NWP%Nlevels), & x_Nwp_Rtm(NWP%Nlevels)) Wvmr_Nwp = 0.0 k_Nwp_Rtm = 0.0 x_Nwp_Rtm = 0.0 end subroutine CREATE_TEMP_NWP_VECTORS !==================================================================== ! subroutine Name: DESTROY_TEMP_NWP_VECTORS ! ! Function: ! destroy and initialize NWP vectors used for RTM calcs ! !==================================================================== subroutine DESTROY_TEMP_NWP_VECTORS() deallocate(Wvmr_Nwp,k_Nwp_Rtm, x_Nwp_Rtm) end subroutine DESTROY_TEMP_NWP_VECTORS !==================================================================== ! subroutine Name: MAP_NWP_RTM ! ! Function: ! develops the mapping of the NWP and RTM profile Levels ! !==================================================================== subroutine MAP_NWP_RTM(Num_Levels_Nwp_Profile, & Press_Profile_Nwp, & Num_Levels_Rtm_Profile, & Press_Profile_Rtm) integer, intent(in):: Num_Levels_Nwp_Profile integer, intent(in):: Num_Levels_Rtm_Profile real, intent(in), dimension(:):: Press_Profile_Nwp real, intent(in), dimension(:):: Press_Profile_Rtm integer:: k integer:: k_temp !--- map NWP Levels to RTM Levels do k = 1, Num_Levels_Nwp_Profile !--- locate the nwp Level within the Rtm Levels !--- P_Std_Nwp(k) should fall between Rtm Levels k_temp and k_temp +1 call LOCATE(Press_Profile_Rtm, Num_Levels_Rtm_Profile, Press_Profile_Nwp(k), k_temp) k_Nwp_Rtm(k) = min(Num_Levels_Rtm_Profile-1,max(1,k_temp)) !-- compute linear weighting factor x_Nwp_Rtm(k) = (Press_Profile_Nwp(k) - Press_Profile_Rtm(k_Nwp_Rtm(k))) / & (Press_Profile_Rtm(k_Nwp_Rtm(k)+1) - Press_Profile_Rtm(k_Nwp_Rtm(k))) end do !--- map RTM Levels to NWP Levels do k = 1, Num_Levels_Rtm_Profile !--- locate the Rtm Level within the Nwp Levels !--- Press_Profile_Rtm(k) should fall between Nwp Levels k_temp and k_temp +1 call LOCATE(Press_Profile_Nwp, Num_Levels_Nwp_Profile, Press_Profile_Rtm(k), k_temp) k_Rtm_Nwp(k) = min(Num_Levels_Nwp_Profile-1,max(1,k_temp)) !-- compute linear weighting factor x_Rtm_Nwp(k) = (Press_Profile_Rtm(k) - Press_Profile_Nwp(k_Rtm_Nwp(k))) / & (Press_Profile_Nwp(k_Rtm_Nwp(k)+1) - Press_Profile_Nwp(k_Rtm_Nwp(k))) end do end subroutine MAP_NWP_RTM !==================================================================== ! subroutine Name: CONVERT_ATMOS_PROF_NWP_RTM ! ! Description: ! This routine interpolate the NWP profiles to profiles with the ! vertical spacing defined by the RTM model. It operates on profiles ! stored in this module ! ! INPUTS: ! ! Highest_Level_Rtm_Nwp - highest Rtm Level that is below the highest nwp Level ! Lowest_Level_Rtm_Nwp - lowest Rtm Level that is above the lowest nwp Level ! Sfc_Level_Rtm - lowest Rtm Level above the surface ! P_near_Sfc_Nwp - lowest standard nwp Level above surface pressure ! !==================================================================== subroutine CONVERT_ATMOS_PROF_NWP_RTM(Num_Levels_NWP_Profile, & Surface_Level_Nwp, & Surface_Elevation_Nwp, & Air_Temperature_Nwp, & Surface_Rh_Nwp, & Surface_Pressure_Nwp, & Press_Profile_Nwp, & T_Profile_Nwp, & Z_Profile_Nwp, & Wvmr_Profile_Nwp, & Ozmr_Profile_Nwp, & Num_Levels_Rtm_Profile, & Press_Profile_Rtm, & T_Profile_Rtm, & Z_Profile_Rtm, & Wvmr_Profile_Rtm, & Ozmr_Profile_Rtm, & T_Std_Profile_Rtm, & Wvmr_Std_Profile_Rtm, & Ozmr_Std_Profile_Rtm) integer, intent(in):: Num_Levels_Nwp_Profile integer(kind=int1), intent(in):: Surface_Level_Nwp real, intent(in):: Surface_Elevation_Nwp real, intent(in):: Air_Temperature_Nwp real, intent(in):: Surface_Rh_Nwp real, intent(in):: Surface_Pressure_Nwp real, intent(in), dimension(:):: Press_Profile_Nwp real, intent(in), dimension(:):: T_Profile_Nwp real, intent(in), dimension(:):: Z_Profile_Nwp real, intent(in), dimension(:):: Wvmr_Profile_Nwp real, intent(in), dimension(:):: Ozmr_Profile_Nwp integer, intent(in):: Num_Levels_Rtm_Profile real, intent(in), dimension(:):: Press_Profile_Rtm real, intent(out), dimension(:):: T_Profile_Rtm real, intent(out), dimension(:):: Z_Profile_Rtm real, intent(out), dimension(:):: Wvmr_Profile_Rtm real, intent(out), dimension(:):: Ozmr_Profile_Rtm real, intent(in), dimension(:):: T_Std_Profile_Rtm real, intent(in), dimension(:):: Wvmr_Std_Profile_Rtm real, intent(in), dimension(:):: Ozmr_Std_Profile_Rtm integer:: k integer:: Lowest_Level_Rtm_Nwp integer:: Highest_Level_Rtm_Nwp integer:: Sfc_Level_Rtm real:: dT_dP_near_Sfc real:: dWvmr_dP_near_Sfc real:: dZ_dP_near_Sfc real:: P_near_Sfc_Nwp real:: Wvmr_Sfc real:: es real:: e real:: T_Offset !--- initialize indices Lowest_Level_Rtm_Nwp = Num_Levels_Rtm_Profile Highest_Level_Rtm_Nwp = 1 Sfc_Level_Rtm = Num_Levels_Rtm_Profile P_Near_Sfc_Nwp = Press_Profile_Nwp(Surface_Level_Nwp) !--- make Wvmr at sfc es = VAPOR(Air_Temperature_Nwp) e = Surface_Rh_Nwp * es / 100.0 Wvmr_Sfc = 1000.0*0.622 * (e / (Surface_Pressure_Nwp - e)) !(g/kg) !--- determine some critical Levels in the Rtm profile do k = 1, Num_Levels_Rtm_Profile if (Press_Profile_Rtm(k) > Press_Profile_Nwp(1)) then Highest_Level_Rtm_Nwp = k exit endif enddo do k = Num_Levels_Rtm_Profile,1,-1 if (Press_Profile_Rtm(k) < Surface_Pressure_Nwp) then Sfc_Level_Rtm = k exit endif enddo do k = Num_Levels_Rtm_Profile,1,-1 if (Press_Profile_Rtm(k) < P_Near_Sfc_Nwp) then Lowest_Level_Rtm_Nwp = k exit endif enddo !--- compute T and Wvmr lapse rate near the surface dT_dP_near_Sfc = 0.0 dZ_dP_near_Sfc = 0.0 dWvmr_dP_near_Sfc = 0.0 if (Surface_Pressure_Nwp /= Press_Profile_Nwp(Num_Levels_Nwp_Profile)) then dT_dP_near_Sfc = & (T_Profile_Nwp(Surface_Level_Nwp) - Air_Temperature_Nwp)/ & (Press_Profile_Nwp(Surface_Level_Nwp) - Surface_Pressure_Nwp) dWvmr_dP_near_Sfc = & (Wvmr_Profile_Nwp(Surface_Level_Nwp) - Wvmr_Sfc)/ & (Press_Profile_Nwp(Surface_Level_Nwp) - Surface_Pressure_Nwp) else dT_dP_near_Sfc = & (T_Profile_Nwp(Surface_Level_Nwp-1) - Air_Temperature_Nwp)/ & (Press_Profile_Nwp(Surface_Level_Nwp-1) - Surface_Pressure_Nwp) dWvmr_dP_near_Sfc = & (Wvmr_Profile_Nwp(Surface_Level_Nwp-1) - Wvmr_Sfc)/ & (Press_Profile_Nwp(Surface_Level_Nwp-1) - Surface_Pressure_Nwp) end if if (Press_Profile_Nwp(Surface_Level_Nwp-1) /= Press_Profile_Nwp(Surface_Level_Nwp)) then dZ_dP_near_Sfc = & (Z_Profile_Nwp(Surface_Level_Nwp-1) - Z_Profile_Nwp(Surface_Level_Nwp))/ & (Press_Profile_Nwp(Surface_Level_Nwp-1) - Press_Profile_Nwp(Surface_Level_Nwp)) end if ! --- compute temperature offset between standard and NWP profiles at top !--- this will be added to the standard profiles T_Offset = T_Profile_Nwp(1) - T_Std_Profile_Rtm(Highest_Level_Rtm_Nwp) !--- for Rtm Levels above the highest nwp Level, use Rtm standard !values do k = 1,Highest_Level_Rtm_Nwp-1 Z_Profile_Rtm(k) = Z_Profile_Nwp(1) T_Profile_Rtm(k) = T_Std_Profile_Rtm(k) + T_Offset Wvmr_Profile_Rtm(k) = Wvmr_Std_Profile_Rtm(k) Ozmr_Profile_Rtm(k) = Ozmr_Std_Profile_Rtm(k) end do !--- Rtm Levels within standard nwp Levels above the surface do k = Highest_Level_Rtm_Nwp, Lowest_Level_Rtm_Nwp T_Profile_Rtm(k) = T_Profile_Nwp(k_Rtm_Nwp(k)) + x_Rtm_Nwp(k) * & (T_Profile_Nwp(k_Rtm_Nwp(k)+1) - T_Profile_Nwp(k_Rtm_Nwp(k))) Z_Profile_Rtm(k) = Z_Profile_Nwp(k_Rtm_Nwp(k)) + x_Rtm_Nwp(k) * & (Z_Profile_Nwp(k_Rtm_Nwp(k)+1) - Z_Profile_Nwp(k_Rtm_Nwp(k))) Wvmr_Profile_Rtm(k) = Wvmr_Profile_Nwp(k_Rtm_Nwp(k)) + x_Rtm_Nwp(k) * & (Wvmr_Profile_Nwp(k_Rtm_Nwp(k)+1) - Wvmr_Profile_Nwp(k_Rtm_Nwp(k))) Ozmr_Profile_Rtm(k) = Ozmr_Profile_Nwp(k_Rtm_Nwp(k)) + x_Rtm_Nwp(k) * & (Ozmr_Profile_Nwp(k_Rtm_Nwp(k)+1) - Ozmr_Profile_Nwp(k_Rtm_Nwp(k))) end do !--- Rtm Levels that are below the lowest nwp Level but above the surface do k = Lowest_Level_Rtm_Nwp+1, Sfc_Level_Rtm T_Profile_Rtm(k) = Air_Temperature_Nwp + dT_dP_near_Sfc * & (Press_Profile_Rtm(k) - Surface_Pressure_Nwp) Wvmr_Profile_Rtm(k) = Wvmr_Sfc + dWvmr_dP_near_Sfc * & (Press_Profile_Rtm(k) - Surface_Pressure_Nwp) Z_Profile_Rtm(k) = Surface_Elevation_Nwp + dZ_dP_near_Sfc * & (Press_Profile_Rtm(k) - Surface_Pressure_Nwp) Ozmr_Profile_Rtm(k) = Ozmr_Profile_Nwp(Num_Levels_Nwp_Profile) end do !--- Rtm Levels below the surface do k = Sfc_Level_Rtm +1, Num_Levels_Rtm_Profile T_Profile_Rtm(k) = Air_Temperature_Nwp Z_Profile_Rtm(k) = Surface_Elevation_Nwp + & dZ_dP_near_Sfc * (Press_Profile_Rtm(k) - Surface_Pressure_Nwp) Wvmr_Profile_Rtm(k) = Wvmr_Sfc Ozmr_Profile_Rtm(k) = Ozmr_Std_Profile_Rtm(k) end do !--- if using NCEP reanalysis which has no ozone profile, use default if (NWP_PIX%Nwp_Opt == 2) then Ozmr_Profile_Rtm = Ozmr_Std_Profile_Rtm end if end subroutine CONVERT_ATMOS_PROF_NWP_RTM !==================================================================== ! subroutine Name: COMPUTE_CLEAR_RAD_PROFILES_RTM ! ! Function: ! Computes clear sky radiance profiles ! !==================================================================== subroutine COMPUTE_CLEAR_RAD_PROFILES_RTM() integer:: Chan_Idx integer:: Lev_Idx real:: T_mean real:: B_mean real:: B_Level real:: Opd_Layer real:: Trans_Layer real:: Trans_Total !--- upwelling profiles Rad_Atm_Prof = Missing_Value_Real4 Rad_BB_Cloud_Prof = Missing_Value_Real4 do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (Ch(Chan_Idx) %Obs_Type /= THERMAL_OBS_TYPE .and. & Ch(Chan_Idx)%Obs_Type /= MIXED_OBS_TYPE) cycle if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%NO) cycle Rad_Atm_Prof(1,Chan_Idx) = 0.0 Rad_BB_Cloud_Prof(1,Chan_Idx) = 0.0 do Lev_Idx = 2, NLevels_Rtm T_mean = 0.5*(T_Prof_Rtm(Lev_Idx-1) + T_Prof_Rtm(Lev_Idx)) B_mean = PLANCK_RAD_FAST(Chan_Idx,T_mean) Rad_Atm_Prof(Lev_Idx,Chan_Idx) = Rad_Atm_Prof(Lev_Idx-1,Chan_Idx) + & (Trans_Atm_Prof(Lev_Idx-1,Chan_Idx) - Trans_Atm_Prof(Lev_Idx,Chan_Idx)) * B_mean B_Level = PLANCK_RAD_FAST(Chan_Idx,T_Prof_Rtm(Lev_Idx)) Rad_BB_Cloud_Prof(Lev_Idx,Chan_Idx) = Rad_Atm_Prof(Lev_Idx,Chan_Idx) + & (Trans_Atm_Prof(Lev_Idx,Chan_Idx) * B_Level) end do end do !--- downwelling profiles Rad_Atm_Dwn_Prof = Missing_Value_Real4 do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (Chan_Idx == 31 .OR. Chan_Idx == 38) then if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%NO) cycle Trans_Total = 1.0 Rad_Atm_Dwn_Prof(1,Chan_Idx) = 0.0 do Lev_Idx = 2, Nlevels_Rtm T_mean = 0.5*(T_Prof_Rtm(Lev_Idx) + T_Prof_Rtm(Lev_Idx-1)) B_mean = PLANCK_RAD_FAST(Chan_Idx,T_mean) Opd_Layer = -1.0*log(Trans_Atm_Prof(Lev_Idx,Chan_Idx)/Trans_Atm_Prof(Lev_Idx-1,Chan_Idx)) Opd_Layer = max(0.0,Opd_Layer) Trans_Layer = exp(-1.0*Opd_Layer) Rad_Atm_Dwn_Prof(Lev_Idx,Chan_Idx) = Trans_Total * Rad_Atm_Dwn_Prof(Lev_Idx-1,Chan_Idx) + & (1.0-Trans_Layer) * B_mean Trans_Total = Trans_Total * Trans_Layer end do else cycle endif end do end subroutine COMPUTE_CLEAR_RAD_PROFILES_RTM !==================================================================== ! subroutine Name: GET_PIXEL_NWP_RTM ! ! Function: ! Calculates the PFAAST transmittance profiles for a given segment. ! !==================================================================== subroutine GET_PIXEL_NWP_RTM(Line_Idx_Min,Num_Lines) integer, intent(in):: Line_Idx_Min integer, intent(in):: Num_Lines integer:: Elem_Idx integer:: Line_Idx integer:: Sfc_Level_Idx integer:: Lon_Idx integer:: Lat_Idx integer:: Lon_Idx_x integer:: Lat_Idx_x integer:: Lon_Idx_Prev integer:: Lat_Idx_Prev integer:: Zen_Idx_Prev real:: Lat_x real:: Lon_x real:: Prof_Weight real:: Satzen_Mid_Bin integer:: Zen_Idx integer:: Error_Status integer:: Chan_Idx integer:: Chan_Idx_For_Pfaast real(kind=real4),save :: Segment_Time_Point_Seconds_Temp = 0 real(kind=real4) :: Start_Time_Point_Hours_Temp real(kind=real4) :: End_Time_Point_Hours_Temp !-------------------------------------------------------------------------- ! Compute Gamma_Trans_Factor for radiance bias !-------------------------------------------------------------------------- call COMPUTE_GAMMA_FACTOR() !################### TEST ! Gamma_Trans_Factor = 0.7 !################### TEST Lon_Idx_Prev = 0 Lat_Idx_Prev = 0 Zen_Idx_Prev = 0 !--- loop over pixels in segment line_loop: do Line_Idx = Line_Idx_Min, Num_Lines + Line_Idx_Min - 1 element_loop: do Elem_Idx = 1, Image%Number_Of_Elements !--- check for bad scans if (Bad_Pixel_Mask(Elem_Idx,Line_Idx) == sym%YES) then cycle endif !--- check for space views if (Geo%Space_Mask(Elem_Idx,Line_Idx) == sym%YES) then cycle endif !--- compute viewing zenith bin for Rtm calculation Zen_Idx_Rtm(Elem_Idx,Line_Idx) = & max(1,min(Rtm_Nvzen,ceiling(Geo%Coszen(Elem_Idx,Line_Idx)/Rtm_Vza_Binsize))) !--- store cell and angular indices Lon_Idx = NWP_PIX%I_Nwp(Elem_Idx,Line_Idx) Lat_Idx = NWP_PIX%J_Nwp(Elem_Idx,Line_Idx) Lon_Idx_x = NWP_PIX%I_Nwp_X(Elem_Idx,Line_Idx) Lat_Idx_x = NWP_PIX%J_Nwp_X(Elem_Idx,Line_Idx) Lon_x = NWP_PIX%Lon_Nwp_Fac(Elem_Idx,Line_Idx) Lat_x = NWP_PIX%Lat_Nwp_Fac(Elem_Idx,Line_Idx) Zen_Idx = Zen_Idx_Rtm(Elem_Idx,Line_Idx) !--- if this is the first time this nwp cell is being processed do the following if (Rtm(Lon_Idx,Lat_Idx)%Flag == 0) then !--- allocate Rtm arrays call ALLOCATE_RTM_CELL(Lon_Idx,Lat_Idx,Rtm_Nvzen) !--- compute mixing ratio profile call COMPUTE_WVMR_PROFILE_NWP(NWP%P_Std, & NWP%T_Prof(:,Lon_Idx,Lat_Idx), & NWP%Rh_Prof(:,Lon_Idx,Lat_Idx), & Wvmr_Nwp) !--- compute tpw profiles call COMPUTE_TPW_PROFILE_NWP(NWP%P_Std, & Wvmr_Nwp, & NWP%Tpw_Prof(:,Lon_Idx,Lat_Idx)) !--- convert the atmospheric profiles from nwp to Rtm pressure coords call CONVERT_ATMOS_PROF_NWP_RTM(NWP%NLevels, & NWP%Sfc_Level(Lon_Idx,Lat_Idx), & NWP%Zsfc(Lon_Idx,Lat_Idx), & NWP%Tmpair(Lon_Idx,Lat_Idx), & NWP%Rhsfc(Lon_Idx,Lat_Idx), & NWP%Psfc(Lon_Idx,Lat_Idx), & NWP%P_Std, & NWP%T_Prof(:,Lon_Idx,Lat_Idx), & NWP%Z_Prof(:,Lon_Idx,Lat_Idx), & Wvmr_Nwp, & NWP%Ozone_Prof(:,Lon_Idx,Lat_Idx), & NLevels_Rtm, & P_Std_Rtm, & T_Prof_Rtm, & Z_Prof_Rtm, & Wvmr_Prof_Rtm, & Ozmr_Prof_Rtm, & T_Std_Rtm, & Wvmr_Std_Rtm, & Ozmr_Std_Rtm) !--- compute tpw profiles call COMPUTE_TPW_PROFILE_NWP(P_Std_Rtm, & Wvmr_Prof_Rtm, & Tpw_Prof_Rtm) !--- store in Rtm structures Rtm(Lon_Idx,Lat_Idx)%T_Prof = T_Prof_Rtm Rtm(Lon_Idx,Lat_Idx)%Z_Prof = Z_Prof_Rtm Rtm(Lon_Idx,Lat_Idx)%Wvmr_Prof = Wvmr_Prof_Rtm Rtm(Lon_Idx,Lat_Idx)%Ozmr_Prof = Ozmr_Prof_Rtm Rtm(Lon_Idx,Lat_Idx)%Tpw_Prof = Tpw_Prof_Rtm !--- find key Levels (sfc, tropo, Inversions) for future processing call FIND_RTM_LEVELS(Lon_Idx,Lat_Idx) !---- compute inversion level call INVERSION_PROFILE(Rtm(Lon_Idx,Lat_Idx)%T_Prof,Rtm(Lon_Idx,Lat_Idx)%Inver_Prof) !--- set flag so not done again Rtm(Lon_Idx,Lat_Idx)%Flag = 1 end if !--- determine if RTM needs to be called for this pixel, if not do the following if (Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%Flag == 0) then !check for existing RTM results call ALLOCATE_GLOBAL_RTM_STRUCTURE_ELEMENT(Lon_Idx,Lat_Idx,Zen_Idx) !--- determine the zenith angle for the RTM Satzen_Mid_Bin = acos((Zen_Idx-1)*Rtm_Vza_Binsize + Rtm_Vza_Binsize/2.0) / DTOR !--- store structures into RTM arguments T_Prof_Rtm = Rtm(Lon_Idx,Lat_Idx)%T_Prof Tpw_Prof_Rtm = Rtm(Lon_Idx,Lat_Idx)%Tpw_Prof Wvmr_Prof_Rtm = Rtm(Lon_Idx,Lat_Idx)%Wvmr_Prof Ozmr_Prof_Rtm = Rtm(Lon_Idx,Lat_Idx)%Ozmr_Prof !-------------------------------------------------------------- ! Call Sensor Specific Fast IR RTM !-------------------------------------------------------------- !--- decide if this can be skipped if ((Lon_Idx_Prev == 0) .or. (abs(Lon_Idx-Lon_Idx_Prev) >= Ilon_Stride) .or. & (Lat_Idx_Prev == 0) .or. (abs(Lat_Idx-Lat_Idx_Prev) >= Ilat_Stride) .or. & (Zen_Idx_Prev == 0) .or. (abs(Zen_Idx - Zen_Idx_Prev) >= Ivza_Stride)) then Start_Time_Point_Hours_Temp = COMPUTE_TIME_HOURS() !-------------------------------------------------------------- ! Call PFAAST Routines for IR channel Transmission Profiles !-------------------------------------------------------------- do Chan_Idx = Chan_Idx_Min,Chan_Idx_Max if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%NO) cycle if (ch(Chan_Idx)%Obs_Type == SOLAR_OBS_TYPE) cycle if (ch(Chan_Idx)%Obs_Type == LUNAR_OBS_TYPE) cycle Chan_Idx_For_Pfaast = Chan_Idx Sc_Name_Rtm = SENSOR_NAME_FOR_RTM(Sensor%WMO_id,Sensor%Sensor_Name, Chan_Idx) call COMPUTE_TRANSMISSION_PFAAST( & trim(Ancil_Data_Dir) & , T_Prof_rtm & , Wvmr_Prof_Rtm & , Ozmr_Prof_Rtm & , Satzen_Mid_Bin & , Sc_Name_Rtm & , Chan_Idx_For_Pfaast & , Trans_Prof_Rtm & , Use_Modis_Channel_Equivalent = .true. ) !---- Copy the output to appropriate channel's tranmission vector Trans_Atm_Prof(:,Chan_Idx) = Trans_Prof_Rtm end do !-------------------------------------------------------------- ! compute transmission profiles for solar channels !-------------------------------------------------------------- do Chan_Idx = Chan_Idx_Min,Chan_Idx_Max if (Chan_Idx >= 20 .and. Chan_Idx /= 26) cycle call SOLAR_TRANS(Chan_Idx,Satzen_Mid_Bin,Error_Status) !---- Copy the output to appropriate channel's tranmission vector Trans_Atm_Prof(:,Chan_Idx) = Trans_Prof_Rtm end do End_Time_Point_Hours_Temp = COMPUTE_TIME_HOURS() Segment_Time_Point_Seconds_Temp = Segment_Time_Point_Seconds_Temp + & 60.0*60.0*(End_Time_Point_Hours_Temp - Start_Time_Point_Hours_Temp) !---------------------------------------------------------------------- ! Manipulate Transmission Profiles !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! Apply Gamma Scaling !---------------------------------------------------------------------- do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%YES .and. Gamma_Trans_Factor(Chan_Idx)/=1.0) then Trans_Atm_Prof(:,Chan_Idx) = Trans_Atm_Prof(:,Chan_Idx) ** Gamma_Trans_Factor(Chan_Idx) end if end do !---------------------------------------------------------------------- ! Compute Tranmissions along solar and total (two-way) paths !---------------------------------------------------------------------- do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (Chan_Idx > 21 .and. Chan_Idx /= 26) cycle Trans_Atm_Total_Prof(:,Chan_Idx) = Trans_Atm_Prof(:,Chan_Idx) ** & ((Geo%Coszen(Elem_Idx,Line_Idx)+Geo%Cossolzen(Elem_Idx,Line_Idx))/Geo%Cossolzen(Elem_Idx,Line_Idx)) Trans_Atm_Solar_Prof(:,Chan_Idx) = Trans_Atm_Prof(:,Chan_Idx) ** & (Geo%Coszen(Elem_Idx,Line_Idx)/Geo%Cossolzen(Elem_Idx,Line_Idx)) end do Lon_Idx_Prev = Lon_Idx Lat_Idx_Prev = Lat_Idx Zen_Idx_Prev = Zen_Idx Trans_Atm_Prof_Prev = Trans_Atm_Prof Trans_Atm_Solar_Prof_Prev = Trans_Atm_Solar_Prof Trans_Atm_Total_Prof_Prev = Trans_Atm_Total_Prof else print *, "Skipped an RTM Call" Trans_Atm_Prof = Trans_Atm_Prof_Prev Trans_Atm_Solar_Prof = Trans_Atm_Solar_Prof_Prev Trans_Atm_Total_Prof = Trans_Atm_Total_Prof_Prev end if !skip rtm if statement !--- compute profiles of radiance (atm and bb cloud) call COMPUTE_CLEAR_RAD_PROFILES_RTM() !--- copy local rtm profiles back into global rtm structure call COPY_LOCAL_RTM_TO_GLOBAL_RTM_STRUCTURE(Lon_Idx,Lat_Idx,Zen_Idx) !--- set mask to indicate this bin or this cell has been computed Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%Flag = 1 end if !end check for previous RTM run !----------------------------------------------------------------------------- ! compute clear TOA radiance !----------------------------------------------------------------------------- !---- compute the surface level for this pixel based on high-res elevation Sfc_Level_Rtm_Pixel(Elem_Idx,Line_Idx) = Rtm(Lon_Idx,Lat_Idx)%Sfc_Level if ((Sfc%Land(Elem_Idx,Line_Idx) == sym%LAND) .and. (Sfc%Zsfc(Elem_Idx,Line_Idx) /= Missing_Value_Real4)) then call LOCATE(Rtm(Lon_Idx,Lat_Idx)%Z_Prof,NLevels_Rtm,Sfc%Zsfc(Elem_Idx,Line_Idx),Sfc_Level_Rtm_Pixel(Elem_Idx,Line_Idx)) endif Sfc_Level_Rtm_Pixel(Elem_Idx,Line_Idx) = max(1,min(NLevels_Rtm-1,Sfc_Level_Rtm_Pixel(Elem_Idx,Line_Idx))) !--- use this scalar for visual clarity Sfc_Level_Idx = Sfc_Level_Rtm_Pixel(Elem_Idx,Line_Idx) !--- if an sst analysis is available, use that if ((Sfc%Land_Mask(Elem_Idx,Line_Idx) == sym%NO) .and. & (Sfc%Snow(Elem_Idx,Line_Idx) == sym%NO_SNOW) .and. & (Use_Sst_Anal == sym%YES) .and. & (Sst_Anal(Elem_Idx,Line_Idx) > 270.0 )) then NWP_PIX%Tsfc(Elem_Idx,Line_Idx) = Sst_Anal(Elem_Idx,Line_Idx) end if !-- vertical interp weight Prof_Weight = (Sfc%Zsfc(Elem_Idx,Line_Idx) - Rtm(Lon_Idx,Lat_Idx)%Z_Prof(Sfc_Level_Idx)) / & (Rtm(Lon_Idx,Lat_Idx)%Z_Prof(Sfc_Level_Idx+1) - Rtm(Lon_Idx,Lat_Idx)%Z_Prof(Sfc_Level_Idx)) !--- constrain - important when high res topo differs from low res nwp topo Prof_Weight = max(0.0,min(1.0,Prof_Weight)) !--- call routine to compute radiative transfer terms such as Rad_Atm !--- Trans_Atm, and clear-sky radinace and brightness temperature !--- map global to local to allow for efficient looping call COMPUTE_CHANNEL_RT(Sfc_Level_Idx,Prof_Weight,Lon_Idx,Lat_Idx,Elem_Idx,Line_Idx,Zen_Idx) !--- compute Ch20 Emissivities call COMPUTE_CH20_EMISSIVITY(Elem_Idx,Line_Idx) !--- compute Emissivity at tropopause call COMPUTE_TROPOPAUSE_EMISSIVITIES(Elem_Idx,Line_Idx,Lon_Idx,Lat_Idx,Zen_Idx) !--- compute split-window beta ratio at tropopause call COMPUTE_BETA_RATIOES(Elem_Idx,Line_Idx) end do element_loop end do line_loop return end subroutine GET_PIXEL_NWP_RTM !==================================================================== ! subroutine Name: FIND_RTM_LEVELS ! ! Function: ! Finds various key Levels in the RTM (tropopause, etc), for each RTM gridcell ! !==================================================================== subroutine FIND_RTM_LEVELS(Lon_Idx,Lat_Idx) integer (kind=int4), intent(in):: Lon_Idx, Lat_Idx integer (kind=int4):: k !--- find surface Level - this is closest nwp Level above the actual surface Rtm(Lon_Idx,Lat_Idx)%Sfc_Level = NLevels_Rtm do k = NLevels_Rtm,1,-1 if (P_Std_Rtm(k) < NWP%Psfc(Lon_Idx,Lat_Idx)) then Rtm(Lon_Idx,Lat_Idx)%Sfc_Level = k exit end if end do !-------------------------------------------------------------------- !--- find tropopause Level based on tropopause pressure !--- tropopause is between tropopause_Level and tropopaue_Level + 1 !-------------------------------------------------------------------- do k = 1, Rtm(Lon_Idx,Lat_Idx)%Sfc_Level-1 if ((P_Std_Rtm(k) <= NWP%P_Trop(Lon_Idx,Lat_Idx)) .and. & (P_Std_Rtm(k+1) > NWP%P_Trop(Lon_Idx,Lat_Idx))) then Rtm(Lon_Idx,Lat_Idx)%Tropo_Level = k end if end do !--- check if tropopause Level found if (Rtm(Lon_Idx,Lat_Idx)%Tropo_Level == 0) then print *, EXE_PROMPT, "Error, tropopause Level not found" end if do k = 1, Rtm(Lon_Idx,Lat_Idx)%Sfc_Level-1 if ((P_Std_Rtm(k) <= 850.0) .and. & (P_Std_Rtm(k+1) > 850.0)) then Rtm(Lon_Idx,Lat_Idx)%Level850 = k endif enddo do k = 1, Rtm(Lon_Idx,Lat_Idx)%Sfc_Level-1 if ((P_Std_Rtm(k) <= 440.0) .and. & (P_Std_Rtm(k+1) > 440.0)) then Rtm(Lon_Idx,Lat_Idx)%Level440 = k endif enddo !--------------------------------------------------------------------- ! find Inversion Level - highest Level Inversion below trop !--------------------------------------------------------------------- Rtm(Lon_Idx,Lat_Idx)%Inversion_Level = 0 if (Rtm(Lon_Idx,Lat_Idx)%Tropo_Level > 0 .and. Rtm(Lon_Idx,Lat_Idx)%Sfc_Level > 0) then do k = Rtm(Lon_Idx,Lat_Idx)%Tropo_Level,Rtm(Lon_Idx,Lat_Idx)%Sfc_Level-1 if ((Rtm(Lon_Idx,Lat_Idx)%T_Prof(k) - Rtm(Lon_Idx,Lat_Idx)%T_Prof(k+1) > delta_t_Inversion) .and. & (P_Std_Rtm(k) >= p_Inversion_min)) then Rtm(Lon_Idx,Lat_Idx)%Inversion_Level = k exit endif enddo endif end subroutine FIND_RTM_LEVELS !==================================================================== ! subroutine Name: COMPUTE_WVMR_PROFILE_NWP ! ! Function: ! Computes the NWP WVMR profile ! ! Input: Lon_Idx = Longitude Index for NWP cell ! Lat_Idx = Latitude Index for NWP cell ! Press_Profile = pressure profile (hPa) ! Temp_Profile = temperature profile (K) ! Rh_Profile = relative humidity profile (%) ! ! Output: Wvmr_Profile = water vapor mixing ratio profile (g/kg) ! Tpw Profile = profile of Tpw from level to space (g/m^2) ! !==================================================================== subroutine COMPUTE_WVMR_PROFILE_NWP(Press_Profile, & Temp_Profile, & Rh_Profile, & Wvmr_Profile) real, intent(in), dimension(:):: Temp_Profile real, intent(in), dimension(:):: Press_Profile real, intent(in), dimension(:):: Rh_Profile real, intent(out), dimension(:):: Wvmr_Profile integer:: Lev_Idx real:: e real:: es integer:: Nlevels Nlevels = size(Press_Profile) !--- make Wvmr_Profile for use in RTM do Lev_Idx = 1, NLevels es = VAPOR(Temp_Profile(Lev_Idx)) e = Rh_Profile(Lev_Idx) * es / 100.0 Wvmr_Profile(Lev_Idx) = 1000.0*0.622 * (e / (Press_Profile(Lev_Idx) - e)) !(g/kg) end do end subroutine COMPUTE_WVMR_PROFILE_NWP !==================================================================== ! subroutine Name: COMPUTE_TPW_PROFILE_NWP ! ! Function: ! Computes the NWP TPW profile ! ! Input: Lon_Idx = Longitude Index for NWP cell ! Lat_Idx = Latitude Index for NWP cell ! Press_Profile = pressure profile (hPa) ! Temp_Profile = temperature profile (K) ! Rh_Profile = relative humidity profile (%) ! ! Output: Wvmr_Profile = water vapor mixing ratio profile (g/kg) ! Tpw Profile = profile of Tpw from level to space (g/m^2) ! !==================================================================== subroutine COMPUTE_TPW_PROFILE_NWP(Press_Profile, & Wvmr_Profile, & Tpw_Profile) real, intent(in), dimension(:):: Press_Profile real, intent(in), dimension(:):: Wvmr_Profile real, intent(out), dimension(:):: Tpw_Profile integer:: Lay_Idx real :: w_mean real :: u_layer integer:: Nlevels Nlevels = size(Press_Profile) !--- make tpw profile for use in atmospheric correction Tpw_Profile(1) = 0.0 do Lay_Idx = 1, NLevels-1 !layer index w_mean = 0.5*(Wvmr_Profile(Lay_Idx+1)+Wvmr_Profile(Lay_Idx)) / 1000.0 !(kg/kg) u_layer = (10.0/g)*(Press_Profile(Lay_Idx+1)-Press_Profile(Lay_Idx))*w_mean Tpw_Profile(Lay_Idx+1) = Tpw_Profile(Lay_Idx) + u_layer end do end subroutine COMPUTE_TPW_PROFILE_NWP !==================================================================== ! subroutine Name: ALLOCATE_RTM ! ! Function: ! Subroutine that allocates memory for nlon x nlat Rtm structure. ! !==================================================================== subroutine ALLOCATE_RTM(nx,ny) integer (kind=int4), intent(in) :: nx, ny integer :: Alloc_Status allocate(Rtm(nx,ny),stat=Alloc_Status) if (Alloc_Status /= 0) then print "(a,'Not enough memory to allocate Rtm_Params structure.')",EXE_PROMPT stop end if end subroutine ALLOCATE_RTM !==================================================================== ! subroutine Name: DEALLOCATE_RTM ! ! Function: ! Subroutine that deallocates memory used for nlon x nlat Rtm structure. ! !==================================================================== subroutine DEALLOCATE_RTM() integer :: Alloc_Status deallocate(Rtm,stat=Alloc_Status) if (Alloc_Status /= 0) then print "(a,'Error deallocating Rtm_Params structure.')",EXE_PROMPT stop end if end subroutine DEALLOCATE_RTM !==================================================================== ! subroutine Name: ALLOCATE_RTM_CELL ! ! Function: ! Subroutine to allocate memory for the RTM structure. ! !==================================================================== subroutine ALLOCATE_RTM_CELL(Lon_Idx,Lat_Idx,Nvza) integer (kind=int4), intent(in) :: Lon_Idx, Lat_Idx integer (kind=int4), intent(in) :: Nvza integer :: Alloc_Status allocate(Rtm(Lon_Idx, Lat_Idx)%d(Nvza),stat=Alloc_Status) allocate(Rtm(Lon_Idx, Lat_Idx)%T_Prof(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx, Lat_Idx)%Z_Prof(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx, Lat_Idx)%Wvmr_Prof(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx, Lat_Idx)%Ozmr_Prof(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx, Lat_Idx)%Tpw_Prof(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx, Lat_Idx)%Inver_Prof(NLevels_Rtm),stat=Alloc_Status) if (Alloc_Status /= 0) then print "(a,'Not enough memory to allocate Rtm_Params structure.')",EXE_PROMPT stop endif Rtm(Lon_Idx,Lat_Idx)%Flag = 1 end subroutine ALLOCATE_RTM_CELL !==================================================================== ! subroutine Name: DEALLOCATE_RTM_CELL ! ! Function: ! Subroutine to deallocates memory for the RTM structure. ! !==================================================================== subroutine DEALLOCATE_RTM_CELL(Lon_Idx,Lat_Idx) integer (kind=int4), intent(in) :: Lon_Idx integer (kind=int4), intent(in) :: Lat_idx integer :: Alloc_Status deallocate(Rtm(Lon_Idx,Lat_Idx)%d,stat=Alloc_Status) deallocate(Rtm(Lon_Idx,Lat_Idx)%T_Prof,stat=Alloc_Status) deallocate(Rtm(Lon_Idx,Lat_Idx)%Z_Prof,stat=Alloc_Status) deallocate(Rtm(Lon_Idx,Lat_Idx)%Wvmr_Prof,stat=Alloc_Status) deallocate(Rtm(Lon_Idx,Lat_Idx)%Ozmr_Prof,stat=Alloc_Status) deallocate(Rtm(Lon_Idx,Lat_Idx)%Tpw_Prof,stat=Alloc_Status) deallocate(Rtm(Lon_Idx,Lat_Idx)%Inver_Prof,stat=Alloc_Status) if (Alloc_Status /= 0) then print "(a,'Error deallocating Rtm_Params structure.')",EXE_PROMPT stop endif Rtm(Lon_Idx,Lat_Idx) % Flag = 0 Rtm(Lon_Idx,Lat_Idx) % Sfc_Level = 0 Rtm(Lon_Idx,Lat_Idx) % Inversion_Level = 0 Rtm(Lon_Idx,Lat_Idx) % Tropo_Level = 0 Rtm(Lon_Idx,Lat_Idx) % Level440 = 0 Rtm(Lon_Idx,Lat_Idx) % Level850 = 0 end subroutine DEALLOCATE_RTM_CELL !==================================================================== ! subroutine Name: DEALLOCATE_RTM_VARS ! ! Function: ! Subroutine to deallocate memory for the RTM profile arrays. ! !==================================================================== subroutine DEALLOCATE_RTM_VARS(Lon_Idx,Lat_Idx,Zen_Idx) integer (kind=int4), intent(in) :: Lon_Idx integer (kind=int4), intent(in) :: Lat_Idx integer (kind=int4), intent(in) :: Zen_Idx integer :: Chan_Idx integer :: Alloc_Status do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (allocated(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Profile)) & deallocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Profile, stat=Alloc_Status) if (allocated(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Solar_Profile)) & deallocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Solar_Profile, stat=Alloc_Status) if (allocated(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Total_Profile)) & deallocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Total_Profile, stat=Alloc_Status) if (allocated(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Profile)) & deallocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Profile, stat=Alloc_Status) if (allocated(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Dwn_Profile)) & deallocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Dwn_Profile, stat=Alloc_Status) if (allocated(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_BB_Cloud_Profile)) & deallocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_BB_Cloud_Profile, stat=Alloc_Status) if (Alloc_Status /= 0) then print "(a,'Error destroying allocate channel Rtm profile.')",EXE_PROMPT stop endif enddo !--- set flag to zero (so we don't deallocate again) Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%Flag = 0 end subroutine DEALLOCATE_RTM_VARS !-------------------------------------------------------------- ! Compute Gamma Factor for Radiance Bias Adjustment ! ! based on satellite number, channel number and nwp source, ! set the value of Gamma_Trans_Factor to reduce RTM bias ! ! nwp flag (1=gfs,2=ncep reanalysis,3=cfsr) !-------------------------------------------------------------- subroutine COMPUTE_GAMMA_FACTOR() !--- initialize to unity Gamma_Trans_Factor = 1.0 !--- GOES-10 if (Sensor%WMO_Id == 254) then if (NWP_PIX%Nwp_Opt == 3) then Gamma_Trans_Factor(20) = 1.25 Gamma_Trans_Factor(27) = 0.79 Gamma_Trans_Factor(31) = 1.15 Gamma_Trans_Factor(31) = 1.05 endif endif !--- GOES-11 if (Sensor%WMO_Id == 255) then if (NWP_PIX%Nwp_opt == 3) then Gamma_Trans_Factor(20) = 1.35 Gamma_Trans_Factor(27) = 0.74 Gamma_Trans_Factor(31) = 1.05 Gamma_Trans_Factor(32) = 1.05 endif endif !--- GOES-12 if (Sensor%WMO_Id == 256) then if (NWP_PIX%Nwp_Opt == 1 .or. NWP_PIX%Nwp_Opt == 3) then !repeat of cfsr Gamma_Trans_Factor(20) = 1.45 Gamma_Trans_Factor(27) = 0.79 Gamma_Trans_Factor(31) = 1.15 Gamma_Trans_Factor(33) = 1.15 end if end if !--- GOES-13 if (Sensor%WMO_Id == 257) then if (NWP_PIX%Nwp_Opt == 1 .or. NWP_PIX%Nwp_Opt == 3) then Gamma_Trans_Factor(20) = 1.00 Gamma_Trans_Factor(27) = 0.794 Gamma_Trans_Factor(31) = 1.075 Gamma_Trans_Factor(33) = 1.064 end if end if !--- GOES-15 if (Sensor%WMO_Id == 259) then Gamma_Trans_Factor(20) = 1.55 Gamma_Trans_Factor(27) = 0.728 Gamma_Trans_Factor(31) = 1.05 Gamma_Trans_Factor(33) = 1.075 end if !--- MET-09 if (Sensor%WMO_Id == 56) then Gamma_Trans_Factor(20) = 1.55 Gamma_Trans_Factor(27) = 0.96 Gamma_Trans_Factor(29) = 1.35 Gamma_Trans_Factor(31) = 0.95 Gamma_Trans_Factor(32) = 0.95 Gamma_Trans_Factor(33) = 1.11 endif !--- MTSAT-02 if (Sensor%WMO_Id == 172) then if (NWP_PIX%Nwp_Opt == 1) then Gamma_Trans_Factor(20) = 1.25 Gamma_Trans_Factor(27) = 0.87 Gamma_Trans_Factor(31) = 1.15 Gamma_Trans_Factor(32) = 1.15 end if end if !--- COMS-1 if (Sensor%WMO_Id == 810) then Gamma_Trans_Factor(20) = 1.25 Gamma_Trans_Factor(27) = 0.936 Gamma_Trans_Factor(31) = 1.05 Gamma_Trans_Factor(32) = 1.075 end if !--- NOAA-19 HIRS if (Sensor%WMO_Id == 223) then Gamma_Trans_Factor(27) = 0.9083 Gamma_Trans_Factor(28) = 0.9500 Gamma_Trans_Factor(33) = 1.1250 Gamma_Trans_Factor(34) = 1.0438 Gamma_Trans_Factor(35) = 1.0833 Gamma_Trans_Factor(36) = 1.0917 end if end subroutine COMPUTE_GAMMA_FACTOR !-------------------------------------------------------------------------------------------------- !> subroutine NAME: SENSOR_NAME_FOR_RTM !! !! Description: !! Knowing the WMO Satellite Identification Number !! !-------------------------------------------------------------------------------------------------- function SENSOR_NAME_FOR_RTM ( wmo_id, sensorname, Chan_Idx ) result ( Sensor_Name_Rtm) integer, intent(in) :: wmo_id character (len =*) , intent(in) :: sensorname integer, intent(in) :: Chan_Idx character (len =20 ) :: Sensor_Name_Rtm integer :: i select case(WMO_Id) case(4) !METOP-A Sensor_Name_Rtm = 'AVHRR-METOPA' case(3) !METOP-B Sensor_Name_Rtm = 'AVHRR-METOPB' case(5) !METOP-C Sensor_Name_Rtm = 'AVHRR-METOPC' case(55) !MSG-8 Sensor_Name_Rtm = 'SEVIRI-MSG08' case(56) !MSG-9 Sensor_Name_Rtm = 'SEVIRI-MSG09' case(57) !MSG-10 Sensor_Name_Rtm = 'SEVIRI-MSG10' case(70) !MSG-11 Sensor_Name_Rtm = 'SEVIRI-MSG11' case(171) !MTSAT-1R Sensor_Name_Rtm = 'MTSAT-1' case(172) !MTSAT-2 Sensor_Name_Rtm = 'MTSAT-2' case(173) !AHI-8 Sensor_Name_Rtm = 'AHI8' case(174) !AHI-9 Sensor_Name_Rtm = 'AHI9' case(200) !NOAA-8 Sensor_Name_Rtm = 'AVHRR-NOAA08' case(201) !NOAA-9 Sensor_Name_Rtm = 'AVHRR-NOAA09' case(202) !NOAA-10 Sensor_Name_Rtm = 'AVHRR-NOAA10' case(203) !NOAA-11 Sensor_Name_Rtm = 'AVHRR-NOAA11' case(204) !NOAA-12 Sensor_Name_Rtm = 'AVHRR-NOAA12' case(205) !NOAA-14 Sensor_Name_Rtm = 'AVHRR-NOAA14' case(206) !NOAA-15 Sensor_Name_Rtm = 'AVHRR-NOAA15' case(207) !NOAA-16 Sensor_Name_Rtm = 'AVHRR-NOAA16' case(208) !NOAA-17 Sensor_Name_Rtm = 'AVHRR-NOAA17' case(209) !NOAA-18 Sensor_Name_Rtm = 'AVHRR-NOAA18' case(223) !NOAA-19 Sensor_Name_Rtm = 'AVHRR-NOAA19' case(224) !VIIRS - SNPP Sensor_Name_Rtm = 'VIIRS-SNPP' case(225) !VIIRS NOAA-20 Sensor_Name_Rtm = 'VIIRS-N20' case(252) !GOES-8 Sensor_Name_Rtm = 'GOES-8' case(253) !GOES-9 Sensor_Name_Rtm = 'GOES-9' case(254) !GOES-10 Sensor_Name_Rtm = 'GOES-10' case(255) !GOES-11 Sensor_Name_Rtm = 'GOES-11' case(256) !GOES-12 Sensor_Name_Rtm = 'GOES-12' case(257) !GOES-13 Sensor_Name_Rtm = 'GOES-13' case(258) !GOES-14 Sensor_Name_Rtm = 'GOES-14' case(259) !GOES-15 Sensor_Name_Rtm = 'GOES-15' case(270) !GOES-16 Sensor_Name_Rtm = 'GOES-16' case(271) !GOES-17 Sensor_Name_Rtm = 'GOES-17' case(706) !NOAA-6 Sensor_Name_Rtm = 'AVHRR-NOAA06' case(707) !NOAA-7 Sensor_Name_Rtm = 'AVHRR-NOAA07' case(708) !NOAA-5 Sensor_Name_Rtm = 'AVHRR-TIROSN' case(783) !MODIS Sensor_Name_Rtm = 'MODIS-TERRA' case(784) !MODIS Sensor_Name_Rtm = 'MODIS-AQUA' case(510) !FY2A Sensor_Name_Rtm ='FY2-1' case(514) !FY2D Sensor_Name_Rtm ='FY2-2' case(515) !FY2E Sensor_Name_Rtm ='FY2-3' case(523) !FY3D Sensor_Name_Rtm = 'FY3-D' case(530) ! FY4-A Sensor_Name_Rtm ='FY4-A' case(810) !COMS Sensor_Name_Rtm ='COMS-1' case default print*,'sensor for WMO number not found in RT Utils ', WMO_id print*,'stopping ... Please fix this in rt_utils.F90' print*,' better tell andi.walther@ssec.wisc.edu' stop end select if (trim ( Sensorname) == 'AVHRR-IFF' .or. & trim ( Sensorname) == 'AVHRR-FUSION') then ! sensor for channels 21:30 and 33:36 is HIRS if ( any ( Chan_Idx == [ (i,i=21,30,1) , 33,34,35,36] ) ) then ! - for this IFF Sensor_Name_Rtm is initially set to AVHRR- select case(WMO_Id) case(4) !METOP-A Sensor_Name_Rtm = 'HIRS-METOPA' case(3) !METOP-B Sensor_Name_Rtm = 'HIRS-METOPB' case(5) !METOP-C Sensor_Name_Rtm = 'HIRS-METOPC' case(200) !NOAA-8 Sensor_Name_Rtm = 'HIRS-NOAA08' case(201) !NOAA-9 Sensor_Name_Rtm = 'HIRS-NOAA09' case(202) !NOAA-10 Sensor_Name_Rtm = 'HIRS-NOAA10' case(203) !NOAA-11 Sensor_Name_Rtm = 'HIRS-NOAA11' case(204) !NOAA-12 Sensor_Name_Rtm = 'HIRS-NOAA12' case(205) !NOAA-14 Sensor_Name_Rtm = 'HIRS-NOAA14' case(206) !NOAA-15 Sensor_Name_Rtm = 'HIRS-NOAA15' case(207) !NOAA-16 Sensor_Name_Rtm = 'HIRS-NOAA16' case(208) !NOAA-17 Sensor_Name_Rtm = 'HIRS-NOAA17' case(209) !NOAA-18 Sensor_Name_Rtm = 'HIRS-NOAA18' case(223) !NOAA-19 Sensor_Name_Rtm = 'HIRS-NOAA19' case(706) !NOAA-6 Sensor_Name_Rtm = 'HIRS-NOAA06' case(707) !NOAA-7 Sensor_Name_Rtm = 'HIRS-NOAA07' case(708) !NOAA-5 Sensor_Name_Rtm = 'HIRS-TIROSN' case default print*,'sensor for WMO number not found in RT Utils for AVHRR-IFF ', WMO_id print*,'stopping ... Please fix this in rt_utils.F90' print*,' better tell andi.walther@ssec.wisc.edu' stop end select end if end if if (trim ( Sensorname) == 'VIIRS-IFF') then ! sensor for channels 27:28 and 33:36 is CRISP this is similar to MODIS-AQUA if ( any ( Chan_Idx == [27,28, 33,34,35,36] ) ) Sensor_Name_Rtm = 'MODIS-AQUA' end if if ( trim (sensorname ) == 'VIIRS-NASA' ) then ! - check what is with 31,32 if ( any ( Chan_Idx == [23,24,25,27,28,30,33,34,35,36] ) ) Sensor_Name_Rtm = 'MODIS-AQUA' end if end function SENSOR_NAME_FOR_RTM !-------------------------------------------------------------------------------------------------- ! set up values for the solar rtm calculations for this sensor !-------------------------------------------------------------------------------------------------- subroutine SOLAR_TRANS(Chan_Idx,Zen_Ang,Error_Status) integer, intent(in):: Chan_Idx real, intent(in):: Zen_Ang integer, intent(out):: Error_Status real, dimension(3):: Tau_H2O_Coef real:: Tau_O2_Column real:: Tau_CO2_Column real:: Tau_CH4_Column real:: Tau_O3_Column real:: Tau_O2 real:: Tau_CO2 real:: Tau_CH4 real:: Tau_O3 real:: Tau_H2O real:: Tpw real:: mu real:: Tau_Gas integer:: Lev_Idx Error_Status = 1 if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%NO) return ! if (Chan_Idx >= 20 .and. Chan_Idx /= 26 .and. Chan_Idx/= 44) return if (ch(Chan_Idx)%Obs_Type /= SOLAR_OBS_TYPE .or. ch(Chan_Idx)%Obs_Type /= LUNAR_OBS_TYPE) return Trans_Prof_Rtm = 1.0 mu = cos(Zen_Ang*DTOR) !--- initialize Tau_H2O_Coef = Solar_Rtm%Tau_H2O_Coef(Chan_Idx,:) Tau_O2_Column = Solar_Rtm%Tau_O2(Chan_Idx) Tau_O3_Column = Solar_Rtm%Tau_O3(Chan_Idx) Tau_CO2_Column = Solar_Rtm%Tau_CO2(Chan_Idx) Tau_CH4_Column = Solar_Rtm%Tau_CH4(Chan_Idx) !--- loop through layers and fill profile do Lev_Idx = 1, Nlevels_Rtm Tpw = Tpw_Prof_Rtm(Lev_Idx) Tau_H2O = Tau_H2O_Coef(1) + Tau_H2O_Coef(2)*Tpw_Prof_Rtm(Lev_Idx) + Tau_H2O_Coef(3)*Tpw_Prof_Rtm(Lev_Idx)**2 Tau_CO2 = Tau_Co2_Column * P_Std_Rtm(Lev_Idx) / P_Std_Rtm(Nlevels_Rtm) Tau_Ch4 = Tau_Ch4_Column * P_Std_Rtm(Lev_Idx) / P_Std_Rtm(Nlevels_Rtm) Tau_O2 = Tau_O2_Column * P_Std_Rtm(Lev_Idx) / P_Std_Rtm(Nlevels_Rtm) Tau_O3 = Tau_O3_Column * P_Std_Rtm(Lev_Idx) / P_Std_Rtm(Nlevels_Rtm) Tau_Gas = max(0.0,Tau_H2O + Tau_O3 + Tau_O2 + Tau_CO2 + Tau_CH4) Trans_Prof_Rtm(Lev_Idx) = exp(-Tau_Gas / mu) end do Error_Status = 0 end subroutine SOLAR_TRANS !============================================================================== ! !============================================================================== subroutine SETUP_SOLAR_RTM ( WMO_Id ) integer, intent(in):: WMO_Id real, parameter:: alpha = 0.5 real, parameter:: Tau_Aer_1um = 0.0 !0.10 !------------------------------------------------------------------------ ! initialize gas coefficients with default values !------------------------------------------------------------------------ include 'default_solar_terms.inc' !---- define channel aerosol properties Solar_Rtm%Tau_Aer = 0.0 Solar_Rtm%Wo_Aer = 0.8 Solar_Rtm%G_Aer = 0.6 ! tau(lambda) = tau(lambda=1) * lambda**(-alpha) Solar_Rtm%Tau_Aer(1) = Tau_Aer_1um * (0.65**(-1.0*alpha)) Solar_Rtm%Tau_Aer(2) = Tau_Aer_1um * (0.86**(-1.0*alpha)) Solar_Rtm%Tau_Aer(3) = Tau_Aer_1um * (0.47**(-1.0*alpha)) Solar_Rtm%Tau_Aer(5) = Tau_Aer_1um * (0.55**(-1.0*alpha)) Solar_Rtm%Tau_Aer(6) = Tau_Aer_1um * (1.60**(-1.0*alpha)) Solar_Rtm%Tau_Aer(7) = Tau_Aer_1um * (2.15**(-1.0*alpha)) Solar_Rtm%Tau_Aer(8) = Tau_Aer_1um * (0.41**(-1.0*alpha)) Solar_Rtm%Tau_Aer(9) = Tau_Aer_1um * (0.44**(-1.0*alpha)) Solar_Rtm%Tau_Aer(10) = Tau_Aer_1um * (0.49**(-1.0*alpha)) Solar_Rtm%Tau_Aer(11) = Tau_Aer_1um * (0.53**(-1.0*alpha)) Solar_Rtm%Tau_Aer(12) = Tau_Aer_1um * (0.55**(-1.0*alpha)) Solar_Rtm%Tau_Aer(13) = Tau_Aer_1um * (0.67**(-1.0*alpha)) Solar_Rtm%Tau_Aer(14) = Tau_Aer_1um * (0.68**(-1.0*alpha)) Solar_Rtm%Tau_Aer(15) = Tau_Aer_1um * (0.75**(-1.0*alpha)) Solar_Rtm%Tau_Aer(16) = Tau_Aer_1um * (0.87**(-1.0*alpha)) Solar_Rtm%Tau_Aer(17) = Tau_Aer_1um * (0.90**(-1.0*alpha)) Solar_Rtm%Tau_Aer(18) = Tau_Aer_1um * (0.93**(-1.0*alpha)) Solar_Rtm%Tau_Aer(19) = Tau_Aer_1um * (0.94**(-1.0*alpha)) Solar_Rtm%Tau_Aer(20) = Tau_Aer_1um * (3.75**(-1.0*alpha)) Solar_Rtm%Tau_Aer(26) = Tau_Aer_1um * (1.38**(-1.0*alpha)) select case(WMO_Id) include 'avhrr_solar_terms.inc' include 'goes_ip_solar_terms.inc' include 'viirs_solar_terms.inc' include 'seviri_solar_terms.inc' include 'modis_solar_terms.inc' include 'abi_solar_terms.inc' include 'ahi_solar_terms.inc' include 'coms_solar_terms.inc' include 'mtsat_solar_terms.inc' ! case(174) ! AHI-9 VALUES NEED TO BE UPDATED case default call MESG("WARNING: Solar transmission terms are not available, using default", level = verb_lev % DEFAULT) end select end subroutine SETUP_SOLAR_RTM !------------------------------------------------------------------------------ ! Routine to compute some needed radiative transfer terms for the IR channels ! ! Input: Chan_Idx - number of the channel being used ! Sfc_Idx - level of the surface in the profiles ! Profile_Weight - interpolation weight for estimated the surface in the ! profiles ! Sfc_Emiss - emissivity of the surface for this channel ! Sfc_Temp - the temperature of the surface ! Rad_Atm_Profile - profile of radiance emitted from level to space ! Trans_Atm_Profile - profile of transmissio from level to space ! ! Output: Rad_Atm - total radiance due to atmospheric emission ! Trans_Atm - total tranmission due to atmosphere ! Rad_Atm_Sfc - total radiance due both atmosphere and surface at TOA ! Bt_Atm_Sfc - Rad_Atm_Sfc expressed as a brightness temperature !------------------------------------------------------------------------------ subroutine COMPUTE_CHANNEL_ATM_SFC_RAD_BT( & Chan_Idx, & Sfc_Idx, & Profile_Weight, & Sfc_Emiss, & Sfc_Temp, & Rad_Atm_Profile, & Trans_Atm_Profile, & Rad_Atm, & Trans_Atm, & Rad_Atm_Sfc, & Bt_Atm_Sfc) integer, intent(in):: Chan_Idx integer, intent(in):: Sfc_Idx real, intent(in):: Profile_Weight real, intent(in):: Sfc_Emiss real, intent(in):: Sfc_Temp real, intent(in), dimension(:):: Rad_Atm_Profile real, intent(in), dimension(:):: Trans_Atm_Profile real, intent(out):: Rad_Atm real, intent(out):: Trans_Atm real, intent(out):: Rad_Atm_Sfc real, intent(out):: Bt_Atm_Sfc real:: Sfc_Rad Sfc_Rad = Sfc_Emiss * PLANCK_RAD_FAST(Chan_Idx,Sfc_Temp) Rad_Atm = Rad_Atm_Profile(Sfc_Idx) + & (Rad_Atm_Profile(Sfc_Idx+1) - Rad_Atm_Profile(Sfc_Idx)) * Profile_Weight Trans_Atm = Trans_Atm_Profile(Sfc_Idx) + & (Trans_Atm_Profile(Sfc_Idx+1) - Trans_Atm_Profile(Sfc_Idx)) * Profile_Weight Rad_Atm_Sfc = Rad_Atm + Trans_Atm * Sfc_Rad Bt_Atm_Sfc = PLANCK_TEMP_FAST(Chan_Idx,Rad_Atm_Sfc) end subroutine COMPUTE_CHANNEL_ATM_SFC_RAD_BT !------------------------------------------------------------------------------ ! Routine to compute some needed radiative transfer terms for the IR channels ! ! Input: Sfc_Idx - level of the surface in the profiles ! Profile_Weight - interpolation weight for estimated the surface in the ! profiles ! Rad_Atm_Dwn_Profile - profile of radiance emitted from level to space ! ! Output: Rad_Atm_Dwn_Sfc - total radiance due to atmospheric emission at surface !------------------------------------------------------------------------------ subroutine COMPUTE_CHANNEL_ATM_DWN_SFC_RAD( & Sfc_Idx, & Profile_Weight, & Rad_Atm_Dwn_Profile, & Rad_Atm_Dwn_Sfc) integer, intent(in):: Sfc_Idx real, intent(in):: Profile_Weight real, intent(in), dimension(:):: Rad_Atm_Dwn_Profile real, intent(out):: Rad_Atm_Dwn_Sfc Rad_Atm_Dwn_Sfc = Rad_Atm_Dwn_Profile(Sfc_Idx) + & (Rad_Atm_Dwn_Profile(Sfc_Idx+1) - Rad_Atm_Dwn_Profile(Sfc_Idx)) * Profile_Weight end subroutine COMPUTE_CHANNEL_ATM_DWN_SFC_RAD !---------------------------------------------------------------------------------------- ! This routine computes radiative transfer terms such as Rad_Atm ! Trans_Atm, and clear-sky radinace and brightness temperature ! ! Input: Sfc_Level_Idx - level just above the surface in the profiles ! Prof_Weight - interpolation weight for interpolating to surface level ! Lon_Idx = longitude index of NWP cell ! Lat_Idx = latitude index of NWP cell ! Zen_Idx = zenith angle index of RTM profile ! ! Output: (note this passed through global arrays) ! Rad_Atm_ChX_Rtm = Radiance Emitted by Atmosphere in Channel X ! Trans_Atm_ChX_Rtm = Transmission by Atmosphere in Channel X ! Rad_Clear_ChX_Rtm = Radiance at TOA for clear skies (atm + sfc) ! Bt_Clear_ChX_Rtm = BT at TOA for clear skies (atm + sfc) ! ! NOTE: These pixels are not smoothed and are blocky at NWP scale. Should fix !---------------------------------------------------------------------------------------- subroutine COMPUTE_CHANNEL_RT(Sfc_Level_Idx,Prof_Weight,Lon_Idx,Lat_Idx,Elem_Idx,Line_Idx,Zen_Idx) integer, intent(in):: Sfc_Level_Idx real, intent(in):: Prof_Weight integer, intent(in):: Lon_Idx integer, intent(in):: Lat_Idx integer, intent(in):: Elem_Idx integer, intent(in):: Line_Idx integer, intent(in):: Zen_Idx integer:: Chan_Idx real:: Sfc_Ref real:: Rad_Ch20_Temp real:: Rad_Clear_Ch20_Solar_Rtm real:: Bt_Clear_Ch20_Solar_Rtm !-------------------------------------------------------------- ! Solar-Only channels, 1,2,6,7,DNB(44) !-------------------------------------------------------------- do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (Ch(Chan_Idx)%Obs_Type /= SOLAR_OBS_TYPE .and. & Ch(Chan_Idx)%Obs_Type /= LUNAR_OBS_TYPE) cycle select case (Chan_Idx) case (1,2,5,6,7,44) if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%YES) then if (allocated( Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Total_Profile )) then Ch(Chan_Idx)%Trans_Atm_Total(Elem_Idx,Line_Idx) = & Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%Ch(Chan_Idx)%Trans_Atm_Total_Profile(Sfc_Level_Idx) + & (Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%Ch(Chan_Idx)%Trans_Atm_Total_Profile(Sfc_Level_Idx+1) - & Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%Ch(Chan_Idx)%Trans_Atm_Total_Profile(Sfc_Level_Idx)) * Prof_Weight end if end if end select end do !-------------------------------------------------------------- ! IR-only channels, 20-38 (except 26), 42, 43 !-------------------------------------------------------------- !--- upwelling do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%NO) cycle if (Ch(Chan_Idx)%Obs_Type /= THERMAL_OBS_TYPE .and. & Ch(Chan_Idx)%Obs_Type /= MIXED_OBS_TYPE ) cycle call COMPUTE_CHANNEL_ATM_SFC_RAD_BT( & Chan_Idx, & Sfc_Level_Idx, & Prof_Weight, & Ch(Chan_Idx)%Sfc_Emiss(Elem_Idx,Line_Idx), & NWP_PIX%Tsfc(Elem_Idx,Line_Idx), & Rtm(Lon_Idx,Lat_Idx)%D(Zen_Idx)%Ch(Chan_Idx)%Rad_Atm_Profile, & Rtm(Lon_Idx,Lat_Idx)%D(Zen_Idx)%Ch(Chan_Idx)%Trans_Atm_Profile, & Ch(Chan_Idx)%Rad_Atm(Elem_Idx,Line_Idx), & Ch(Chan_Idx)%Trans_Atm(Elem_Idx,Line_Idx), & Ch(Chan_Idx)%Rad_Toa_Clear(Elem_Idx,Line_Idx), & Ch(Chan_Idx)%Bt_Toa_Clear(Elem_Idx,Line_Idx)) end do !--- downwelling (only channel 31/38) do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (Chan_Idx == 31 .OR. Chan_Idx == 38) then if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%NO) cycle call COMPUTE_CHANNEL_ATM_DWN_SFC_RAD( & Sfc_Level_Idx, & Prof_Weight, & Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%Ch(Chan_Idx)%Rad_Atm_Dwn_Profile, & Ch(Chan_Idx)%Rad_Atm_Dwn_Sfc(Elem_Idx,Line_Idx)) else cycle endif end do !-------------------------------------------------------------- !-- Add Solar to Ch20 clear variables !-------------------------------------------------------------- if (Sensor%Chan_On_Flag_Default(20) == sym%YES) then !--- add in solar component - does not account for glint Trans_Atm_Ch20_Solar_Rtm(Elem_Idx,Line_Idx) = & Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(20)%Trans_Atm_Solar_Profile(Sfc_Level_Idx) + & (Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(20)%Trans_Atm_Solar_Profile(Sfc_Level_Idx+1) - & Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(20)%Trans_Atm_Solar_Profile(Sfc_Level_Idx)) * & Prof_Weight Trans_Atm_Ch20_Solar_Total_Rtm(Elem_Idx,Line_Idx) = & Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(20)%Trans_Atm_Total_Profile(Sfc_Level_Idx) + & (Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(20)%Trans_Atm_Total_Profile(Sfc_Level_Idx+1) - & Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(20)%Trans_Atm_Total_Profile(Sfc_Level_Idx)) * & Prof_Weight Rad_Clear_Ch20_Solar_Rtm = ch(20)%Rad_Toa_Clear(Elem_Idx,Line_Idx) Bt_Clear_Ch20_Solar_Rtm = ch(20)%Bt_Toa_Clear(Elem_Idx,Line_Idx) if (Geo%Cossolzen(Elem_Idx,Line_Idx) >= 0.0) then Sfc_Ref = 1.0 - ch(20)%Sfc_Emiss(Elem_Idx,Line_Idx) Rad_Ch20_Temp = ch(20)%Rad_Toa_Clear(Elem_Idx,Line_Idx) + & Trans_Atm_Ch20_Solar_Total_Rtm(Elem_Idx,Line_Idx) * & Sfc_Ref * (Geo%Cossolzen(Elem_Idx,Line_Idx)*Solar_Ch20_Nu / pi) Rad_Clear_Ch20_Solar_Rtm = Rad_Ch20_Temp Bt_Clear_Ch20_Solar_Rtm = PLANCK_TEMP_FAST(20,Rad_Ch20_Temp) end if ch(20)%Rad_Toa_Clear(Elem_Idx,Line_Idx) = Rad_Clear_Ch20_Solar_Rtm ch(20)%Bt_Toa_Clear(Elem_Idx,Line_Idx) = Bt_Clear_Ch20_Solar_Rtm end if end subroutine COMPUTE_CHANNEL_RT !------------------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------------------- subroutine COMPUTE_CH20_EMISSIVITY(Elem_Idx,Line_Idx) integer, intent(in):: Elem_Idx integer, intent(in):: Line_Idx real:: Rad_Ch20_Temp real:: Ch20_Sfc_Rad if ((Sensor%Chan_On_Flag_Default(20) == sym%YES) .and. (Sensor%Chan_On_Flag_Default(31)==sym%YES)) then Ch20_Sfc_Rad = ch(20)%Sfc_Emiss(Elem_Idx,Line_Idx) * PLANCK_RAD_FAST(20,NWP_PIX%Tsfc(Elem_Idx,Line_Idx)) Rad_Ch20_Temp = PLANCK_RAD_FAST(20,ch(31)%Bt_Toa_Clear(Elem_Idx,Line_Idx)) Ch(20)%Emiss_Rel_11um_Clear(Elem_Idx,Line_Idx) = ch(20)%Rad_Toa_Clear(Elem_Idx,Line_Idx) / Rad_Ch20_Temp end if if ((Sensor%Chan_On_Flag_Default(20) == sym%YES) .and. (Sensor%Chan_On_Flag_Default(38)==sym%YES)) then Ch20_Sfc_Rad = ch(20)%Sfc_Emiss(Elem_Idx,Line_Idx) * PLANCK_RAD_FAST(20,NWP_PIX%Tsfc(Elem_Idx,Line_Idx)) Rad_Ch20_Temp = PLANCK_RAD_FAST(20,ch(38)%Bt_Toa_Clear(Elem_Idx,Line_Idx)) Ch(20)%Emiss_Rel_10_4um_Clear(Elem_Idx,Line_Idx) = ch(20)%Rad_Toa_Clear(Elem_Idx,Line_Idx) / Rad_Ch20_Temp end if end subroutine COMPUTE_CH20_EMISSIVITY !------------------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------------------- subroutine COMPUTE_TROPOPAUSE_EMISSIVITIES(Elem_Idx,Line_Idx,Lon_Idx,Lat_Idx,Zen_Idx) integer, intent(in):: Elem_Idx integer, intent(in):: Line_Idx integer, intent(in):: Lon_Idx integer, intent(in):: Lat_Idx integer, intent(in):: Zen_Idx integer:: Chan_Idx integer:: Lev_Bnd integer:: dim1 integer:: dim2 dim1 = Image%Number_Of_Elements dim2 = Image%Number_Of_Lines_Per_Segment Lev_Bnd = Rtm(Lon_Idx,Lat_Idx)%Tropo_Level !--- check for missing tropopause level if (Rtm(Lon_Idx,Lat_Idx)%Tropo_Level == 0) then Bad_Pixel_Mask(Elem_Idx,Line_Idx) = sym%YES return end if do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max select case (Chan_Idx) case(20,27,29,31,32,33,37,38) if (Sensor%Chan_On_Flag_Default(Chan_Idx)==sym%YES) then ch(Chan_Idx)%Emiss_Tropo(Elem_Idx,Line_Idx) = & EMISSIVITY(ch(Chan_Idx)%Rad_Toa(Elem_Idx,Line_Idx), & ch(Chan_Idx)%Rad_Toa_Clear(Elem_Idx,Line_Idx), & Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_BB_Cloud_Profile(Lev_Bnd)) !if (isnan(ch(Chan_Idx)%Emiss_Tropo(Elem_Idx,Line_Idx))) then ! print *, "Nan in emiss_tropo ", chan_idx, ch(Chan_Idx)%Rad_Toa_Clear(Elem_Idx,Line_Idx) !endif end if end select end do end subroutine COMPUTE_TROPOPAUSE_EMISSIVITIES !------------------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------------------- subroutine COMPUTE_BETA_RATIOES(Elem_Idx,Line_Idx) integer, intent(in):: Elem_Idx integer, intent(in):: Line_Idx !--- compute 11 and 12 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(31) == sym%YES .and. & Sensor%Chan_On_Flag_Default(32) == sym%YES) then Beta_11um_12um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(32)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(31)%Emiss_Tropo(Elem_Idx,Line_Idx)) end if !--- compute 10 and 12 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(38) == sym%YES .and. & Sensor%Chan_On_Flag_Default(32) == sym%YES) then Beta_104um_12um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(32)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(38)%Emiss_Tropo(Elem_Idx,Line_Idx)) end if !--- compute 11 and 8.5 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(31) == sym%YES .and. & Sensor%Chan_On_Flag_Default(29) == sym%YES) then Beta_11um_85um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(29)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(31)%Emiss_Tropo(Elem_Idx,Line_Idx)) endif !--- compute 10 and 8.5 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(38) == sym%YES .and. & Sensor%Chan_On_Flag_Default(29) == sym%YES) then Beta_104um_85um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(29)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(38)%Emiss_Tropo(Elem_Idx,Line_Idx)) endif !--- compute 11 and 6.7 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(31) == sym%YES .and. & Sensor%Chan_On_Flag_Default(27) == sym%YES) then Beta_11um_67um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(27)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(31)%Emiss_Tropo(Elem_Idx,Line_Idx)) end if !--- compute 10 and 6.7 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(38) == sym%YES .and. & Sensor%Chan_On_Flag_Default(27) == sym%YES) then Beta_104um_67um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(27)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(38)%Emiss_Tropo(Elem_Idx,Line_Idx)) end if !--- compute 11 and 13.3 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(31) == sym%YES .and. & Sensor%Chan_On_Flag_Default(33) == sym%YES) then Beta_11um_133um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(33)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(31)%Emiss_Tropo(Elem_Idx,Line_Idx)) endif !--- compute 10 and 13.3 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(38) == sym%YES .and. & Sensor%Chan_On_Flag_Default(33) == sym%YES) then Beta_104um_133um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(33)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(38)%Emiss_Tropo(Elem_Idx,Line_Idx)) endif !--- compute 11 and 10 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(31) == sym%YES .and. & Sensor%Chan_On_Flag_Default(38) == sym%YES) then Beta_11um_104um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(38)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(31)%Emiss_Tropo(Elem_Idx,Line_Idx)) endif !--- compute 10 and 11 beta ratio at tropopause if (Sensor%Chan_On_Flag_Default(31) == sym%YES .and. & Sensor%Chan_On_Flag_Default(38) == sym%YES) then Beta_104um_11um_Tropo_Rtm(Elem_Idx,Line_Idx) = BETA_RATIO( & ch(31)%Emiss_Tropo(Elem_Idx,Line_Idx), & ch(38)%Emiss_Tropo(Elem_Idx,Line_Idx)) endif end subroutine COMPUTE_BETA_RATIOES !==================================================================== ! FUNCTION Name: BETA_RATIO ! ! Function: ! Computes the beta ratio for two Emissivities. ! ! Input: Emiss_top - emissivity in the numerator ! Emiss_bot - emissivity in the denominator ! ! Output: Beta - the beta value from the two emissivities ! !==================================================================== function BETA_RATIO(Emiss_top, Emiss_bot) result(beta) real(kind=real4), intent(in) :: Emiss_top real(kind=real4), intent(in) :: Emiss_bot real(kind=real4) :: beta beta = Missing_Value_Real4 if (Emiss_top > 0.0 .and. Emiss_top < 1.0 .and. & Emiss_bot > 0.0 .and. Emiss_bot < 1.0) then beta = alog(1.0 - Emiss_top)/alog(1.0 - Emiss_bot) end if return end function BETA_RATIO !==================================================================== ! Function Name: EMISSIVITY ! ! Function: ! Computes the effective emissivity ! ! Input: Rad_Toa - channel radiance at top of atmosphere(toa) ! Rad_Clear_Tau - channel radiance at toa for clear skies ! Rad_Cloud_BB_Toa - channel radiance at TOA if cloud were a Black-Body ! ! Output: Emiss - the effective cloud emissivity ! !==================================================================== function EMISSIVITY(Radiance_Toa, Radiance_Clear_Toa, Radiance_Cloud_BB_Toa) result(Emiss) real(kind=real4), intent(in) :: Radiance_Toa real(kind=real4), intent(in) :: Radiance_Clear_Toa real(kind=real4), intent(in) :: Radiance_Cloud_BB_Toa real(kind=real4) :: Emiss Emiss = Missing_Value_Real4 if (Radiance_Cloud_BB_Toa /= Radiance_Clear_Toa) then Emiss = (Radiance_Toa - Radiance_Clear_Toa) / & (Radiance_Cloud_BB_Toa - Radiance_Clear_Toa) end if return end function EMISSIVITY !==================================================================== ! Function Name: NADIR_EMISSIVITY ! ! Function: Convert a slant-path emissivity to a nadir value ! ! Input: Emiss - slant path emissivity ! Cos_Zen - cosine of viewing zenith angle ! ! Output: Nadir_Emiss - emissivity if viewed nadir ! !==================================================================== function NADIR_EMISSIVITY(Emiss,Cos_Zen) result(Nadir_Emiss) real(kind=real4), intent(in) :: Emiss real(kind=real4), intent(in) :: Cos_Zen real(kind=real4) :: Nadir_Emiss real:: Tau_Temp Nadir_Emiss = Missing_Value_Real4 if ((Emiss > 0.0) .and. (Emiss < 1.00)) then Tau_Temp = -Cos_Zen*alog(1.0 - Emiss) Nadir_Emiss = 1.0 - exp(-Tau_Temp) end if return end function NADIR_EMISSIVITY !=============================================================================== ! !=============================================================================== subroutine COPY_LOCAL_RTM_TO_GLOBAL_RTM_STRUCTURE(Lon_Idx,Lat_Idx,Zen_Idx) integer, intent(in):: Lon_Idx integer, intent(in):: Lat_Idx integer, intent(in):: Zen_Idx integer:: Chan_Idx do Chan_Idx = Chan_Idx_Min, Chan_Idx_Max if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%NO) cycle Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Total_Profile = Trans_Atm_Total_Prof(:,Chan_Idx) Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Solar_Profile = Trans_Atm_Solar_Prof(:,Chan_Idx) Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Profile = Trans_Atm_Prof(:,Chan_Idx) Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Profile = Rad_Atm_Prof(:,Chan_Idx) Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_BB_Cloud_Profile = Rad_BB_Cloud_Prof(:,Chan_Idx) end do Chan_Idx = 31 if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%YES) then Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Dwn_Profile = Rad_Atm_Dwn_Prof(:,Chan_Idx) end if Chan_Idx = 38 if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%YES) then Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Dwn_Profile = Rad_Atm_Dwn_Prof(:,Chan_Idx) end if end subroutine COPY_LOCAL_RTM_TO_GLOBAL_RTM_STRUCTURE !============================================================================== ! !============================================================================== subroutine ALLOCATE_GLOBAL_RTM_STRUCTURE_ELEMENT(Lon_Idx,Lat_Idx,Zen_Idx) integer, intent(in):: Lon_Idx integer, intent(in):: Lat_Idx integer, intent(in):: Zen_Idx integer:: Alloc_Status integer:: Chan_Idx do Chan_Idx = Chan_Idx_Min,Chan_Idx_Max if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%YES) then allocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Profile(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Solar_Profile(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Trans_Atm_Total_Profile(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Profile(NLevels_Rtm),stat=Alloc_Status) allocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_BB_Cloud_Profile(NLevels_Rtm),stat=Alloc_Status) endif enddo Chan_Idx = 31 if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%YES) then allocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Dwn_Profile(NLevels_Rtm),stat=Alloc_Status) endif Chan_Idx = 38 if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%YES) then allocate(Rtm(Lon_Idx,Lat_Idx)%d(Zen_Idx)%ch(Chan_Idx)%Rad_Atm_Dwn_Profile(NLevels_Rtm),stat=Alloc_Status) endif end subroutine ALLOCATE_GLOBAL_RTM_STRUCTURE_ELEMENT !=============================================================================== ! compute an inversion profile ! a level is considered in a inversion if it is colder than the level above it !=============================================================================== subroutine INVERSION_PROFILE(T_Prof,Inver_Prof) real, dimension(:), intent(in):: T_Prof integer (kind=int1), dimension(:), intent(out):: Inver_Prof integer:: N_Levels integer:: Lev_Idx Inver_Prof = 0 n_levels = size (T_prof) do Lev_Idx = 2, N_Levels if (T_Prof(Lev_idx) < T_Prof(Lev_Idx-1)) Inver_Prof(Lev_Idx) = 1 enddo end subroutine INVERSION_PROFILE !------------------------------------------------------------ ! compute the single scater and aerosol reflectance ! this assumes that the gas is mixed in with scattering ! ! Input: ! scatangle = scattering angle ! coszen = cosine of viewing zenith angle ! cossolzen = cosine of solar zenith angle ! tau_ray = rayleigh optical depth ! tau_gas = gaseous optical depth ! tau_aer = aerosol optical depth ! wo_aer = aerosol single scatter albedo ! g_aer = aerosol asymmetry parameter ! cloud_albedo_zen = cloud or surface albedo from illumination along coszen ! cloud_albedo_solzen = cloud or surface albedo from illumination along ! cosolzen ! note- albedoes are (0-1) ! ! Output: ! Ref_SS = single scatter reflectance from the layer above the ! surface or cloud at the top of atmosphere (0-110%) ! ! Reference: ! Correction of Rayleigh scattering effects in cloud optical thickness retrievals ! Menghua Wang and Michael D. King ! JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 102, NO. D22, PAGES 25,915-25,926, ! NOVEMBER 27, 1997 !----------------------------------------------------------- subroutine COMPUTE_CLEAR_SKY_SCATTER(Tau_Aer, & wo_Aer, & g_Aer, & Tau_Ray, & Tau_gas, & scatangle, & Coszen, & Cossolzen, & Cloud_Albedo_View, & Cloud_Albedo_Sun, & Ref_ss) real, intent(in):: Tau_Aer real, intent(in):: wo_Aer real, intent(in):: g_Aer real, intent(in):: Tau_Ray real, intent(in):: Tau_Gas real, intent(in):: Scatangle real, intent(in):: Coszen real, intent(in):: Cossolzen real, intent(in):: Cloud_Albedo_View real, intent(in):: Cloud_Albedo_Sun real, intent(out):: Ref_ss real:: Airmass real:: P_Aer real:: P_Ray real:: Tau_Total real:: Tau_Scat_Total real:: Trans_Total real:: Tau_Iso_Total real:: Trans_Iso_Total_View real:: Trans_Iso_Total_Sun real:: Tau_Iso_Scat_Total real:: mu real:: Pf real:: wo real:: Ref_ss_a real:: Ref_ss_b real:: Ref_ss_c !--- compute cosine of scattering angle mu = cos(scatangle*dtor) !-- compute Rayleigh phase function Airmass = 1.0/Coszen + 1.0/Cossolzen P_Ray = 0.75*(1.0 + mu**2) !--- compute total transmission Tau_Total = Tau_Aer + Tau_Ray + Tau_gas Trans_Total = exp(-Tau_Total*Airmass) Tau_Iso_Total = (1.0-g_Aer)*Tau_Aer + Tau_Ray + Tau_gas Trans_Iso_Total_View = exp(-Tau_Iso_Total/Coszen) Trans_Iso_Total_Sun = exp(-Tau_Iso_Total/Cossolzen) !--- compute total scattering optical depth Tau_Scat_Total = wo_Aer*Tau_Aer + Tau_Ray Tau_Iso_Scat_Total = wo_Aer*(1.0-g_Aer)*Tau_Aer + Tau_Ray !--- single scatter albedo wo = (wo_Aer*Tau_Aer + Tau_Ray)/ ( Tau_Total ) !aerosol phase function (Henyey-Greenstein) P_Aer = (1.0 - g_Aer**2)/( (1.0 + g_Aer**2 - 2.0*g_Aer*mu)**(1.5) ) !--- compute effective phase function Pf = P_aer if (Tau_Scat_Total > 0.0) then Pf = (wo_Aer*Tau_Aer*P_Aer + Tau_Ray*P_Ray)/(Tau_Scat_Total) endif !--- compute single scatter reflectance (0-100%) !--- single scattering in the layer without hitting surface Ref_ss_a = wo*Pf/(4.0*Airmass*Coszen*Cossolzen) * (1.0 - Trans_Total ) !--- single scattering in the layer without hitting surface Ref_ss_b = (Tau_Iso_Scat_Total / (2.0*Cossolzen)) * & Trans_Iso_Total_View * Cloud_Albedo_View Ref_ss_c = (Tau_Iso_Scat_Total / (2.0*Coszen)) * & Trans_Iso_Total_Sun * Cloud_Albedo_Sun Ref_ss = 100.0*(Ref_ss_a + Ref_ss_b + Ref_ss_c) end subroutine COMPUTE_CLEAR_SKY_SCATTER !------------------------------------------------------------ ! atmospheric correction for reflectance channels ! ! note, original optical depths taken from ! Ignatov and Stowe, 2002, vol 59, JAS, pg 313-334 ! ! AKH nows runs modtran to generate these numbers !----------------------------------------------------------- subroutine ATMOS_CORR(Line_Idx_Min,Num_Lines) integer, intent(in):: Line_Idx_Min integer, intent(in):: Num_Lines integer:: Line_Idx integer:: Elem_Idx integer:: Line_Idx_Max integer:: Elem_Idx_Max integer:: Elem_Idx_Min integer:: Num_Elements integer:: Chan_Idx real:: Ref_ss real:: Albedo_View real:: Albedo_Sun real:: Trans_Total real:: Tau_Total real:: Tau_Ray real:: Tau_Gas real:: Tau_H2O real:: Tau_Aer real:: Wo_Aer real:: G_Aer real:: Tpw_Ac real:: Zc, Pc real:: Source_Zen real:: Cos_Source_Zen real:: Airmass_Factor real:: Scattering_Angle !--- compute gas terms real, parameter :: h_h2o = 2000.0 !m real, parameter :: P_Sfc = 1013.25 Elem_Idx_Min = 1 Num_Elements = Image%Number_Of_Elements Elem_Idx_Max = Elem_Idx_Min + Num_Elements - 1 Line_Idx_Max = Line_Idx_Min + Num_Lines - 1 !---------- loop over scan lines line_loop: do Line_Idx = Line_Idx_Min, Line_Idx_Max element_loop: do Elem_Idx = Elem_Idx_Min, Elem_Idx_Max !--- check for bad individual pixels if (Bad_Pixel_Mask(Elem_Idx,Line_Idx) == sym%YES) cycle channel_loop: do Chan_Idx = 1,Nchan_Clavrx if (Sensor%Chan_On_Flag_Default(Chan_Idx) == sym%NO) cycle if (Chan_Idx /= 1 .and. Chan_Idx /= 2 .and. Chan_Idx /= 3 .and. & Chan_Idx /= 4 .and. Chan_Idx /= 5 .and. Chan_Idx /= 6 .and. & Chan_Idx /= 7 .and. Chan_Idx /= 26 .and. Chan_Idx /= 44) cycle !--- check for valid data if (Chan_Idx /= 44 .and. ch(Chan_Idx)%Ref_Toa(Elem_Idx,Line_Idx) == Missing_Value_Real4) cycle if (Chan_Idx == 44) then if (ch(Chan_Idx)%Ref_Lunar_Toa(Elem_Idx,Line_Idx) == Missing_Value_Real4) cycle endif !--- set source angle Source_Zen = Geo%SolZen(Elem_Idx,Line_Idx) Scattering_Angle = Geo%Scatangle(Elem_Idx,Line_Idx) if (Chan_Idx == 44) then Source_Zen = Geo%LunZen(Elem_Idx,Line_Idx) Scattering_Angle = Geo%Scatangle_Lunar(Elem_Idx,Line_Idx) endif !--- check for appropriate illumination if (Source_Zen >= 90.0) cycle !--- compute gas terms Tpw_Ac = NWP_PIX%Tpw(Elem_Idx,Line_Idx) if (Zc_Opaque_Cloud(Elem_Idx,Line_Idx) /= MISSING_VALUE_REAL4) then Zc = max(0.0,Zc_Opaque_Cloud(Elem_Idx,Line_Idx)) Tpw_Ac = Tpw_Ac * exp(-Zc/h_h2o) endif Tau_H2O = Solar_Rtm%Tau_H2O_Coef(Chan_Idx,1) + Solar_Rtm%Tau_H2O_Coef(Chan_Idx,2)*Tpw_Ac + & Solar_Rtm%Tau_H2O_Coef(Chan_Idx,3)*(Tpw_Ac**2) Tau_Gas = max(0.0,Tau_H2O) + Solar_Rtm%Tau_O3(Chan_Idx) + Solar_Rtm%Tau_O2(Chan_Idx) & + Solar_Rtm%Tau_CO2(Chan_Idx) + Solar_Rtm%Tau_CH4(Chan_Idx) Tau_Ray = Solar_Rtm%Tau_Ray(Chan_Idx) Tau_Aer = Solar_Rtm%Tau_Aer(Chan_Idx) if (Pc_Opaque_Cloud(Elem_Idx,Line_Idx) /= MISSING_VALUE_REAL4) then Pc = min(1000.0,max(0.0,Pc_Opaque_Cloud(Elem_Idx,Line_Idx))) Tau_Ray = Tau_Ray * Pc_Opaque_Cloud(Elem_Idx,Line_Idx) / P_Sfc Tau_Aer = Tau_Aer * (Pc_Opaque_Cloud(Elem_Idx,Line_Idx) / P_Sfc)**2 endif Wo_Aer = Solar_Rtm%Wo_Aer(Chan_Idx) G_Aer = Solar_Rtm%G_Aer(Chan_Idx) !------------------------------------------------------------------------ ! select gas and surface reflectance parameters !------------------------------------------------------------------------ select case (Chan_Idx) case(1) if (ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) / 100.0 else Albedo_View = Ch1_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = Ch1_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif case(2) if (ch(2)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = ch(2)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) / 100.0 else Albedo_View = Ch2_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = Ch2_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif case(3) if (ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) / 100.0 else Albedo_View = Ch1_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = Ch1_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif case(4) if (ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) / 100.0 else Albedo_View = Ch1_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = Ch1_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif case(5) if (ch(5)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = ch(5)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) / 100.0 else Albedo_View = Ch6_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 !Note there is no Ch5_Sfc_Alb_Umd endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = Ch5_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif case(6) if (ch(6)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = ch(6)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) / 100.0 else Albedo_View = Ch6_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = Ch6_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif case(7) if (ch(7)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = ch(7)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) / 100.0 else Albedo_View = Ch6_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 !Note there is no Ch7_Sfc_Alb_Umd endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = Ch7_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif case(26) if (ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) / 100.0 else Albedo_View = Ch1_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 !Note there is no Ch26_Sfc_Alb_Umd endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = Ch1_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) / 100.0 endif case(44) !DNB - use mean of ch1 and ch2 for sfc reflectance if (ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4 & .and. ch(2)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx) /= Missing_Value_Real4) then Albedo_View = 0.5*(ch(1)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx)+ch(2)%Sfc_Ref_White_Sky(Elem_Idx,Line_Idx)) / 100.0 else Albedo_View = 0.5*(Ch1_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) + Ch2_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx))) / 100.0 endif if (Sfc%Snow(Elem_Idx,Line_Idx) /= sym%NO_SNOW) then Albedo_View = 0.5*(Ch1_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx)) & + Ch2_Snow_Sfc_Alb_Umd(Sfc%Sfc_Type(Elem_Idx,Line_Idx))) / 100.0 endif end select !-- assume lambertian Albedo_Sun = Albedo_View Cos_Source_Zen = Missing_Value_Real4 if (Source_Zen >= 0.0 .and. Source_Zen < 90.0) Cos_Source_Zen = cos(Source_Zen*DTOR) Airmass_Factor = 1.0 / Geo%CosZen(Elem_Idx,Line_Idx) + & 1.0 / Cos_Source_Zen !--- compute atmospheric scattering call COMPUTE_CLEAR_SKY_SCATTER(Tau_Aer, & Wo_Aer, & G_Aer, & Tau_Ray, & Tau_Gas, & Scattering_Angle, & Geo%Coszen(Elem_Idx,Line_Idx), & Cos_Source_Zen, & Albedo_View, & Albedo_Sun, & Ref_ss) !--- compute total transmission for combining terms Tau_Total = Tau_Aer + Tau_Ray + Tau_Gas Trans_Total = exp(-Tau_Total*Airmass_Factor) if (Chan_Idx /= 44) then !--- compute atmospherically corrected reflectance (at sfc level)s ch(Chan_Idx)%Ref_Sfc(Elem_Idx,Line_Idx) = (ch(Chan_Idx)%Ref_Toa(Elem_Idx,Line_Idx) - Ref_ss) / Trans_Total !--- compute top of clear-sky atmosphere reflectance ch(Chan_Idx)%Ref_Toa_Clear(Elem_Idx,Line_Idx) = Ref_ss + Trans_Total*100.0*Albedo_View else !--- compute atmospherically corrected reflectance (at sfc level)s ch(Chan_Idx)%Ref_Lunar_Sfc(Elem_Idx,Line_Idx) = (ch(Chan_Idx)%Ref_Lunar_Toa(Elem_Idx,Line_Idx) - Ref_ss) / Trans_Total !--- compute top of clear-sky atmosphere reflectance ch(Chan_Idx)%Ref_Lunar_Toa_Clear(Elem_Idx,Line_Idx) = Ref_ss + Trans_Total*100.0*Albedo_View endif end do channel_loop end do element_loop end do line_loop end subroutine ATMOS_CORR ! !--- end of module ! end module RT_UTILITIES_MOD