diff --git a/src/gsi/cads.f90 b/src/gsi/cads.f90 new file mode 100644 index 0000000000..6d7e584ef1 --- /dev/null +++ b/src/gsi/cads.f90 @@ -0,0 +1,2230 @@ +module cads +!$$$ module documentation block +! +! module: cads +! prgmmr: Jung +! +! abstract: module containing subroutines for the cloud and aerosol detection software +! +! program history log: +! +! +! +! subroutines included: +! +! +! remarks: variable definitions +! +! +!$$$ end documentation block + + + use kinds, only: i_kind, r_kind + implicit none + save + +! set default to private + private +! set routines to public + public :: cloud_aerosol_detection + public :: cads_setup_cloud + public :: Cloud_Detect_Type + public :: cads_imager_calc + + public :: M__Sensor,N__Num_Bands,N__GradChkInterval,N__Band_Size,N__Bands,N__Window_Width, & + N__Window_Bounds,R__BT_Threshold,R__Grad_Threshold,R__Window_Grad_Threshold, L__Do_Quick_Exit, & + L__Do_CrossBand, N__BandToUse,L__Do_Imager_Cloud_Detection, N__Num_Imager_Chans, & + N__Num_Imager_Clusters,N__Imager_Chans,R__Stddev_Threshold,R__Coverage_Threshold, & + R__FG_Departure_Threshold + + INTEGER(i_kind) :: M__Sensor ! Unique ID for sensor + INTEGER(i_kind) :: N__Num_Bands ! Number of channel bands + INTEGER(i_kind), POINTER :: N__GradChkInterval(:) ! Window width used in gradient calculation + INTEGER(i_kind), POINTER :: N__Band_Size(:) ! Number of channels in each band + INTEGER(i_kind), POINTER :: N__Bands(:,:) ! Channel lists + INTEGER(i_kind), POINTER :: N__Window_Width(:) ! Smoothing filter window widths per band + INTEGER(i_kind), POINTER :: N__Window_Bounds(:,:) ! Channels in the spectral window gradient check + INTEGER(i_kind), POINTER :: N__BandToUse(:) ! Band number assignment for each channel + LOGICAL :: L__Do_Quick_Exit ! On/off switch for the Quick Exit scenario + LOGICAL :: L__Do_CrossBand ! On/off switch for the cross-band method + REAL(r_kind), POINTER :: R__BT_Threshold(:) ! BT threshold for cloud contamination + REAL(r_kind), POINTER :: R__Grad_Threshold(:) ! Gradient threshold for cloud contamination + REAL(r_kind), POINTER :: R__Window_Grad_Threshold(:) ! Threshold for window gradient check in QE + + LOGICAL :: L__Do_Imager_Cloud_Detection ! On/off switch for the imager cloud detection + INTEGER(i_kind) :: N__Num_Imager_Chans ! No. of imager channels + INTEGER(i_kind) :: N__Num_Imager_Clusters ! No. of clusters to be expected + INTEGER(i_kind),POINTER :: N__Imager_Chans(:) ! List of imager channels + REAL(r_kind),POINTER :: R__Stddev_Threshold(:) ! St. Dev. threshold, one for each imager channel + REAL(r_kind) :: R__Coverage_Threshold ! Threshold for fractional coverage of a cluster + REAL(r_kind) :: R__FG_Departure_Threshold ! Threshold for imager FG departure + + +! set passed variables to public + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and MeteoFrance. +! +! Copyright 2020, EUMETSAT, All Rights Reserved. + +! * CADS_Module * +! A. Collard ECMWF 01/02/06 + +! * PURPOSE * +! ----------- +! Sets up structures to be used in processing of advanced IR sounders. + +! * MODIFICATIONS * +! ----------------- +! 01/02/06 A.Collard 1.0 Original export version. +! 17/11/09 R.Eresmaa 1.1 Include parameters of the Quick Exit / +! long-wave window gradient check. +! 11/11/11 R.Eresmaa 1.2 Add processing capability for CrIS. +! 03/12/13 R.Eresmaa 2.0 Add imager-assisted cloud detection. +! 10/11/15 R.Eresmaa 2.2 Changed instrument ID naming convention. +! Changed aerosol detection parameters. +! 20/12/16 R.Eresmaa 2.3 Remove aerosol detection parameters. +! 05/02/19 R.Eresmaa 2.4 Explicit KIND specifications. +! 16/04/20 R.Eresmaa 3.0 Combine cloud and aerosol detection, rename. +! Include aerosol type recognition. +! Include land sensitivity parameters. +! Include trace gas detection. Rename. + + + INTEGER(i_kind), PARAMETER :: INST_ID_AIRS = 11 + INTEGER(i_kind), PARAMETER :: INST_ID_IASI = 16 + INTEGER(i_kind), PARAMETER :: INST_ID_CRIS = 27 + INTEGER(i_kind), PARAMETER :: INST_ID_IRS = 57 + INTEGER(i_kind), PARAMETER :: INST_ID_IASING = 59 + INTEGER(i_kind), PARAMETER :: INST_ID_IKFS2 = 94 + INTEGER(i_kind), PARAMETER :: INST_ID_HIRAS = 97 + INTEGER(i_kind), PARAMETER :: INST_ID_GIIRS = 98 + + INTEGER(i_kind), PARAMETER :: JP__MIN_SENSOR_INDEX = INST_ID_AIRS + INTEGER(i_kind), PARAMETER :: JP__MAX_SENSOR_INDEX = INST_ID_GIIRS + + TYPE Aerosol_Detect_Type + INTEGER(i_kind) :: M__Sensor ! Unique ID for sensor + INTEGER(i_kind) :: N__Num_Aerosol_Tests ! Number of aerosol detection tests + INTEGER(i_kind), POINTER :: N__Num_Regression(:) ! Number of conversion coefficients for AOD + INTEGER(i_kind), POINTER :: N__Num_Aerosol_Chans(:) ! Number of aerosol detection channels + INTEGER(i_kind), POINTER :: N__Aerosol_Chans(:,:) ! List of aerosol detection channels + INTEGER(i_kind) :: N__Mean_Aerosol_Chans ! Boxcar averaging window width + REAL(r_kind), POINTER :: R__Aerosol_TBD(:,:) ! Aerosol detection thresholds + REAL(r_kind), POINTER :: R__coef_AOD(:,:) ! Coefficients for conversion to AOD + REAL(r_kind) :: R__Rank_Thres_Coeff(3) ! Coefficients to restrict rejections to affected channels + REAL(r_kind) :: R__Unclassified_Thres ! Rejection threshold for unclassified aerosol + REAL(r_kind) :: R__Land_Fraction_Thres ! Threshold for land fraction in FOV + END TYPE Aerosol_Detect_Type + + TYPE Cloud_Detect_Type + INTEGER(i_kind) :: M__Sensor ! Unique ID for sensor + INTEGER(i_kind) :: N__Num_Bands ! Number of channel bands + INTEGER(i_kind), POINTER :: N__GradChkInterval(:) ! Window width used in gradient calculation + INTEGER(i_kind), POINTER :: N__Band_Size(:) ! Number of channels in each band + INTEGER(i_kind), POINTER :: N__Bands(:,:) ! Channel lists + INTEGER(i_kind), POINTER :: N__Window_Width(:) ! Smoothing filter window widths per band + INTEGER(i_kind), POINTER :: N__Window_Bounds(:,:) ! Channels in the spectral window gradient check + INTEGER(i_kind), POINTER :: N__BandToUse(:) ! Band number assignment for each channel + LOGICAL :: L__Do_Quick_Exit ! On/off switch for the Quick Exit scenario + LOGICAL :: L__Do_CrossBand ! On/off switch for the cross-band method + REAL(r_kind), POINTER :: R__BT_Threshold(:) ! BT threshold for cloud contamination + REAL(r_kind), POINTER :: R__Grad_Threshold(:) ! Gradient threshold for cloud contamination + REAL(r_kind), POINTER :: R__Window_Grad_Threshold(:) ! Threshold for window gradient check in QE + + LOGICAL :: L__Do_Imager_Cloud_Detection ! On/off switch for the imager cloud detection + INTEGER(i_kind) :: N__Num_Imager_Chans ! No. of imager channels + INTEGER(i_kind) :: N__Num_Imager_Clusters ! No. of clusters to be expected + INTEGER(i_kind),POINTER :: N__Imager_Chans(:) ! List of imager channels + REAL(r_kind),POINTER :: R__Stddev_Threshold(:) ! St. Dev. threshold, one for each imager channel + REAL(r_kind) :: R__Coverage_Threshold ! Threshold for fractional coverage of a cluster + REAL(r_kind) :: R__FG_Departure_Threshold ! Threshold for imager FG departure + END TYPE Cloud_Detect_Type + + TYPE Land_Sensitivity_Type + INTEGER(r_kind) :: M__Sensor ! Unique ID for sensor + REAL(r_kind) :: R__Land_Fraction_Thres ! Threshold on land fraction + REAl(r_kind) :: R__Level_Thres ! Threshold on normalized channel height assignment + END TYPE Land_Sensitivity_Type + + TYPE Trace_Gas_Detect_Type + INTEGER(i_kind) :: M__Sensor ! Unique ID for sensor + INTEGER(i_kind) :: N__Num_Trace_Gas_Checks ! Number of trace gases to be checked + INTEGER(i_kind),POINTER :: N__Num_Tracer_Channels(:) ! Number of gas-sensitive channels + INTEGER(i_kind),POINTER :: N__Tracer_Channels(:,:) ! Gas-sensitive channels + INTEGER(i_kind),POINTER :: N__Num_Control_Channels(:) ! Number of control channels + INTEGER(i_kind),POINTER :: N__Control_Channels(:,:) ! Control channels + INTEGER(i_kind),POINTER :: N__Num_Flagged_Channels(:) ! Number of affected channels + INTEGER(i_kind),POINTER :: N__Flagged_Channels(:,:) ! Affected channels + REAL(r_kind),POINTER :: R__D_Obs_Threshold(:) ! Observed Tb difference threshold + REAL(r_kind),POINTER :: R__D_Dep_Threshold(:) ! Departure difference threshold + END TYPE Trace_Gas_Detect_Type + + + TYPE(Aerosol_Detect_Type) :: & + S__CADS_Setup_Aerosol(JP__Min_Sensor_Index:JP__Max_Sensor_Index) + + TYPE(Cloud_Detect_Type) :: & + S__CADS_Setup_Cloud(JP__Min_Sensor_Index:JP__Max_Sensor_Index) + + TYPE(Land_Sensitivity_Type) :: & + S__CADS_Setup_Land(JP__Min_Sensor_Index:JP__Max_Sensor_Index) + + TYPE(Trace_Gas_Detect_Type) :: & + S__CADS_Setup_Trace_Gas(JP__Min_Sensor_Index:JP__Max_Sensor_Index) + + +contains + +SUBROUTINE CADS_Abort(String) + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and MeteoFrance. +! +! Copyright 2020, EUMETSAT, All Rights Reserved. + +! *CADS_Abort* +! R. Eresmaa ECMWF 16/04/20 + +! * PURPOSE * +! ----------- +! Controlled abortion of running CADS when facing exceptions such as +! necessary input files missing or they are corrupt. + +! * INTERFACE * +! ------------- +! *CALL* * CADS_Abort()* from +! CADS_Main, CADS_Setup_Aerosol, CADS_Setup_Cloud, +! CADS_Setup_Land_Sensitivity, or CADS_Setup_Trace_Gas. + + IMPLICIT NONE + CHARACTER(LEN=*) :: String + + WRITE(*,*) String + STOP + +END SUBROUTINE CADS_Abort + +subroutine cloud_aerosol_detection( I__Sensor_ID, I__Num_Chans, I__Chan_ID, & + I__Min_Level, I__Max_Level, Z__BT_Obser, Z__BT_Model, Z__Chan_Height, K__Chan_ID_Imager, & + Z__Cluster_Fraction, Z__BT_in_Cluster, Z__BT_Overall_SDev, Z__BT_Model_Imager, & + I__Flag_Cloud, Z__Cloud_Level ) + +!$$$ subprogram documentation block +! . . . +! subprogram: cloud_aerosol_detection determine clear/cloudy profiles from hyperspectral IR instruments +! +! prgmmr: jung org: cimss date: 2022-10-17 +! +! abstract: determine if a profile is clear/cloudy. If cloudy, determine which channels are affected +! This subroutine is designed for infrared hyperspectral sounders. Current code supports AIRS, IASI and CrIS.a +! This subroutine is based on the Cloud and Aerosol Detection Software Version 3 developed within the context +! of the EUMETSAT and Met Office, UK, by one or more partners within the Numerical Weather Predicion's +! Satellite Application Facilities. A version of this code is operational at ECMWF. +! COPYRIGHT 2020, EUMETSAT, ALL RIGHTS RESERVED. +! +! program history log: +! 2022-10-17 jung Initial coding +! +! input argument list: +! I_Sensor_ID - internal sensor identification. +! I__Num_Chans - number of channels per obs +! I__Chan_ID - array of actual channel numbers +! Z__Longitude - FOV longitude +! Z__Latitude - FOV latitude +! Z__Land_Fraction - FOV land fraction +! I__Min_Level - model tropopause height (start of cloud detection) +! I__Max_Level - model top of boundary layer ( stop of cloud detection) +! Z__BT_Obser - observaton brightness temperature +! Z__BT_Model - model derived brightness temperature +! Z__Chan_Height - model derived height where an opaque cloud influences the radiance. +! also used to re-organize channels +! Z__Cloud_Level - Cloud height assignment +! +! output argument list: +! icloud_layer - model layer where cloud is detected +! +! attributes: +! language: f90 +! machine: +! +!$$$ end documentation block + + use kinds, only: i_kind, r_kind + implicit none + + integer(i_kind), intent(in ) :: I__Sensor_ID + integer(i_kind), intent(in ) :: I__Num_Chans + integer(i_kind),dimension(I__Num_Chans),intent(in ) :: I__Chan_ID + integer(i_kind), intent(in ) :: I__Min_Level !tropopause pressure + integer(i_kind), intent(in ) :: I__Max_Level !boundary layer pressure + real(r_kind), intent(in ) :: Z__BT_Obser(:) !Observation BT + real(r_kind), intent(in ) :: Z__BT_Model(:) !Model derived BT + real(r_kind), intent(in ) :: Z__Chan_Height(:) !Channel height assignmenta + integer(i_kind), intent(in ) :: K__Chan_ID_Imager(:) ! imager channel numbers + real(r_kind), intent(in ) :: Z__Cluster_Fraction(:) + real(r_kind), intent(in ) :: Z__BT_in_Cluster(:,:) + real(r_kind), intent(in ) :: Z__BT_Overall_SDev(:) + real(r_kind), intent(in ) :: Z__BT_Model_Imager(:) + real(r_kind), intent( out) :: Z__Cloud_Level ! cloud height assignment + integer(i_kind),dimension(I__Num_Chans),intent( out) :: I__Flag_Cloud ! cloud use flag + +! Interim prodcts + +! Diagnostics: percentages of positive detections +! Input/Output file management + + N__Num_Imager_Chans = S__CADS_Setup_Cloud(I__Sensor_ID) % N__Num_Imager_Chans + N__Num_Imager_Clusters = S__CADS_Setup_Cloud(I__Sensor_ID) % N__Num_Imager_Clusters + + CALL CADS_Detect_Cloud( I__Sensor_ID, I__Num_Chans, I__Chan_ID,I__Min_Level, I__Max_Level, N__Num_Imager_Chans, & + K__Chan_ID_Imager, N__Num_Imager_Clusters, I__Flag_Cloud, Z__BT_Obser, Z__BT_Model, Z__Chan_Height, & + Z__Cluster_Fraction, Z__BT_in_Cluster, Z__BT_Overall_SDev, Z__BT_Model_Imager, Z__Cloud_Level ) + +end subroutine cloud_aerosol_detection + +SUBROUTINE CADS_Setup_Cloud + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and Meteo France. +! +! Copyright 2020, EUMETSAT, All Rights Reserved. + + +! * Cloud detection setup * +! A. Collard ECMWF 01/02/06 + +! * PURPOSE * +! ----------- +! Initialise cloud detection parameters for advanced infrared sounders. + +! * INTERFACE * +! ------------- +! CADS_Setup_Cloud is called from CADS_Main. + +! * METHOD * +! ---------- +! Default values are assigned to the cloud detections setup structure. + +! MODIFICATIONS +! ------------- +! 01/02/06 A.Collard 1.0 Original code. +! 19/10/06 A.Collard 1.1 Use IASI 300 Subset Channels. +! 17/11/09 R.Eresmaa 1.2 Use IASI 366 Subset Channels. +! Include parameters of the Quick Exit / +! long-wave window gradient check parameters. +! 11/11/11 R.Eresmaa 1.3 Default channel list for AIRS bands 3-5 +! modified. +! Processing capability for CrIS added +! assuming a selection of 320 channels. +! 03/12/13 R,Eresmaa 2.0 Imager-assisted cloud detection added for +! IASI. +! Updated setup for CrIS. +! 19/01/15 R.Eresmaa 2.1 Remove unused variable specifications and +! switch aerosol detection on by default for +! AIRS and IASI. +! 10/11/15 R.Eresmaa 2.2 Changed instrument ID naming convention. +! Changed parameters of aerosol detection. +! 20/12/16 R.Eresmaa 2.3 Remove settings for aerosol detection. +! 05/02/19 R.Eresmaa 2.4 Explicit KIND specifications. +! Add HIRAS, GIIRS (IASING + IRS added earlier) +! 16/04/20 R.Eresmaa 3.0 Rename, tidy up. + + use kinds, only: i_kind, r_kind + use gsi_io, only: verbose + IMPLICIT NONE + +! Local variables + + CHARACTER(LEN=6) :: CL__InstrumentName + CHARACTER(LEN=20) :: CL__Cloud_Detection_File + + INTEGER(i_kind) :: J, J__Sensor ! Loop variables + INTEGER(i_kind) :: INIU1, IOS + +!----------------------- +! Namelist variables +!----------------------- + +! N.B. Max_Bands must be greater than 5 + INTEGER(i_kind), PARAMETER :: JP__Max_Bands = 8 + INTEGER(i_kind), PARAMETER :: JP__Max_Channels = 8461 + + INTEGER(i_kind) :: M__Sensor + INTEGER(i_kind) :: N__Num_Bands + INTEGER(i_kind) :: N__GradChkInterval(JP__Max_Bands) + INTEGER(i_kind) :: N__Band_Size(JP__Max_Bands) + INTEGER(i_kind) :: N__Bands(JP__Max_Channels,JP__Max_Bands) + INTEGER(i_kind) :: N__Window_Width(JP__Max_Bands) + INTEGER(i_kind) :: N__Window_Bounds(JP__Max_Bands,2) + REAL(r_kind) :: R__BT_Threshold(JP__Max_Bands) + REAL(r_kind) :: R__Grad_Threshold(JP__Max_Bands) + REAL(r_kind) :: R__Window_Grad_Threshold(JP__Max_Bands) + LOGICAL :: L__Do_Quick_Exit + LOGICAL :: L__Do_CrossBand + INTEGER(i_kind) :: N__BandToUse(JP__Max_Bands) + +! Imager-based cloud detection + LOGICAL :: L__Do_Imager_Cloud_Detection + INTEGER(i_kind) :: N__Num_Imager_Chans + INTEGER(i_kind) :: N__Num_Imager_Clusters + INTEGER(i_kind) :: N__Imager_Chans(JP__Max_Bands) + REAL(r_kind) :: R__Stddev_Threshold(JP__Max_Bands) + REAL(r_kind) :: R__Coverage_Threshold + REAL(r_kind) :: R__FG_Departure_Threshold + +! Namelist + NAMELIST / Cloud_Detect_Coeffs / M__Sensor, N__Num_Bands, & + N__Band_Size, N__Bands, N__Window_Width, N__Window_Bounds, & + N__GradChkInterval, R__BT_Threshold, R__Grad_Threshold, & + R__Window_Grad_Threshold, L__Do_Quick_Exit, & + L__Do_CrossBand, N__BandToUse, & + L__Do_Imager_Cloud_Detection, N__Num_Imager_Chans, & + N__Num_Imager_Clusters, N__Imager_Chans, & + R__Stddev_Threshold, R__Coverage_Threshold, & + R__FG_Departure_Threshold + +!============================================================================ +! Loop through sensors setting up cloud detection +!============================================================================ + + SensorLoop : DO J__Sensor = JP__Min_Sensor_Index, JP__Max_Sensor_Index + +! SELECT CASE (I__Sensor_ID) + SELECT CASE (J__Sensor) + + CASE(INST_ID_AIRS) + !==================== + ! Set up AIRS + !==================== + + CL__InstrumentName='AIRS' + CL__Cloud_Detection_File = 'AIRS_CLDDET.NL' + + N__Num_Bands = 5 + + N__Band_Size(:) = 0 + N__Band_Size(1:N__Num_Bands) =(/138, 36, 54, 23, 65 /) + + N__Bands(:,:)= 0 + + N__Bands(1:N__Band_Size(1),1) = & + (/ 1, 6, 7, 10, 11, 15, 16, 17, 20, 21, & + 22, 24, 27, 28, 30, 36, 39, 40, 42, 51, & + 52, 54, 55, 56, 59, 62, 63, 68, 69, 71, & + 72, 73, 74, 75, 76, 77, 78, 79, 80, 82, & + 83, 84, 86, 92, 93, 98, 99, 101, 104, 105, & + 108, 110, 111, 113, 116, 117, 123, 124, 128, 129, & + 138, 139, 144, 145, 150, 151, 156, 157, 159, 162, & + 165, 168, 169, 170, 172, 173, 174, 175, 177, 179, & + 180, 182, 185, 186, 190, 192, 193, 198, 201, 204, & + 207, 210, 213, 215, 216, 218, 221, 224, 226, 227, & + 232, 239, 248, 250, 251, 252, 253, 256, 257, 261, & + 262, 267, 272, 295, 299, 305, 308, 309, 310, & + 318, 321, 333, 338, 355, 362, 375, 475, & + 484, 497, 528, 587, 672, 787, 791, 843, 870, 914, & + 950 /) + + N__Bands(1:N__Band_Size(2),2) = & + (/ 1003, 1012, 1019, 1024, 1030, 1038, 1048, 1069, 1079, 1082, & + 1083, 1088, 1090, 1092, 1095, 1104, 1111, 1115, 1116, 1119, & + 1120, 1123, 1130, 1138, 1142, 1178, 1199, 1206, 1221, 1237, & + 1252, 1260, 1263, 1266, 1278, 1285 /) + + N__Bands(1:N__Band_Size(3),3) = & + (/ 1290, 1301, 1304, 1329, 1371, 1382, 1415, 1424, 1449, 1455, & + 1466, 1471, 1477, 1479, 1488, 1500, 1519, 1520, 1538, 1545, & + 1565, 1574, 1583, 1593, 1614, 1627, 1636, 1644, 1652, 1669, & + 1674, 1681, 1694, 1708, 1717, 1723, 1740, 1748, 1751, 1756, & + 1763, 1766, 1771, 1777, 1780, 1783, 1794, 1800, 1803, 1806, & + 1812, 1826, 1843, 1852 /) + + N__Bands(1:N__Band_Size(4),4) = & + (/ 1865, 1866, 1867, 1868, 1869, 1872, 1873, 1875, 1876, 1877, & + 1881, 1882, 1883, 1884, 1897, 1901, 1911, 1917, 1918, 1921, & + 1923, 1924, 1928 /) + + N__Bands(1:N__Band_Size(5),5) = & + (/ 1937, 1938, 1939, 1941, 1946, 1947, 1948, 1958, 1971, 1973, & + 1988, 1995, 2084, 2085, 2097, 2098, 2099, 2100, 2101, 2103, & + 2104, 2106, 2107, 2108, 2109, 2110, 2111, 2112, 2113, 2114, & + 2115, 2116, 2117, 2118, 2119, 2120, 2121, 2122, 2123, 2128, & + 2134, 2141, 2145, 2149, 2153, 2164, 2189, 2197, 2209, 2226, & + 2234, 2280, 2318, 2321, 2325, 2328, 2333, 2339, 2348, 2353, & + 2355, 2363, 2370, 2371, 2377 /) + + N__GradChkInterval(:) = 0 + N__GradChkInterval(1:N__Num_Bands) = (/ 5,5,5,5,5 /) + + N__Window_Width(:) = 0 + N__Window_Width(1:N__Num_Bands) = (/ 14,6,8,5,8 /) + + N__Window_Bounds(:,:) = 0 + N__Window_Bounds(1,1) = 475 + N__Window_Bounds(1,2) = 950 + + R__BT_Threshold(:) = 0.0_r_kind + R__BT_Threshold(1:N__Num_Bands) = (/ 0.43_r_kind, 0.5_r_kind, 0.5_r_kind, 0.5_r_kind, 0.5_r_kind/) + + R__Grad_Threshold(:) = 0.0_r_kind + R__Grad_Threshold(1:N__Num_Bands) = (/ 0.02_r_kind, 0.02_r_kind, 0.02_r_kind, 0.02_r_kind, 0.02_r_kind /) + + R__Window_Grad_Threshold(:) = 0.0_r_kind + R__Window_Grad_Threshold(1) = 0.4_r_kind + + L__Do_Quick_Exit = .TRUE. + + + ! This is cross-band: + + L__Do_CrossBand = .TRUE. + + N__BandToUse(:) = 0 + N__BandToUse(1:N__Num_Bands) = (/ 1,1,1,4,1 /) + + + ! This is the setup for imager cloud detection + + L__Do_Imager_Cloud_Detection = .FALSE. + + N__Num_Imager_Chans = 0 + N__Num_Imager_Clusters = 0 + N__Imager_Chans(:) = 0 + + R__Stddev_Threshold(:) = 0.0_r_kind + R__Coverage_Threshold = 0.0_r_kind + R__FG_Departure_Threshold = 0.0_r_kind + + + CASE(INST_ID_IASI) + !==================== + ! Set up IASI + !==================== + + CL__InstrumentName='IASI' + CL__Cloud_Detection_File = 'IASI_CLDDET.NL' + + N__Num_Bands = 5 + + N__Band_Size(:) = 0 + N__Band_Size(1:N__Num_Bands) =(/ 184, 15, 116, 4, 15 /) + + N__Bands(:,:)= 0 + + ! Use the "IASI 366" Subset + N__Bands(1:N__Band_Size(1),1) = & + (/ 16, 38, 49, 51, 55, 57, 59, 61, 63, 66, & + 70, 72, 74, 79, 81, 83, 85, 87, 89, 92, & + 95, 97, 99, 101, 104, 106, 109, 111, 113, 116, & + 119, 122, 125, 128, 131, 133, 135, 138, 141, 144, & + 146, 148, 151, 154, 157, 159, 161, 163, 165, 167, & + 170, 173, 176, 178, 179, 180, 183, 185, 187, 189, & + 191, 193, 195, 197, 199, 201, 203, 205, 207, 210, & + 212, 214, 217, 219, 222, 224, 226, 228, 230, 232, & + 234, 236, 239, 241, 242, 243, 246, 249, 252, 254, & + 256, 258, 260, 262, 265, 267, 269, 271, 272, 273, & + 275, 278, 280, 282, 284, 286, 288, 290, 292, 294, & + 296, 299, 301, 303, 306, 308, 310, 312, 314, 316, & + 318, 320, 323, 325, 327, 329, 331, 333, 335, 341, & + 347, 350, 352, 354, 356, 358, 360, 362, 364, 366, & + 369, 371, 373, 375, 377, 379, 381, 386, 389, 404, & + 407, 410, 414, 416, 426, 428, 432, 434, 445, 457, & + 515, 546, 552, 566, 571, 573, 646, 662, 668, 756, & + 867, 921, 1027, 1090, 1133, 1191, 1194, 1271, 1805, 1884, & + 1946, 1991, 2094, 2239 /) + + N__Bands(1:N__Band_Size(2),2) = & + (/ 1479, 1509, 1513, 1521, 1536, 1574, 1579, 1585, 1587, 1626, & + 1639, 1643, 1652, 1658, 1671 /) + + N__Bands(1:N__Band_Size(3),3) = & + (/ 2119, 2213, 2271, 2321, 2398, 2701, 2741, 2819, 2889, 2907, & + 2910, 2919, 2939, 2944, 2948, 2951, 2958, 2977, 2985, 2988, & + 2991, 2993, 3002, 3008, 3014, 3027, 3029, 3036, 3047, 3049, & + 3053, 3058, 3064, 3069, 3087, 3093, 3098, 3105, 3107, 3110, & + 3127, 3136, 3151, 3160, 3165, 3168, 3175, 3178, 3207, 3228, & + 3244, 3248, 3252, 3256, 3263, 3281, 3303, 3309, 3312, 3322, & + 3375, 3378, 3411, 3438, 3440, 3442, 3444, 3446, 3448, 3450, & + 3452, 3454, 3458, 3467, 3476, 3484, 3491, 3497, 3499, 3504, & + 3506, 3509, 3518, 3527, 3555, 3575, 3577, 3580, 3582, 3586, & + 3589, 3599, 3653, 3658, 3661, 4032, 5368, 5371, 5379, 5381, & + 5383, 5397, 5399, 5401, 5403, 5405, 5455, 5480, 5483, 5485, & + 5492, 5502, 5507, 5509, 5517, 5558 /) + + N__Bands(1:N__Band_Size(4),4) = & + (/ 5988, 5992, 5994, 6003 /) + + N__Bands(1:N__Band_Size(5),5) = & + (/ 6982, 6985, 6987, 6989, 6991, 6993, 6995, 6997, 7267, 7269, & + 7424, 7426, 7428, 7885, 8007 /) + + N__GradChkInterval(:) = 0 + N__GradChkInterval(1:N__Num_Bands) = (/12,5,5,5,5 /) + + N__Window_Width(:) = 0 + N__Window_Width(1:N__Num_Bands) = (/ 10,6,8,5,8 /) + + N__Window_Bounds(:,:) = 0 + N__Window_Bounds(1,1) = 573 + N__Window_Bounds(1,2) = 2239 + + R__BT_Threshold(:) = 0.0_r_kind + R__BT_Threshold(1:N__Num_Bands) = (/ 0.5_r_kind, 0.5_r_kind, 0.5_r_kind, 0.5_r_kind, 0.5_r_kind /) + + R__Grad_Threshold(:) = 0.0_r_kind + R__Grad_Threshold(1:N__Num_Bands) = (/ 0.02_r_kind, 0.02_r_kind, 0.02_r_kind, 0.02_r_kind, 0.02_r_kind /) + + R__Window_Grad_Threshold(:) = 0.0_r_kind + R__Window_Grad_Threshold(1) = 0.4_r_kind + + L__Do_Quick_Exit = .TRUE. + + + ! This is cross-band: + + L__Do_CrossBand = .TRUE. + + N__BandToUse(:) = 0 + N__BandToUse(1:N__Num_Bands) = (/ 1,1,1,1,1 /) + + + ! This is the setup for imager cloud detection + + L__Do_Imager_Cloud_Detection = .TRUE. + + N__Num_Imager_Chans = 2 + N__Num_Imager_Clusters = 7 + + N__Imager_Chans(1:N__Num_Imager_Chans) = (/ 2, 3 /) + + R__Stddev_Threshold(1:N__Num_Imager_Chans) = (/ 0.75_r_kind, 0.80_r_kind /) + + R__Coverage_Threshold = 0.03_r_kind + R__FG_Departure_Threshold = 1.0_r_kind + + + CASE(INST_ID_CRIS) + !==================== + ! Set up CRIS + !==================== + + CL__InstrumentName='CRIS' + CL__Cloud_Detection_File = 'CRIS_CLDDET.NL' + + N__Num_Bands = 5 + + N__Band_Size(:) = 0 + + N__Band_Size(1:N__Num_Bands) =(/ 137, 123, 76, 12, 6 /) + + N__Bands(:,:)= 0 + + ! Use the "CRIS 300" Subset + N__Bands(1:N__Band_Size(1),1) = & + (/ 1, 5, 9, 13, 17, 18, 19, 20, 21, 22, & + 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, & + 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, & + 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, & + 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, & + 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, & + 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, & + 83, 84, 85, 86, 87, 88, 91, 92, 93, 94, & + 95, 96, 97, 99, 101, 105, 107, 109, 111, 113, & + 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, & + 125, 133, 135, 137, 139, 141, 144, 147, 161, 173, & + 177, 181, 185, 195, 210, 221, 225, 229, 249, 257, & + 269, 273, 293, 301, 317, 333, 349, 369, 409, 433, & + 457, 481, 501, 549, 701, 705, 709 /) + + N__Bands(1:N__Band_Size(2),2) = & + (/ 3, 6, 7, 8, 10, 12, 14, 15, 16, 89, & + 90, 102, 103, 104, 106, 108, 110, 114, 126, 127, & + 129, 132, 134, 138, 140, 143, 145, 146, 148, 149, & + 150, 151, 153, 155, 156, 157, 158, 159, 162, 163, & + 164, 165, 166, 169, 170, 171, 172, 175, 180, 189, & + 200, 201, 205, 206, 214, 217, 218, 226, 228, 230, & + 231, 233, 236, 237, 240, 241, 245, 248, 252, 264, & + 265, 281, 285, 297, 324, 327, 361, 378, 389, 392, & + 400, 473, 493, 500, 503, 511, 527, 528, 529, 530, & + 531, 534, 538, 542, 544, 545, 547, 550, 553, 555, & + 590, 594, 598, 602, 606, 610, 614, 618, 622, 626, & + 645, 649, 653, 657, 661, 665, 685, 702, 703, 704, & + 706, 707, 713 /) + + N__Bands(1:N__Band_Size(3),3) = & + (/ 717, 725, 728, 729, 730, 731, 732, 733, 734, 735, & + 736, 741, 749, 757, 765, 773, 781, 789, 794, 797, & + 805, 806, 815, 822, 829, 839, 845, 853, 861, 868, & + 869, 872, 877, 885, 887, 893, 898, 900, 909, 912, & + 915, 917, 921, 929, 933, 941, 949, 957, 963, 965, & + 973, 975, 978, 981, 989, 991, 993, 996, 1005, 1014, & + 1025, 1029, 1037, 1042, 1053, 1061, 1073, 1077, 1085, 1093, & + 1101, 1109, 1117, 1125, 1133, 1141 /) + + N__Bands(1:N__Band_Size(4),4) = & + (/ 1149, 1157, 1164, 1165, 1173, 1181, 1189, 1197, 1205, 1213, & + 1221, 1251 /) + + N__Bands(1:N__Band_Size(5),5) = & + (/ 1189, 1197, 1205, 1213, 1221, 1251 /) + + + N__GradChkInterval(:) = 0 + N__GradChkInterval(1:N__Num_Bands) = (/ 5,5,5,3,3 /) + + N__Window_Width(:) = 0 + N__Window_Width(1:N__Num_Bands) = (/ 6,6,8,3,3 /) + + N__Window_Bounds(:,:) = 0 + N__Window_Bounds(1,1) = 229 + N__Window_Bounds(1,2) = 549 + + R__BT_Threshold(:) = 0.0_r_kind + R__BT_Threshold(1:N__Num_Bands) = (/ 0.5_r_kind, 0.5_r_kind, 0.5_r_kind, 0.5_r_kind, 0.5_r_kind /) + + R__Grad_Threshold(:) = 0.0_r_kind + R__Grad_Threshold(1:N__Num_Bands) = (/ 0.02_r_kind, 0.02_r_kind, 0.02_r_kind, 0.02_r_kind, 0.02_r_kind /) + + R__Window_Grad_Threshold(:) = 0.0_r_kind + R__Window_Grad_Threshold(1) = 0.4_r_kind + + L__Do_Quick_Exit = .TRUE. + + + ! This is cross-band: + + L__Do_CrossBand = .TRUE. + + N__BandToUse(:) = 0 + N__BandToUse(1:N__Num_Bands) = (/ 1,1,1,1,1 /) + + + ! This is the setup for imager cloud detection + + L__Do_Imager_Cloud_Detection = .FALSE. + + N__Num_Imager_Chans = 0 + N__Num_Imager_Clusters = 0 + N__Imager_Chans(:) = 0 + + R__Stddev_Threshold(:) = 0.0_r_kind + R__Coverage_Threshold = 0.0_r_kind + R__FG_Departure_Threshold = 0.0_r_kind + + + CASE(INST_ID_IRS) + !==================== + ! Set up IRS + !==================== + + CL__InstrumentName='IRS' + CL__Cloud_Detection_File = 'IRS_CLDDET.NL' + + N__Num_Bands = 1 + + N__Band_Size(:) = 0 + + N__Band_Size(1:N__Num_Bands) =(/ 138 /) + + N__Bands(:,:)= 0 + + N__Bands(1:N__Band_Size(1),1) = & + (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, & + 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, & + 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, & + 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, & + 41, 42, 43, 44, 45, 46, 48, 53, 54, 55, & + 56, 57, 58, 60, 61, 62, 63, 65, 70, 74, & + 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, & + 85, 86, 87, 89, 90, 91, 92, 93, 94, 95, & + 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, & + 106, 107, 108, 109, 118, 119, 131, 145, 163, 169, & + 177, 180, 190, 195, 199, 209, 215, 221, 231, 237, & + 252, 262, 268, 281, 289, 298, 312, 322, 328, 341, & + 347, 359, 375, 384, 390, 404, 412, 421, 648, 656, & + 667, 678, 686, 692, 709, 750, 792, 808 /) + + N__GradChkInterval(:) = 0 + N__GradChkInterval(1:N__Num_Bands) = (/ 12 /) + + N__Window_Width(:) = 0 + N__Window_Width(1:N__Num_Bands) = (/ 10 /) + + N__Window_Bounds(:,:) = 0 + N__Window_Bounds(1,1) = 131 + N__Window_Bounds(1,2) = 808 + + R__BT_Threshold(:) = 0.0_r_kind + R__BT_Threshold(1:N__Num_Bands) = (/ 0.4_r_kind /) + + R__Grad_Threshold(:) = 0.0_r_kind + R__Grad_Threshold(1:N__Num_Bands) = (/ 0.02_r_kind /) + + R__Window_Grad_Threshold(:) = 0.0_r_kind + R__Window_Grad_Threshold(1) = 0.4_r_kind + + L__Do_Quick_Exit = .TRUE. + + + ! This is cross-band: + + L__Do_CrossBand = .TRUE. + + N__BandToUse(:) = 0 + N__BandToUse(1:N__Num_Bands) = (/ 1 /) + + + ! This is the setup for imager cloud detection + + L__Do_Imager_Cloud_Detection = .FALSE. + + N__Num_Imager_Chans = 0 + N__Num_Imager_Clusters = 0 + N__Imager_Chans(:) = 0 + + R__Stddev_Threshold(:) = 0.0_r_kind + R__Coverage_Threshold = 0.0_r_kind + R__FG_Departure_Threshold = 0.0_r_kind + + + CASE(INST_ID_IASING) + !==================== + ! Set up IASING + !==================== + + CL__InstrumentName='IASING' + CL__Cloud_Detection_File = 'IASING_CLDDET.NL' + + N__Num_Bands = 1 + + N__Band_Size(:) = 0 + + N__Band_Size(1:N__Num_Bands) =(/ 254 /) + + N__Bands(:,:)= 0 + + N__Bands(1:N__Band_Size(1),1) = & + (/ 31, 75, 97, 101, 109, 113, 117, 121, 125, 131, & + 139, 143, 147, 157, 161, 165, 169, 173, 177, 183, & + 189, 193, 197, 201, 207, 211, 217, 221, 225, 231, & + 237, 243, 249, 255, 261, 265, 269, 275, 281, 287, & + 291, 295, 301, 307, 313, 317, 321, 325, 329, 333, & + 339, 345, 351, 355, 357, 359, 365, 369, 373, 377, & + 381, 385, 389, 393, 397, 401, 403, 405, 407, 409, & + 411, 413, 415, 417, 419, 421, 423, 425, 427, 429, & + 431, 433, 435, 437, 439, 441, 443, 445, 447, 449, & + 451, 453, 455, 457, 459, 461, 463, 465, 467, 469, & + 471, 473, 475, 477, 479, 481, 483, 485, 487, 489, & + 491, 493, 495, 497, 499, 501, 503, 505, 507, 509, & + 511, 513, 515, 517, 519, 521, 523, 525, 527, 529, & + 531, 533, 535, 537, 539, 541, 543, 545, 547, 549, & + 551, 553, 555, 557, 559, 561, 563, 565, 567, 569, & + 571, 573, 575, 577, 579, 581, 583, 585, 587, 589, & + 591, 593, 595, 597, 601, 603, 605, 607, 609, 611, & + 613, 615, 617, 619, 621, 623, 625, 627, 629, 631, & + 633, 635, 637, 639, 641, 643, 645, 647, 649, 651, & + 653, 655, 657, 659, 661, 663, 665, 667, 669, 681, & + 693, 699, 703, 707, 711, 715, 719, 723, 727, 731, & + 737, 741, 745, 749, 753, 757, 761, 771, 777, 807, & + 813, 819, 827, 831, 851, 855, 863, 867, 889, 913, & + 1029, 1091, 1103, 1131, 1141, 1145, 1291, 1323, 1335, 1511, & + 1733, 1841, 2053, 2179, 2265, 2381, 2387, 2541, 3609, 3767, & + 3891, 3981, 4187, 4477 /) + + N__GradChkInterval(:) = 0 + N__GradChkInterval(1:N__Num_Bands) = (/ 25 /) + + N__Window_Width(:) = 0 + N__Window_Width(1:N__Num_Bands) = (/ 20 /) + + N__Window_Bounds(:,:) = 0 + N__Window_Bounds(1,1) = 1145 + N__Window_Bounds(1,2) = 4477 + + R__BT_Threshold(:) = 0.0_r_kind + R__BT_Threshold(1:N__Num_Bands) = (/ 0.27_r_kind /) + + R__Grad_Threshold(:) = 0.0_r_kind + R__Grad_Threshold(1:N__Num_Bands) = (/ 0.02_r_kind /) + + R__Window_Grad_Threshold(:) = 0.0_r_kind + R__Window_Grad_Threshold(1) = 0.4_r_kind + + L__Do_Quick_Exit = .TRUE. + + + ! This is cross-band: + + L__Do_CrossBand = .TRUE. + + N__BandToUse(:) = 0 + N__BandToUse(1:N__Num_Bands) = (/ 1 /) + + ! This is the setup for imager cloud detection + + L__Do_Imager_Cloud_Detection = .FALSE. + + N__Num_Imager_Chans = 0 + N__Num_Imager_Clusters = 0 + N__Imager_Chans(:) = 0 + + R__Stddev_Threshold(:) = 0.0_r_kind + R__Coverage_Threshold = 0.0_r_kind + R__FG_Departure_Threshold = 0.0_r_kind + + + END SELECT + + !------------------------------------------------------------------ + ! Open and read file containing cloud detection setup for the + ! current instrument + !------------------------------------------------------------------ + + INIU1=107 + OPEN(INIU1,STATUS='OLD',FORM='FORMATTED', & + FILE=TRIM(CL__Cloud_Detection_File), IOSTAT=IOS) + IF (IOS == 0) THEN + READ(INIU1,nml=Cloud_Detect_Coeffs,IOSTAT=IOS) + IF (IOS == 0) THEN + if ( verbose ) WRITE(*,'(3X,A)') TRIM(CL__InstrumentName) // & + ' CLOUD DETECTION FILE READ OK' + ELSE + CALL CADS_Abort('PROBLEM READING '//TRIM(CL__InstrumentName)//& + 'CLOUD DETECTION FILE') + ENDIF + CLOSE(INIU1) + ELSE + if ( verbose ) WRITE(*,'(3X,A)') 'NO '//TRIM(CL__InstrumentName) // & + ' CLOUD DETECTION FILE : Using Default Values' + ENDIF + + IF (MAXVAL(N__Band_Size(:)) > JP__Max_Channels) & + CALL CADS_Abort('Too many channels specified in cloud '//& + 'detection - increase JP__Max_Channels') + + + M__Sensor = J__SENSOR + + !------------------------------------------------------------------ + ! Set up the S__CADS_Setup_Cloud structure for current sensor + !------------------------------------------------------------------ + + S__CADS_Setup_Cloud(J__SENSOR) % M__SENSOR = M__Sensor + + S__CADS_Setup_Cloud(J__SENSOR) % N__Num_Bands = N__Num_Bands + + ALLOCATE( S__CADS_Setup_Cloud(J__SENSOR) % N__Band_Size(N__Num_Bands) ) + + S__CADS_Setup_Cloud(J__SENSOR) % N__Band_Size(:) = & + N__Band_Size(1:N__Num_Bands) + + ALLOCATE(S__CADS_Setup_Cloud(J__SENSOR) % N__Bands & + (MAXVAL(N__Band_Size(:)), N__Num_Bands)) + + S__CADS_Setup_Cloud(J__SENSOR) % N__Bands(:,:) = 0 + + DO J = 1, N__Num_Bands + S__CADS_Setup_Cloud(J__SENSOR) % N__Bands(1:N__Band_Size(J),J) = & + N__Bands(1:N__Band_Size(J),J) + ENDDO + + ALLOCATE( S__CADS_Setup_Cloud(J__SENSOR) % N__Window_Width(N__Num_Bands) ) + + S__CADS_Setup_Cloud(J__SENSOR) % N__Window_Width(:) = & + N__Window_Width(1:N__Num_Bands) + + ALLOCATE( S__CADS_Setup_Cloud(J__SENSOR) % R__BT_Threshold(N__Num_Bands) ) + S__CADS_Setup_Cloud(J__SENSOR) % R__BT_Threshold(:) = & + R__BT_Threshold(1:N__Num_Bands) + + ALLOCATE(S__CADS_Setup_Cloud(J__SENSOR) % R__Grad_Threshold(N__Num_Bands)) + S__CADS_Setup_Cloud(J__SENSOR) % R__Grad_Threshold(:) = & + R__Grad_Threshold(1:N__Num_Bands) + + ALLOCATE(S__CADS_Setup_Cloud(J__SENSOR) % & + R__Window_Grad_Threshold(N__Num_Bands)) + + S__CADS_Setup_Cloud(J__SENSOR) % R__Window_Grad_Threshold(:) = & + R__Window_Grad_Threshold(1:N__Num_Bands) + + ALLOCATE(S__CADS_Setup_Cloud(J__SENSOR) % N__GradChkInterval(N__Num_Bands)) + S__CADS_Setup_Cloud(J__SENSOR) % N__GradChkInterval(:) = & + N__GradChkInterval(1:N__Num_Bands) + + ALLOCATE(S__CADS_Setup_Cloud(J__SENSOR) % N__Window_Bounds(N__Num_Bands,2)) + S__CADS_Setup_Cloud(J__SENSOR) % N__Window_Bounds(:,:) = & + N__Window_Bounds(1:N__Num_Bands,:) + + S__CADS_Setup_Cloud(J__SENSOR) % L__Do_Quick_Exit = L__Do_Quick_Exit + + + !------------- + ! Cross Band + !------------- + + S__CADS_Setup_Cloud(J__SENSOR) % L__Do_CrossBand = L__Do_CrossBand + + ALLOCATE( S__CADS_Setup_Cloud(J__SENSOR) % N__BandToUse(N__Num_Bands) ) + S__CADS_Setup_Cloud(J__SENSOR) % N__BandToUse(:) = & + N__BandToUse(1:N__Num_Bands) + + + !------------- + ! Imager cloud detection + !------------- + + S__CADS_Setup_Cloud(J__SENSOR) % L__Do_Imager_Cloud_Detection = & + L__Do_Imager_Cloud_Detection + + S__CADS_Setup_Cloud(J__SENSOR) % N__Num_Imager_Chans = & + N__Num_Imager_Chans + + S__CADS_Setup_Cloud(J__SENSOR) % N__Num_Imager_Clusters = & + N__Num_Imager_Clusters + + ALLOCATE( S__CADS_Setup_Cloud(J__SENSOR) % & + N__Imager_Chans(N__Num_Imager_Chans)) + S__CADS_Setup_Cloud(J__SENSOR) % N__Imager_Chans(:) = & + N__Imager_Chans(1:N__Num_Imager_Chans) + + ALLOCATE( S__CADS_Setup_Cloud(J__SENSOR) % & + R__Stddev_Threshold(N__Num_Imager_Chans)) + S__CADS_Setup_Cloud(J__SENSOR) % R__Stddev_Threshold(:) = & + R__Stddev_Threshold(1:N__Num_Imager_Chans) + + S__CADS_Setup_Cloud(J__SENSOR) % R__Coverage_Threshold = & + R__Coverage_Threshold + + S__CADS_Setup_Cloud(J__SENSOR) % R__FG_Departure_Threshold = & + R__FG_Departure_Threshold + + ENDDO SensorLoop + +END SUBROUTINE CADS_SETUP_CLOUD + +SUBROUTINE CADS_Detect_Cloud( K__Sensor, K__NChans, K__ChanID, K__Minlev, K__Maxlev, & + K__Num_Imager_Chans, K__Chan_ID_Imager, K__Num_Imager_Clusters, & + K__Cloud_Flag, P__ObsBTs, P__ModelBTs, P__Chan_Level, P__Cluster_Fraction,& + P__BT_in_Cluster, P__BT_Overall_SDev, P__BT_Model_Imager, Z__Cloud_Level ) + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and MeteoFrance. + +! Copyright 2020, EUMETSAT, All Rights Reserved. + +! * CADS_Detect_Cloud * +! Phil Watts ECMWF 21/01/02 + +! * PURPOSE * +! ----------- +! Flag the presence or otherwise of cloud contamination in AIRS/IASI +! channels using a rank-sorted/model difference method. Currently +! only a digital filter is supported. + + +! * INTERFACE * +! ------------- +! *CALL* * CADS_Detect_Cloud( )* (from CADS_Main) +! WHERE K__Sensor : Satellite sensor (AIRS/IASI/CrIS) +! K__NChans : Number of channels +! K__ChanID : Channel indices of input channels +! K__Minlev : Highest allowed starting point for the cloud search +! K__Maxlev : Lowest allowed starting point in the initial cloud search +! K__Num_Imager_Chans : Number of collocated imager channels +! K__Chan_ID_Imager : Collocated imager channel indices +! K__Num_Imager_Clusters : Number of collocated clusters +! K__Cloud_Flag : Cloud flag by channel; 0=clear, 1=cloudy +! P__ObsBTs : Potentially cloud-affected observed BTs +! P__ModelBTs : Clear background brightness temperatures (BTs) +! P__Chan_Level : Channel height assignments +! P__Cluster_Fraction : Fractional coverage of each cluster within FOV +! P__BT_in_Cluster : Cluster-mean brightness temperature (BT) on each channel +! P__BT_Overall_SDev : Overall BT standard deviation on each channel +! P__BT_Model_Imager : Forward-modelled BT on each channel +! Z__Cloud_Level : Cloud height assignment + +! * EXTERNALS * +! ------------- +! CADS_Detect_Cloud_Imager, CADS_Detect_Cloud_Heapsort, +! CADS_Detect_Cloud_Smooth, CADS_Detect_Cloud_Scenario, +! CADS_Detect_Cloud_Separator + +! * MODIFICATIONS * +! ----------------- +! A.Collard 1.0 01/02/06 Original export version +! A.Collard 1.0.1 03/05/06 Allow for missing channels +! A.Collard 1.0.2 04/05/06 Allow cross-band cloud detection +! A.Collard 1.0.3 15/01/07 Initialise with automatic cross-band for +! all channels from band 1 for IASI +! R.Eresmaa 1.1 17/11/09 Include parameters of the Quick Exit / +! long-wave window gradient check. +! Pass K__Chan_Low to CF_DIGITAL to allow +! detecting cirrus in case of compensating +! humidity bg error in PBL. +! R.Eresmaa 1.2 11/11/11 Modify the cross-band option to be based +! on the lowest clear channel rather than +! on the highest cloud-contaminated one +! R.Eresmaa 2.0 27/11/13 Add input cloud flag based on collocated +! imager data +! R.Eresmaa 2.1 13/01/15 Make array size specifications implicit. +! R.Eresmaa 2.2 10/11/15 Instrument ID naming convention made +! consistent with RTTOV. +! Changed setting of the aerosol flag. +! R.Eresmaa 2.2.1 13/11/15 Don't allow flagging missing channels clear +! through the cross-band option. +! R.Eresmaa 2.3 20/12/16 Remove the call to aerosol detection. +! R.Eresmaa 2.4 05/02/19 Explicit KIND specifications. +! R.Eresmaa 3.0 16/04/20 Move the call to imager-based detection here. + + use kinds, only: i_kind, r_kind + use gsi_io, only: verbose + IMPLICIT NONE + +!* 0.1 Global arrays + INTEGER(i_kind), INTENT(IN) :: K__Sensor ! Sensor + INTEGER(i_kind), INTENT(IN) :: K__NChans ! No. of channels + INTEGER(i_kind), INTENT(IN) :: K__ChanID(:) ! Channel IDs + INTEGER(i_kind), INTENT(IN) :: K__Minlev ! Highest starting point for cloud search + INTEGER(i_kind), INTENT(IN) :: K__Maxlev ! Lowest starting point in the initial search + INTEGER(i_kind), INTENT(IN) :: K__Num_Imager_Chans ! No. of imager channels + INTEGER(i_kind), INTENT(IN) :: K__Chan_ID_Imager(:) ! Imager channel IDs + INTEGER(i_kind), INTENT(IN) :: K__Num_Imager_Clusters ! No. of imager clusters + INTEGER(i_kind), INTENT(OUT) :: K__Cloud_Flag(:) ! Output cloud flags + REAL(r_kind), INTENT(IN) :: P__ObsBTs(:) ! Observed BTs + REAL(r_kind), INTENT(IN) :: P__ModelBTs(:) ! Model clear BTs + REAL(r_kind), INTENT(IN) :: P__Chan_Level(:) ! Channel height assignments + REAL(r_kind), INTENT(IN) :: P__Cluster_Fraction(:) ! Cluster coverages + REAL(r_kind), INTENT(IN) :: P__BT_in_Cluster(:,:) ! Mean BT in cluster / channel + REAL(r_kind), INTENT(IN) :: P__BT_Overall_Sdev(:) ! St.Dev of imager BT in FOV + REAL(r_kind), INTENT(IN) :: P__BT_Model_Imager(:) ! Model-based estimate of imager BT + REAL(r_kind), INTENT(OUT) :: Z__Cloud_Level ! Cloud hight assignment + +!* 0.2 local variables + INTEGER(i_kind) :: IST,ICOUNT,J,I_K,JBAND,JBAND2 + INTEGER(i_kind) :: I__Imager_Flag ! Preliminary cloud flag from collocated imager data + +!* 0.3 Local variables - band splitting details + INTEGER(i_kind), POINTER :: I__Bands(:,:) ! Channel bands + INTEGER(i_kind), POINTER :: I__Band_Size(:) ! Number of channels per band + INTEGER(i_kind), POINTER :: I__BandToUse(:) ! Cross-band definitions + INTEGER(i_kind) :: I__Num_Bands ! Number of bands + INTEGER(i_kind) :: I__NumFoundChans ! Number of usable channels + INTEGER(i_kind) :: I__BandNumber(K__NChans) ! Channel band indicator + INTEGER(i_kind) :: I__WindowBounds(2) ! Boundary of window + INTEGER(i_kind) :: I__Window_Chans(2) ! Boundary of long-wave window + INTEGER(i_kind), ALLOCATABLE :: I__INDEX(:) ! Channel ranking within a band + INTEGER(i_kind), ALLOCATABLE :: IDCHAN(:) ! Overall channel ranking + INTEGER(i_kind), ALLOCATABLE :: I__Cloud_Flag(:) ! Rank-sorted output cloud flags + INTEGER(i_kind) :: I__Scenario_Index ! 1--Quick Exit, 2--Warm Start, 3--Cold Start + INTEGER(i_kind) :: I__Start_Channel ! Final starting channel in the cloud search + + LOGICAL :: LL__Do_CrossBand + +! Input array projections (handling one detection band at a time) + REAL(r_kind), ALLOCATABLE :: Z__DBT(:) ! Original departures + REAL(r_kind), ALLOCATABLE :: Z__Smooth_DBT(:) ! Smoothed departures + REAL(r_kind), ALLOCATABLE :: Z__LEVEL(:) ! Channel height assignments + +!* 0.4 Local variables - digital filter parameters + INTEGER(i_kind) :: I__CHAN_HIGH ! Channel at K__Minlev + INTEGER(i_kind) :: I__CHAN_LOW ! Channel at K__Maxlev + INTEGER(i_kind) :: I__FirstCloudyChannel ! Highest cloud-affected channel + INTEGER(i_kind) :: I__LastClearChannel ! Lowest clear channel + INTEGER(i_kind),POINTER :: I__Window_Width(:) ! Box-car filter width + INTEGER(i_kind),POINTER :: I__GradChkInterval(:) ! Gradient-check interval + +!====================================================================== + + +! Get correct processing parameters for this sensor: + I__Num_Bands = S__CADS_Setup_Cloud(K__Sensor) % N__Num_Bands + I__Band_Size => S__CADS_Setup_Cloud(K__Sensor) % N__Band_Size + I__Bands => S__CADS_Setup_Cloud(K__Sensor) % N__Bands + I__Window_Width => S__CADS_Setup_Cloud(K__Sensor) % N__Window_Width + I__BandToUse => S__CADS_Setup_Cloud(K__Sensor) % N__BandToUse + LL__Do_CrossBand = S__CADS_Setup_Cloud(K__Sensor) % L__Do_CrossBand + I__GradChkInterval => S__CADS_Setup_Cloud(K__Sensor) % N__GradChkInterval + + +! Initialise + K__Cloud_Flag(:)=1 ! intialise ALL channels to cloudy + + +! Imager-based cloud detection + I__Imager_Flag=0 ! Default assumption: no cloud affecting collocated imager data + CALL CADS_Detect_Cloud_Imager( K__Sensor, K__Num_Imager_Chans, K__Chan_ID_Imager, K__Num_Imager_Clusters, & + I__Imager_Flag, P__Cluster_Fraction, P__BT_in_Cluster, P__BT_Overall_SDev, P__BT_Model_Imager ) + +! If using cross-band, set up an array indicating which channels correspond +! to which bands in K__ChanID + IF (LL__Do_CrossBand) THEN + I__BandNumber(:)=-1 ! Initialise + DO JBAND = 1, I__Num_Bands + DO I_K=1,K__NChans + IF (ANY(I__BANDS(:,JBAND) == K__ChanID(I_K))) & + I__BandNumber(I_K)=JBand + ENDDO + ENDDO + ENDIF + + +!1 Loop over bands + Band_Loop: DO JBAND = 1, I__Num_Bands + + ! Don't bother doing the cloud detection if we're just going to use + ! the results from another band anyway: + IF (LL__Do_CrossBand) THEN + IF (.NOT.(ANY(I__BandToUse(:) == JBAND))) CYCLE + ENDIF + + ALLOCATE (Z__DBT(I__Band_Size(JBAND))) + Z__DBT(:) = 0.0_r_kind + + ALLOCATE (Z__LEVEL(I__Band_Size(JBAND))) + Z__LEVEL(:) = REAL(K__Maxlev) + + ALLOCATE (I__Cloud_Flag(I__Band_Size(JBAND))) + ALLOCATE (I__INDEX(I__Band_Size(JBAND))) + + ALLOCATE (IDCHAN(I__Band_Size(JBAND))) + IDCHAN(:) = 1 + + + I__WindowBounds(:) = & + S__CADS_Setup_Cloud(K__Sensor) % N__Window_Bounds(JBand,:) + +!1.1 find channels within current band -------------------------------------- + I__NumFoundChans = 0 + I__Window_Chans(:) = -1 + + DO J=1,I__Band_Size(JBAND) + DO I_K=1,K__NChans + IF (K__ChanID(I_K) == I__BANDS(J,JBAND)) THEN +! IF (P__ObsBTs(I_K) < 0. .OR. P__ModelBTs(I_K) < 0.) CYCLE + IF (P__ObsBTs(I_K) < 60.0_r_kind .OR. P__ModelBTs(I_K) < 60.0_r_kind) CYCLE ! Missing channels are set to 50.0K + I__NumFoundChans = I__NumFoundChans + 1 + Z__DBT(I__NumFoundChans)=P__ObsBTs(I_K)-P__ModelBTs(I_K) + Z__LEVEL(I__NumFoundChans)=P__Chan_Level(I_K) + I__INDEX(I__NumFoundChans)=I__NumFoundChans + IDCHAN(I__NumFoundChans)=I_K + IF (K__ChanID(I_K) == I__WindowBounds(1)) & + I__Window_Chans(1) = I__NumFoundChans + IF (K__ChanID(I_K) == I__WindowBounds(2)) & + I__Window_Chans(2) = I__NumFoundChans + ENDIF + ENDDO + ENDDO + IF ( I__NumFoundChans == 0 ) THEN + if (verbose) WRITE(*,*) & + '**CADS_Detect_Cloud - WARNING: ' // & + 'CHANNELS NOT FOUND CYCLING BAND: **', JBAND + IF (ALLOCATED(Z__DBT)) DEALLOCATE (Z__DBT) + IF (ALLOCATED(Z__LEVEL)) DEALLOCATE (Z__LEVEL) + IF (ALLOCATED(I__Cloud_Flag)) DEALLOCATE (I__Cloud_Flag) + IF (ALLOCATED(I__INDEX)) DEALLOCATE (I__INDEX) + IF (ALLOCATED(IDCHAN)) DEALLOCATE (IDCHAN) + CYCLE Band_Loop + ENDIF + +!---------------------------------------------------------------------------- + IST=0 + ICOUNT=I__NumFoundChans + I__Cloud_Flag(:)=1 + +!2. Sort according to channel height assignments + CALL CADS_Detect_Cloud_Heapsort(I__NumFoundChans,Z__Level,I__Index) + +!2.1 Find I__CHAN_LOW - lowest channel considered in the initial cloud search + J=1 + DO WHILE (J < I__NumFoundChans .AND. Z__Level(I__Index(J)) < REAL(K__Maxlev)) + J=J+1 + ENDDO + + IF (J == I__NumFoundChans) THEN + I__CHAN_LOW = I__NumFoundChans-1 + ELSE + I__CHAN_LOW = J + ENDIF + IF(I__CHAN_LOW <= 1)I__CHAN_LOW=1 + +!2.1a Find I__CHAN_HIGH - highest allowed channel for starting the cloud search + J=1 + DO WHILE (J < I__NumFoundChans .AND. Z__Level(I__Index(J)) < REAL(K__Minlev)) + J=J+1 + ENDDO + I__CHAN_HIGH=J + + +! Smoothing + ALLOCATE (Z__Smooth_DBT(I__NumFoundChans)) + Z__Smooth_DBT(:) = 0.0_r_kind + + CALL CADS_Detect_Cloud_Smooth( I__NumFoundChans, I__Window_Width(JBAND), Z__DBT(I__INDEX(1:I__NumFoundChans)), & + Z__Smooth_DBT(1:I__NumFoundChans) ) + + +!3. Choice of cloud detection scenario + + CALL CADS_Detect_Cloud_Scenario( K__Sensor, JBAND, I__NumFoundChans, I__GradChkInterval(JBAND), I__Index(1:I__NumFoundChans), & + I__CHAN_HIGH, I__CHAN_LOW, I__Window_Chans, I__Imager_Flag, I__Scenario_Index, I__Start_Channel, Z__Smooth_DBT(1:I__NumFoundChans)) + + +!4. Identify the separation between clear/cloudy channels + + CALL CADS_Detect_Cloud_Separator( K__Sensor, JBAND, I__NumFoundChans, I__GradChkInterval(JBAND), I__Index(1:I__NumFoundChans), & + I__Cloud_Flag, I__FirstCloudyChannel, I__LastClearChannel, I__Scenario_Index, I__Start_Channel, Z__Smooth_DBT(1:I__NumFoundChans)) + + K__Cloud_Flag(IDCHAN(1:I__NumFoundChans)) = & + I__Cloud_Flag(1:I__NumFoundChans) + + ! Set cloud level for cross-band: + IF (I__FirstCloudyChannel == 0) THEN ! FOV is completely clear + Z__Cloud_Level = 1.e20_r_kind ! Large value + ELSE + Z__Cloud_Level = P__Chan_Level(IDCHAN(I__LastClearChannel)) + ENDIF + + ! Automatically do cross band cloud detection for all + ! interferometer channels (whether assigned a band or not) if + ! JBand == 1. This can be over-ridden for the other bands. + + IF (K__Sensor /= INST_ID_AIRS .AND. JBand == 1) & + WHERE(P__Chan_Level(:) < Z__Cloud_Level) K__Cloud_Flag(:) = 0 + + CrossBand : IF (LL__Do_CrossBand) THEN + ! Cross Band: + ! Loop through bands applying cloud detection to those that take their + ! cloud detection information from the current band JBAND. + DO JBand2 = 1, I__Num_Bands + IF (I__BandToUse(JBand2) == JBand) THEN + WHERE(P__Chan_Level(:) < Z__Cloud_Level .AND. & + I__BandNumber == JBand2 .AND. & + P__OBSBTs(:)>0.0_r_kind ) K__Cloud_Flag(:) = 0 + ENDIF + ENDDO + ENDIF CrossBand + +! Deallocate arrays + IF (ALLOCATED(Z__DBT)) DEALLOCATE (Z__DBT) + IF (ALLOCATED(Z__Smooth_DBT)) DEALLOCATE (Z__Smooth_DBT) + IF (ALLOCATED(Z__LEVEL)) DEALLOCATE (Z__LEVEL) + IF (ALLOCATED(I__Cloud_Flag)) DEALLOCATE (I__Cloud_Flag) + IF (ALLOCATED(I__INDEX)) DEALLOCATE (I__INDEX) + IF (ALLOCATED(IDCHAN)) DEALLOCATE (IDCHAN) + + ENDDO Band_Loop + +! Nullify pointers + NULLIFY(I__Band_Size, I__Bands, I__Window_Width, I__BandToUse) + +END SUBROUTINE CADS_Detect_Cloud + +SUBROUTINE CADS_Detect_Cloud_Imager( K__Sensor, K__Nchans, K__Chanid, K__Nclust, K__Cloud_Flag, P__Cl_Fraction, & + P__Cl_Mean, P__Ov_Stddev, P__FG_BT ) + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and MeteoFrance. +! +! Copyright 2020, EUMETSAT, All Rights Reserved. + +! *CADS_Detect_Cloud_Imager* +! R.Eresmaa ECMWF 12/02/13 + +! * PURPOSE * +! ----------- +! Provide additional information for the cloud detection by making use +! of collocated imager data, such as AVHRR collocated with IASI. + +! * INTERFACE * +! ------------- +! *CALL* * CADS_Detect_Cloud_Imager( )* (from CADS_Detect_Cloud) +! WHERE K__Sensor : Satellite sensor id +! K__Nchans : Number of channels received as input +! K__Chanid : Provided channel IDs +! K__Nclust : Highest possible number of clusters +! K__Cloud_Flag : Output cloud flag (0-7, 0=clear) +! P__Cl_Fraction : Fractional coverage of each cluster within FOV +! P__Cl_Mean : Cluster-mean brightness temperature (BT) on each +! channel +! P__Ov_Stddev : Overall BT standard deviation on each channel +! P__FG_BT : Forward-modelled BT on each channel + +! * METHOD * +! ---------- +! A preliminary indicator of presence of clouds in the sounder +! field-of-view (FOV) is derived using statistical radiance information +! within collocated clusters of imager pixels. + +! * MODIFICATIONS * +! ----------------- +! 03/12/13 R.Eresmaa 2.0 Original export version. +! 19/01/15 R.Eresmaa 2.1 Make array size specifications implicit. +! Verify that channels intended to be used +! are received as input. +! 05/02/19 R.Eresmaa 2.4 Explicit kind specifications. +! 16/04/20 R.Eresmaa 3.0 Rename and tidy up. + + use kinds, only: i_kind, r_kind + IMPLICIT NONE + +!* Global arrays + INTEGER(i_kind), INTENT(IN) :: K__Sensor ! Sensor id + INTEGER(i_kind), INTENT(IN) :: K__Nchans ! No. of channels + INTEGER(i_kind), INTENT(IN) :: K__Chanid(:) ! Channel IDs + INTEGER(i_kind), INTENT(IN) :: K__Nclust ! No. of clusters + INTEGER(i_kind), INTENT(OUT) :: K__Cloud_Flag ! Output cloud flag + REAL(r_kind), INTENT(IN) :: P__Cl_Fraction(:) ! Cluster fractions + REAL(r_kind), INTENT(IN) :: P__Cl_Mean(:,:) ! Cluster-mean BTs + REAL(r_kind), INTENT(IN) :: P__Ov_Stddev(:) ! Overall BT st.devs. + REAL(r_kind), INTENT(IN) :: P__FG_BT(:) ! First guess BT + +!* Local variables - Setup of the imager cloud detection + INTEGER(i_kind) :: I__Num_Imager_Chans ! No. of used channels + INTEGER(i_kind), POINTER :: I__Imager_Chans(:) ! List of used channels + REAL(r_kind), POINTER :: Z__Stddev_Threshold(:) ! Homogeneity thresholds + REAL(r_kind) :: Z__Coverage_Threshold ! Coverage threshold + REAL(r_kind) :: Z__FG_Departure_Threshold ! FG departure threshold + +!* Additional local variables + INTEGER(i_kind) :: I, J, IK, I_Temp_Flag, ICOUNT + INTEGER(i_kind) :: I__Chan_Index(K__Nchans) + REAL(r_kind) :: Z__Wsqdev, Z__Intercluster + REAL(r_kind),dimension(K__Nclust) :: Z__Sqdev + + + +!* 1.0 Initialize cloud flags as clear + + K__Cloud_Flag=0 + + IF (S__CADS_Setup_Cloud(K__Sensor) % L__Do_Imager_Cloud_Detection) THEN + + +!* 1.1 Setup + + I__Num_Imager_Chans = & + S__CADS_Setup_Cloud(K__Sensor) % N__Num_Imager_Chans + I__Imager_Chans => & + S__CADS_Setup_Cloud(K__Sensor) % N__Imager_Chans + Z__Stddev_Threshold => & + S__CADS_Setup_Cloud(K__Sensor) % R__Stddev_Threshold + Z__Coverage_Threshold = & + S__CADS_Setup_Cloud(K__Sensor) % R__Coverage_Threshold + Z__FG_Departure_Threshold = & + S__CADS_Setup_Cloud(K__Sensor) % R__FG_Departure_Threshold + + + +!* 1.2 Channel indexing + I__Chan_Index(:) = 0 + ICOUNT=0 + DO I=1,K__Nchans + IK=0 + DO J=1,I__Num_Imager_Chans + IF (K__Chanid(I)==I__Imager_Chans(J)) THEN + ICOUNT=ICOUNT+1 + IK=ICOUNT + EXIT + ENDIF + ENDDO + I__Chan_Index(I)=IK + ENDDO + + +!* 2.0 Compute squared first guess departures for each cluster + + DO J=1,K__Nclust + Z__Sqdev(J) = 0.0_r_kind + DO I=1,K__Nchans + IF (I__Chan_Index(I)==0) CYCLE + Z__Sqdev(J) = Z__Sqdev(J) + (P__Cl_Mean(I,J)-P__FG_BT(I))**2 + ENDDO + ENDDO + +!* 2.1 Homogeneity check: Do not diagnose presence of cloud if BT +! standard deviation falls below given threshold on at least one +! channel. + + I_Temp_Flag=1 + DO I=1,K__Nchans + IF (I__Chan_Index(I)==0) CYCLE + IF (P__Ov_Stddev(I)Z__Sqdev(J) .OR. Z__Intercluster>Z__Sqdev(IK)) THEN + K__Cloud_Flag=K__Cloud_Flag+2 + Exit Consistency_Check + ENDIF + ENDDO + ENDDO Consistency_Check + + +!* 2.3 First guess departure check: Do not diagnose presence of cloud +! if fraction-weighted first guess departure falls below given +! threshold. + + Z__Wsqdev = SUM(P__Cl_Fraction(:)*Z__Sqdev(:)) + IF (Z__Wsqdev>=Z__FG_Departure_Threshold) K__Cloud_Flag=K__Cloud_Flag+1 + + ENDIF ! L__Do_Imager_Cloud_Detection + +END SUBROUTINE CADS_Detect_Cloud_Imager + + +SUBROUTINE CADS_Detect_Cloud_Heapsort(N, A, K_Index) + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and MeteoFrance. +! +! Copyright 2020, EUMETSAT, All Rights Reserved. + +! * CADS_Detect_Cloud_Heapsort * +! A.Collard ECMWF 01/02/06 + +! * PURPOSE * +! ----------- +! Basic heapsort algorithm. + +! * INTERFACE * +! ------------- +! *CALL* * CADS_Detect_Cloud_Heapsort( )* (from CADS_Detect_Cloud) +! WHERE N : Length of input array +! A : Real input array +! K_Index : Output ranked array + +! * MODIFICATIONS * +! ----------------- +! 16/05/06 A.Collard 1.0 Original version. +! 05/02/19 R.Eresmaa 2.4 Explicit KIND specifications +! 16/04/20 R.Eresmaa 3.0 Rename as part of the big clean for CADS V3 + + + use kinds, only: i_kind, r_kind + IMPLICIT NONE + +! Subroutine arguments + INTEGER(i_kind), INTENT(IN) :: N + REAL(r_kind), INTENT(IN) :: A(:) + INTEGER(i_kind), INTENT(INOUT) :: K_Index(:) + + INTEGER(i_kind) :: I,J,RIGHT,LEFT,IDX + REAL(r_kind) :: TMP + +!------------------------------------------ + + IF (N <= 1) RETURN + LEFT = N/2+1 + RIGHT = N + + DO + IF (LEFT > 1) THEN + LEFT = LEFT - 1 + IDX = K_Index(LEFT) + ELSE + IDX = K_Index(RIGHT) + K_Index(RIGHT) = K_Index(1) + RIGHT = RIGHT - 1 + IF (RIGHT == 1) THEN + K_Index(1) = IDX + EXIT + ENDIF + ENDIF + TMP = A(IDX) + I = LEFT + J = 2*LEFT + DO WHILE (J <= RIGHT) + IF (J < RIGHT) THEN + IF (A(K_Index(J)) < A(K_Index(J+1))) J = J + 1 + ENDIF + IF (TMP < A(K_Index(J))) THEN + K_Index(I) = K_Index(J) + I = J + J = 2*J + ELSE + J = RIGHT + 1 + ENDIF + ENDDO + K_Index(I) = IDX + ENDDO + +END SUBROUTINE CADS_Detect_Cloud_Heapsort + +SUBROUTINE CADS_Detect_Cloud_Smooth(KV,KW,PV,PVA) + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and MeteoFrance. +! +! Copyright 2020, EUMETSAT, All Rights Reserved. + +! * CADS_Detect_Cloud_Smooth * - Boxcar-averaging in a REAL array +! * Phil Watts ECMWF 24/01/02 + +! * PURPOSE * +! ----------- +! Calculate the moving average (smoothing filter) of array +! No error checking supplied. + +! * INTERFACE * +! ------------- +! *CALL* * CADS_Detect_Cloud_Smooth( )* (from CADS_Detect_Cloud) +! WHERE KV : Number of elements in V +! KW : Window width for filter +! PV : Input array to be averaged +! PVA : Averaged array + +! * MODIFICATIONS * +! ----------------- +! 01/02/06 A.Collard 1.0 Original export version. +! 13/01/15 R.Eresmaa 2.1 Make array size specifications implicit. +! 05/02/19 R.Eresmaa 2.4 Explicit KIND specifications. +! 16/04/20 R.Eresmaa 3.0 Rename and tidy up. + + use kinds, only: i_kind, r_kind + IMPLICIT NONE + +!* 0.1 global variables + INTEGER(i_kind), INTENT(IN) :: KV ! length of V + INTEGER(i_kind), INTENT(IN) :: KW ! length of averaging window + REAL(r_kind), INTENT(IN) :: PV(:) ! original array + REAL(r_kind), INTENT(INOUT) :: PVA(:) ! averaged array + +!* 0.2 local variables + INTEGER(i_kind) :: INJ,J,I + + PVA(:)=0.0_r_kind + + DO I = 1,KV ! loop over array elements + INJ=0 + DO J=I-KW/2,I+KW/2,1 ! loop over window + IF (J > 0 .AND. J < (KV+1)) THEN ! if window element exists in + ! original array + INJ=INJ+1 + PVA(I)=PVA(I)+PV(J) ! add value + ENDIF + ENDDO + PVA(I)=PVA(I)/REAL(INJ) ! mean value + ENDDO + +END SUBROUTINE CADS_Detect_Cloud_Smooth + +SUBROUTINE CADS_Detect_Cloud_Scenario( K__Sensor, K__Band, K__NumChans, K__GradChkInterval, K__Index, K__Chan_High, & + K__Chan_Low, K__Chan_Windows, K__Imager_Flag, K__Scen_Index, K__Start_Channel, P__DBT) + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and MeteoFrance. +! +! Copyright 2020, EUMETSAT, All Rights Reserved. + +! * CADS_Detect_Cloud_Scenario * +! PHIL WATTS ECMWF 21/01/02 + +! * PURPOSE * +! ----------- +! Determine which of the three possible scenarios best describes +! the input data. +! Quick Exit - no cloud in the FOV +! Warm Start - warm cloud above relatively colder surface +! Cold Start - cold cloud above relatively warmer surface (most common) + +! * INTERFACE * +! ------------- +! * CALL* * CADS_Detect_Cloud_Scenario( )* (from CADS_Detect_Cloud) +! WHERE K__Sensor : Satellite sensor (AIRS/IASI/CrIS) +! K__Band : Band number +! K__NumChans : Number of channels in this band +! K__GradChkInterval : Gradient-checking interval +! K__Index : Ranking index for the input dBT signal +! K__Chan_High : High channel considered in initial minimum search +! K__Chan_Low : Low channel considered in initial minimum search +! K__Chan_Windows : Two channels defining longwave window +! K__Imager_Flag : Input flag from collocated imager data +! K__Scen_Index : Choice of cloud detection scenario (1, 2, or 3) +! K__Start_Channel : Channel index for the start of final search +! P__DBT : Input dBT signal + +! * MODIFICATIONS * +! ----------------- +! 03/02/06 A.Collard 1.0 Tidy up in preparation for IASI +! 03/05/06 A.Collard 1.0.1 Band size is now passed in (allows for +! missing channels). +! 04/05/06 A.Collard 1.0.2 The index of the first cloudy channel is now +! returned to allow cross-band cloud detection +! 16/02/07 A.Collard 1.0.3 Change to the padding to allow the bottom +! channel to be flagged as clear in a +! non-quickstart situation. +! 16/01/09 A.Collard 1.1 Gradient check on quick exit +! Start channel for cold start moved to highest +! channel where BT threshold exceeded +! 11/11/11 R.Eresmaa 1,2 Index of the lowest clear channel added to +! the output parameters. +! Change of the starting channel is no longer +! allowed in cases where gradient > -threshold. +! 04/12/13 R.Eresmaa 2.0 Allow quick exit only if collocated imager +! data supports hypothesis of a clear FOV +! 13/01/15 R.Eresmaa 2.1 Remove the need to create temporary array in +! the call to MOVINGA. +! the call to MOVINGA. +! 04/02/19 R.Eresmaa 2.4 Explicit KIND specifications. +! 16/04/20 R.Eresmaa 3.0 Divide the previous CF_Digital in two: +! Cloud_Scenario (here) and Cloud_Separator. + + + use kinds, only: i_kind, r_kind + IMPLICIT NONE + +!* 0.1 Global arrays + INTEGER(i_kind), INTENT(IN) :: K__SENSOR ! Sensor + INTEGER(i_kind), INTENT(IN) :: K__Band ! Band number + INTEGER(i_kind), INTENT(IN) :: K__NumChans ! Number of usable channels in band + INTEGER(i_kind), INTENT(IN) :: K__GradChkInterval ! Gradient-check interval + INTEGER(i_kind), INTENT(IN) :: K__INDEX(:) ! Ranking index for dBT + INTEGER(i_kind), INTENT(IN) :: K__Chan_High ! First channel clear of high stratospheric model errors + INTEGER(i_kind), INTENT(IN) :: K__Chan_Low ! Last channel clear of PBL humidity errors + INTEGER(i_kind), INTENT(IN) :: K__Chan_Windows(2) ! Two channels defining long-wave window bounds + INTEGER(i_kind), INTENT(IN) :: K__Imager_Flag ! Input imager cloud flag + INTEGER(i_kind), INTENT(OUT) :: K__Scen_Index ! Choice of scenario + INTEGER(i_kind), INTENT(OUT) :: K__Start_Channel ! Final starting channel + REAL(r_kind), INTENT(IN) :: P__DBT(:) ! Input ranked-smoothed dBT signal + +! Local variables + REAL(r_kind), ALLOCATABLE :: Z__DBT_w_Buffer(:) ! Smoothed-ranked DBT + INTEGER(i_kind) :: I__Buffer ! No. of buffer channels + INTEGER(i_kind) :: I__Start_Channel ! Primary starting channel for cloud search + INTEGER(i_kind) :: I__Start_Channel_Surf ! Secondary starting channel for cloud search + INTEGER(i_kind) :: I__Max_Channel ! Channel corresponding to maximum of the smoothed dBT + INTEGER(i_kind) :: JCH,JMIN(1),JMAX(1),I + + LOGICAL :: LLCOLD, LL__WINDOW_GRAD_CHECK, LL__StartChannelChanged + LOGICAL :: LL__Search_for_Cloud_Top + +! These carry the values in S__CADS_Setup_Cloud + REAL(r_kind) :: Z__BT_Threshold ! Solution contaminated threshold + REAL(r_kind) :: Z__Grad_Threshold ! Gradient threshold at which to stop filter procession + REAL(r_kind) :: Z__Window_Grad_Threshold ! Gradient threshold for window check + + +!============================================================================= + + + Z__BT_Threshold = & + S__CADS_Setup_Cloud(K__SENSOR) % R__BT_Threshold(K__Band) + Z__Grad_Threshold = & + S__CADS_Setup_Cloud(K__SENSOR) % R__Grad_Threshold(K__Band) + Z__Window_Grad_Threshold = & + S__CADS_Setup_Cloud(K__SENSOR) % R__Window_Grad_Threshold(K__Band) + + +!1. Include buffer channels at the start and end of the input smoothed +! departure array + + I__BUFFER = K__GradChkInterval + ALLOCATE(Z__DBT_w_Buffer(-I__Buffer+1:K__NumChans+1)) + + Z__DBT_w_Buffer(1:K__NumChans) = P__DBT(:) + Z__DBT_w_Buffer(-I__BUFFER+1:0) = Z__DBT_w_Buffer(1) + Z__DBT_w_Buffer(K__NumChans+1) = Z__DBT_w_Buffer(K__NumChans) + + +!2. Prepare for the cloud search + +! First define a set of key channels + + JMIN=MINLOC(Z__DBT_w_Buffer(K__Chan_High:K__NumChans)) + I__Start_Channel_Surf = K__Chan_High+JMIN(1)-1 + + JMIN=MINLOC(Z__DBT_w_Buffer(K__Chan_High:K__Chan_Low)) + I__Start_Channel = K__Chan_High+JMIN(1)-1 + +! Look for highest channel with DBT<-BT_Threshold and move I__Start_Channel +! there if higher than current I__Start_Channel: + JCH = I__Start_Channel + StartChanLoop : DO I=K__Chan_High,K__NumChans + IF (Z__DBT_w_Buffer(I) < -Z__BT_Threshold .OR. I == I__Start_Channel) THEN + JCH = I + Exit StartChanLoop + ENDIF + ENDDO StartChanLoop + I__Start_Channel = JCH + +! Do the same with I__Start_Channel_Surf + JCH = I__Start_Channel_Surf + StartChanLoop_Surf : DO I=K__Chan_High,K__NumChans + IF (Z__DBT_w_Buffer(I) < -Z__BT_Threshold .OR. I == I__Start_Channel_Surf) THEN + JCH = I + Exit StartChanLoop_Surf + ENDIF + ENDDO StartChanLoop_Surf + I__Start_Channel_Surf = JCH + +! Find the position of the equivalent maximum departure (for quick exit test) + JMAX=MAXLOC(Z__DBT_w_Buffer(K__Chan_High:K__NumChans)) + I__Max_Channel = K__Chan_High+JMAX(1)-1 + +! Long-wave window gradient check + LL__WINDOW_GRAD_CHECK=.TRUE. + IF (ALL(K__Chan_Windows > 0)) LL__WINDOW_GRAD_CHECK = & + (ABS(Z__DBT_w_Buffer(K__INDEX(K__Chan_Windows(1))) - & + Z__DBT_w_Buffer(K__INDEX(K__Chan_Windows(2)))) & + < Z__Window_Grad_Threshold) + +! Choose scenario to be followed + LL__Search_for_Cloud_Top=.TRUE. + IF (ABS(Z__DBT_w_Buffer(I__Start_Channel_Surf)) < Z__BT_Threshold .AND. & + ABS(Z__DBT_w_Buffer(I__Start_Channel)) < Z__BT_Threshold .AND. & + ABS(Z__DBT_w_Buffer(I__Max_Channel)) < Z__BT_Threshold .AND. & + ABS(Z__DBT_w_Buffer(K__NumChans)) < Z__BT_Threshold .AND. & + LL__WINDOW_GRAD_CHECK .AND. & + K__Imager_Flag==0 .AND. & + S__CADS_Setup_Cloud(K__SENSOR) % L__Do_Quick_Exit) THEN + !Quick exit + LL__Search_for_Cloud_Top=.FALSE. + ELSEIF (ABS(Z__DBT_w_Buffer(I__Start_Channel)) < Z__BT_Threshold .AND. & + Z__DBT_w_Buffer(K__NumChans) > Z__BT_Threshold ) THEN + !Warm cloud start at next-to-bottom channel (allowing one channel for + !gradient calculations). + LLCOLD = .FALSE. + I__Start_Channel = K__NumChans-1 + ELSEIF (Z__DBT_w_Buffer(I__Start_Channel) < -Z__BT_Threshold ) THEN + LLCOLD = .TRUE. + ELSEIF (Z__DBT_w_Buffer(I__Start_Channel) > Z__BT_Threshold ) THEN + LLCOLD = .FALSE. + ELSE + LLCOLD = .TRUE. + ENDIF + + IF (LL__Search_for_Cloud_Top) THEN ! Either cold or warm start + ! (but not quick exit) + + JCH=I__Start_Channel + +! Re-evaluate the choice of scenario: +! If the primary starting channel appears clear, and the secondary +! starting channel is lower, start from the latter. In that case +! re-evaluate whether cold or warm start is more appropriate. + IF (I__Start_Channel /= I__Start_Channel_Surf) THEN + + LL__StartChannelChanged = .FALSE. + IF (LLCOLD .AND. ( (Z__DBT_w_Buffer(JCH-1)-Z__DBT_w_Buffer(JCH+1)) < & + Z__Grad_Threshold .AND. & + Z__DBT_w_Buffer(JCH-K__GradChkInterval)-Z__DBT_w_Buffer(JCH+1) < & + Z__Grad_Threshold .AND. & + ABS(Z__DBT_w_Buffer(JCH)) < Z__BT_Threshold)) THEN + I__Start_Channel = I__Start_Channel_Surf + LL__StartChannelChanged = .TRUE. + ENDIF + + IF (LL__StartChannelChanged) THEN + + IF (ABS(Z__DBT_w_Buffer(I__Start_Channel)) < Z__BT_Threshold .AND. & + Z__DBT_w_Buffer(K__NumChans) > Z__BT_Threshold ) THEN + !Warm cloud start at next-to-bottom channel (allowing one channel for + !gradient calculations). + LLCOLD = .FALSE. + I__Start_Channel = K__NumChans-1 + ELSEIF (Z__DBT_w_Buffer(I__Start_Channel) < -Z__BT_Threshold ) THEN + LLCOLD = .TRUE. + ELSEIF (Z__DBT_w_Buffer(I__Start_Channel) > Z__BT_Threshold ) THEN + LLCOLD = .FALSE. + ELSE + LLCOLD = .TRUE. + ENDIF + JCH = I__Start_Channel + + ENDIF + ENDIF + + IF (LLCOLD) THEN + K__Scen_Index=3 + ELSE + K__Scen_Index=2 + ENDIF + K__Start_Channel = JCH + + ELSE + + K__Scen_Index=1 + K__Start_Channel=0 + + ENDIF ! Search for cloud top + + IF (ALLOCATED(Z__DBT_w_Buffer)) DEALLOCATE(Z__DBT_w_Buffer) + +END SUBROUTINE CADS_Detect_Cloud_Scenario + +SUBROUTINE CADS_Detect_Cloud_Separator( K__Sensor, K__Band, K__NumChans, K__GradChkInterval, K__Index, K__Cloud_Flag, & + K__Cloud_Level, K__Clear_Level, K__Scen_Index, K__Start_Channel, P__DBT) + +! This software was developed within the context of the EUMETSAT +! Satellite Application Facility on Numerical Weather Prediction +! (NWP SAF), under the Cooperation Agreement dated 7 December 2016, +! between EUMETSAT and the Met Office, UK, by one or more partners +! within the NWP SAF. The partners in the NWP SAF are the Met +! Office, ECMWF, DWD and MeteoFrance. +! +! Copyright 2020, EUMETSAT, All Rights Reserved. + +! * CADS_Detect_Cloud_Separator * +! PHIL WATTS ECMWF 21/01/02 + +! * PURPOSE * +! ----------- +! Along the vertically-ranked and smoothed array of departures, find +! the separating point at which all cloud-affected channels are on +! one side and all clear channels are on the other side. + +! * INTERFACE * +! ------------ +! * CALL* * CADS_Detect_Cloud_Separator( )* (from CADS_Detect_Cloud) +! WHERE K__Sensor : Satellite sensor (AIRS/IASI/CrIS) +! K__Band : Band number +! K__NumChans : Number of channels in this band +! K__GradChkInterval : Gradient-checking interval +! K__Index : Ranking index for the input dBT signal +! K__Cloud_Flag : Cloud flag by channel; 0=clear, 1=cloudy +! K__Cloud_Level : Index of the highest cloud-contaminated channel +! K__Clear_Level : Index of the lowest clear channel +! K__Scen_Index : Choice of cloud detection scenario (1, 2, or 3) +! K__Start_Channel : Starting channel for the cloud search +! P__DBT : Input dBT signal + +! MODIFICATIONS +! 03/02/06 A.Collard 1.0 Tidy up in preparation for IASI +! 03/05/06 A.Collard 1.0.1 Band size is now passed in (allows for +! missing channels). +! 04/05/06 A.Collard 1.0.2 The index of the first cloudy channel is now +! returned to allow cross-band cloud detection +! 16/02/07 A.Collard 1.0.3 Change to the padding to allow the bottom +! channel to be flagged as clear in a +! non-quickstart situation. +! 16/01/09 A.Collard 1.1 Gradient check on quick exit +! Start channel for cold start moved to highest +! channel where BT threshold exceeded +! 11/11/11 R.Eresmaa 1,2 Index of the lowest clear channel added to +! the output parameters. +! Change of the starting channel is no longer +! allowed in cases where gradient > -threshold. +! 04/12/13 R.Eresmaa 2.0 Allow quick exit only if collocated imager +! data supports hypothesis of a clear FOV +! 13/01/15 R.Eresmaa 2.1 Remove the need to create temporary array in +! the call to MOVINGA. +! 04/02/19 R.Eresmaa 2.4 Explicit KIND specifications. +! 16/04/20 R.Eresmaa 3.0 Divide the previous CF_Digital in two: +! Cloud_Scenario and Cloud_Separator (here). + + use kinds, only: i_kind, r_kind + IMPLICIT NONE + +!* 0.1 Global arrays + INTEGER(i_kind), INTENT(IN ) :: K__SENSOR ! Sensor + INTEGER(i_kind), INTENT(IN ) :: K__Band ! Band number + INTEGER(i_kind), INTENT(IN ) :: K__NumChans ! Number of usable channels in band + INTEGER(i_kind), INTENT(IN ) :: K__GradChkInterval ! Gradient-check interval + INTEGER(i_kind), INTENT(IN ) :: K__INDEX(:) ! Ranking index for dBT + INTEGER(i_kind), INTENT(INOUT) :: K__Cloud_Flag(:) ! Cloud flags + INTEGER(i_kind), INTENT( OUT) :: K__Cloud_Level ! Index of highest cloudy channel + INTEGER(i_kind), INTENT( OUT) :: K__Clear_Level ! Index of lowest clear channel + INTEGER(i_kind), INTENT(IN ) :: K__Scen_Index ! Choice of scenario + INTEGER(i_kind), INTENT(IN ) :: K__Start_Channel ! Choice of scenario + REAL(r_kind), INTENT(IN ) :: P__DBT(:) ! Input ranked dBT signal + + +! Local variables + REAL(r_kind), ALLOCATABLE :: Z__DBT_w_Buffer(:) ! Smoothed-ranked DBT + INTEGER(i_kind) :: I__Buffer ! No. of buffer channels + INTEGER(i_kind) :: JCH + +! These carry the values in S__CADS_Setup_Cloud + REAL(r_kind) :: Z__BT_Threshold ! Solution contaminated threshold + REAL(r_kind) :: Z__Grad_Threshold ! Gradient threshold at which to stop + ! filter procession + +!============================================================================= + + + Z__BT_Threshold = & + S__CADS_Setup_Cloud(K__SENSOR) % R__BT_Threshold(K__Band) + Z__Grad_Threshold = & + S__CADS_Setup_Cloud(K__SENSOR) % R__Grad_Threshold(K__Band) + + K__Cloud_Flag(:)=1 + +!1. Include buffer channels at the start and end of the input smoothed +! departure array + + I__BUFFER = K__GradChkInterval + ALLOCATE(Z__DBT_w_Buffer(-I__Buffer+1:K__NumChans+1)) + + Z__DBT_w_Buffer(1:K__NumChans) = P__DBT(:) + Z__DBT_w_Buffer(-I__BUFFER+1:0) = Z__DBT_w_Buffer(1) + Z__DBT_w_Buffer(K__NumChans+1) = Z__DBT_w_Buffer(K__NumChans) + + +!2. Search for the lowest non-contaminated channel + + JCH = K__Start_Channel + + SELECT CASE (K__Scen_Index) + + CASE (1) ! Quick Exit + K__Cloud_Level = 0 + + CASE (2) ! Warm Start +! In the case of Warm Start, progress towards higher channels whilst +! -ve difference is decreasing + DO WHILE ( ((Z__DBT_w_Buffer(JCH-1)-Z__DBT_w_Buffer(JCH+1)) < & + -1.0_r_kind * Z__Grad_Threshold .OR. & + (Z__DBT_w_Buffer(JCH-K__GradChkInterval)-Z__DBT_w_Buffer(JCH+1)) < & + -1.0_r_kind * Z__Grad_Threshold .OR. & + ABS(Z__DBT_w_Buffer(JCH)) > Z__BT_Threshold) .AND. JCH > 1 ) + JCH = JCH-1 + ENDDO + K__Cloud_Level = JCH + + CASE (3) ! Cold Start +! In the case of Cold Start, progress towards higher channels whilst +! -ve difference is decreasing + DO WHILE (( (Z__DBT_w_Buffer(JCH-1)-Z__DBT_w_Buffer(JCH+1)) > & + Z__Grad_Threshold .OR. & + (Z__DBT_w_Buffer(JCH-K__GradChkInterval)-Z__DBT_w_Buffer(JCH+1)) > & + Z__Grad_Threshold .OR. & + ABS(Z__DBT_w_Buffer(JCH)) > Z__BT_Threshold) .AND. JCH > 1 ) + JCH = JCH-1 + ENDDO + K__Cloud_Level = JCH + + CASE DEFAULT + RETURN + + END SELECT + +!3. Output channel indices for the highest cloud and lowest clear levels + IF (K__Cloud_Level > 1) THEN + K__Cloud_Flag(K__INDEX(1:K__Cloud_Level-1))=0 + K__Clear_Level=K__INDEX(K__Cloud_Level-1) + K__Cloud_Level=K__INDEX(K__Cloud_Level) + ELSEIF (K__Cloud_Level>0) THEN + K__Clear_Level=K__INDEX(K__Cloud_Level) + K__Cloud_Level=K__INDEX(K__Cloud_Level) + ELSE + K__Cloud_Flag(:)=0 + ENDIF + + IF (ALLOCATED(Z__DBT_w_Buffer)) DEALLOCATE(Z__DBT_w_Buffer) + +END SUBROUTINE CADS_Detect_Cloud_Separator + +subroutine cads_imager_calc(obstype,isis,nobs,nreal,nchanl,nsig,data_s,init_pass,mype, & + imager_cluster_fraction,imager_cluster_bt,imager_chan_stdev, imager_model_bt) + +!$$$ subprogram documentation block +! +! subprogram: cads_imager_calc compute model equivalent to the imager channels used by CADS +! prgmmr: Jung +! +! abstract: accumulate the data necessary to derive the model equivalent brightness temperatures +! used by the cloud and aerosol detection software for the imager cloud tests. +! +! program history log: +! +! +! +! subroutines included: +! +! +! input argument list: +! +! obstype - type of tb observation +! isis - sensor/instrument/satellite id +! nobs - number of observations +! nreal - number of pieces of info (location, time, etc) per obs +! nchanl - number of channels per obs +! nsig - number of model layers +! data_s - array containing input data information for a specific sensor +! init_pass - state of "setup" processing +! mype - mpi task id +! +! output argument list: + +! imager_cluster_fraction - CADS cluster fraction ( dimension 7) +! imager_cluster_bt - avreage brightness temperature of a cluster +! imager_chan stdev - brightness temperature standard deviation of the cluster +! imager_model_bt - model derived brightness temperature +! +! +!$$$ end documentation block + + use kinds, only: i_kind, r_kind + use constants, only: zero + use radiance_mod, only: rad_obs_type + use radinfo, only: jpch_rad, nusis, crtm_coeffs_path, nsigradjac + use crtm_interface, only: init_crtm, call_crtm, destroy_crtm, itime + use obsmod, only: dval_use + use gsi_nstcouplermod, only: nstinfo + + implicit none + + logical, intent(in) :: init_pass + character(len=10), intent(in) :: obstype + character(len=20), intent(in) :: isis + integer(i_kind), intent(in) :: nobs, nreal, nchanl, nsig + integer(i_kind), intent(in) :: mype + real(r_kind),dimension(nreal+nchanl,nobs),intent(in) :: data_s + real(r_kind),dimension(7,nobs), intent(out) :: imager_cluster_fraction + real(r_kind),dimension(2,7,nobs), intent(out) :: imager_cluster_bt + real(r_kind),dimension(2,nobs), intent(out) :: imager_chan_stdev, imager_model_bt + +! local variables + integer(i_kind) :: jc, i, n + integer(i_kind) :: itmp1_cads, itmp2_cads, nchanl_cads, maxinfo, dval_info, cads_info, error_status + integer(i_kind),allocatable,dimension(:) :: ich_cads + logical :: imager_spccoeff, imager_taucoeff + real(r_kind) :: dtime, clw_guess, ciw_guess, rain_guess, snow_guess + real(r_kind) :: trop5, tzbgr, dtsavg, sfc_speed + real(r_kind),dimension(nsig) :: qvp, tvp, qs, prsltmp + real(r_kind),dimension(nsig+1) :: prsitmp + real(r_kind),allocatable,dimension(:) :: tsim_cads, emissivity_cads, chan_level_cads + real(r_kind),allocatable,dimension(:) :: ts_cads, emissivity_k_cads,data_s_cads + real(r_kind),allocatable,dimension(:,:) :: ptau5_cads, temp_cads, wmix_cads, jacobian_cads + character(len=80) :: spc_filename, tau_filename + character(len=20) :: isis_cads + character(len=10) :: obstype_cads + + type(rad_obs_type) :: radmod + + cads_info = 23 + dval_info = 0 + if (dval_use) dval_info = 2 + + itmp1_cads = len(trim(obstype)) + itmp2_cads = len(trim(isis)) + + if ( obstype == 'iasi' ) then + isis_cads = 'avhrr3'//isis(itmp1_cads+1:itmp2_cads) + obstype_cads = 'avhrr' +! nchanl_cads = 3 !channels 3 - 5 + elseif ( obstype == 'cris' .or. obstype == 'cris-fsr' ) then +! isis_cads = 'viirs-m'//isis(itmp1+1:itmp2) When naming convention becomes standarized with CrIS + if ( isis == 'cris-fsr_npp' .or. isis == 'cris_npp' ) then + isis_cads = 'viirs-m_npp' + elseif ( isis == 'cris-fsr_n20' ) then + isis_cads = 'viirs-m_n20' + spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin' + inquire(file=trim(spc_filename), exist=imager_spccoeff) + if ( .not. imager_spccoeff ) isis_cads = 'viirs-m_j1' + elseif ( isis == 'cris-fsr_n21' ) then + isis_cads = 'viirs-m_n21' + spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin' + inquire(file=trim(spc_filename), exist=imager_spccoeff) + if ( .not. imager_spccoeff ) isis_cads = 'viirs-m_j2' + endif + obstype_cads = 'viirs-m' +! nchanl_cads = 5 ! channels 12 - 16 + endif + + spc_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.SpcCoeff.bin' + inquire(file=trim(spc_filename), exist=imager_spccoeff) + tau_filename = trim(crtm_coeffs_path)//trim(isis_cads)//'.TauCoeff.bin' + inquire(file=trim(tau_filename), exist=imager_taucoeff) + +! IF the RTM files exist allocate and setup various arrays for the RTM + if ( imager_spccoeff .and. imager_taucoeff) then + nchanl_cads = 0 + do i=1,jpch_rad + if (trim(isis_cads) == nusis(i)) then + nchanl_cads = nchanl_cads +1 + endif + end do + + allocate( ich_cads(nchanl_cads) ) + jc = 0 + do i=1,jpch_rad + if (trim(isis_cads) == nusis(i)) then + jc = jc +1 + ich_cads(jc) = i + endif + end do + + call init_crtm(init_pass,-99,mype,nchanl_cads,nreal,isis_cads,obstype_cads,radmod) + +! Initialize variables needed for the infrared cloud and aerosol detection software + allocate(data_s_cads(nreal+nchanl_cads),tsim_cads(nchanl_cads),emissivity_cads(nchanl_cads), & + chan_level_cads(nchanl_cads),ptau5_cads(nsig,nchanl_cads),ts_cads(nchanl_cads),emissivity_k_cads(nchanl_cads), & + temp_cads(nsig,nchanl_cads),wmix_cads(nsig,nchanl_cads), jacobian_cads(nsigradjac,nchanl_cads)) + + do n = 1,nobs ! loop to derive imager BTs for CADS +! Extract analysis relative observation time. + dtime = data_s(itime,n) + maxinfo = nreal - cads_info - dval_info - nstinfo + if ( sum(data_s(maxinfo+1:maxinfo+7,n)) > 0.90_r_kind ) then ! imager cluster information exists for this profile + data_s_cads = data_s(1:nreal+nchanl_cads,n) + call call_crtm(obstype_cads,dtime,data_s_cads,nchanl_cads,nreal,ich_cads, & + tvp,qvp,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & + trop5,tzbgr,dtsavg,sfc_speed,tsim_cads,emissivity_cads,chan_level_cads, & + ptau5_cads,ts_cads,emissivity_k_cads,temp_cads,wmix_cads,jacobian_cads,error_status) + +! Transfer imager data to arrays for qc_irsnd + imager_cluster_fraction(1:7,n) = data_s(maxinfo+1:maxinfo+7,n) + imager_cluster_bt(1,1:7,n) = data_s(maxinfo+8:maxinfo+14,n) + imager_cluster_bt(2,1:7,n) = data_s(maxinfo+15:maxinfo+21,n) + imager_chan_stdev(1:2,n) = data_s(maxinfo+22:maxinfo+23,n) + imager_model_bt(1:2,n) = tsim_cads(nchanl_cads-1:nchanl_cads) + endif ! imager information exists + end do ! End loop to derive imager BTs + + call destroy_crtm + deallocate(data_s_cads,tsim_cads,emissivity_cads, ich_cads,chan_level_cads,ptau5_cads,& + ts_cads,emissivity_k_cads, temp_cads,wmix_cads, jacobian_cads) + endif ! RTM files exist + + end subroutine cads_imager_calc + +end module cads diff --git a/src/gsi/constants.f90 b/src/gsi/constants.f90 index c1c5f3f879..291a18ec97 100644 --- a/src/gsi/constants.f90 +++ b/src/gsi/constants.f90 @@ -31,7 +31,6 @@ module constants ! 2016-02-15 Johnson, Y. Wang, X. Wang - define additional constant values for ! radar DA, POC: xuguang.wang@ou.edu ! 2019-09-25 X.Su - put stndrd_atmos_ps constant values -! 2022-03-01 X.Lu & X.Wang - increased max_varname_length for HAFS dual ens. POC: xuguang.wang@ou.edu ! ! Subroutines Included: ! sub init_constants_derived - compute derived constants diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90 index 2305c84340..4bb1191001 100644 --- a/src/gsi/crtm_interface.f90 +++ b/src/gsi/crtm_interface.f90 @@ -977,7 +977,7 @@ end subroutine destroy_crtm subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & h,q,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsl,prsi, & trop5,tzbgr,dtsavg,sfc_speed,& - tsim,emissivity,ptau5,ts, & + tsim,emissivity,chan_level,ptau5,ts, & emissivity_k,temp,wmix,jacobian,error_status,tsim_clr,tcc, & tcwv,hwp_ratio,stability,layer_od,jacobian_aero) !$$$ subprogram documentation block @@ -1097,6 +1097,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind) ,intent( out) :: sfc_speed,dtsavg real(r_kind),dimension(nchanl+nreal) ,intent(in ) :: data_s real(r_kind),dimension(nchanl) ,intent( out) :: tsim,emissivity,ts,emissivity_k + real(r_kind),dimension(nchanl) ,intent( out) :: chan_level character(10) ,intent(in ) :: obstype integer(i_kind) ,intent( out) :: error_status real(r_kind),dimension(nsig,nchanl) ,intent( out) :: temp,ptau5,wmix @@ -1150,6 +1151,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & real(r_kind):: sno00,sno01,sno10,sno11,secant_term real(r_kind):: hwp_total,theta_700,theta_sfc,hs real(r_kind):: dlon,dlat,dxx,dyy,yy,zz,garea + real(r_kind):: radiance, radiance_overcast, radiance_ratio real(r_kind),dimension(0:3):: wgtavg real(r_kind),dimension(nsig,nchanl):: omix real(r_kind),dimension(nsig,nchanl,n_aerosols_jac):: jaero @@ -2217,8 +2219,10 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & end do end if + chan_level = zero + !$omp parallel do schedule(dynamic,1) private(i) & -!$omp private(total_od,k,kk,m,term,ii,cwj) +!$omp private(total_od,k,kk,m,term,ii,cwj,radiance,radiance_overcast,radiance_ratio) do i=1,nchanl ! Zero jacobian and transmittance arrays do k=1,nsig @@ -2228,6 +2232,16 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, & wmix(k,i)=zero end do + radiance=rtsolution(i,1)%radiance + do k=msig, 1, -1 + radiance_overcast = rtsolution(i,1)%upwelling_overcast_radiance(k) + radiance_ratio = abs(radiance_overcast/radiance) + if (radiance_ratio < 0.99_r_kind) then + chan_level(i) = atmosphere(1)%pressure(k) / r10 + exit + endif + enddo + ! Simulated brightness temperatures tsim(i)=rtsolution(i,1)%brightness_temperature diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake index 5a7d29c208..ce74d91c63 100644 --- a/src/gsi/gsi_files.cmake +++ b/src/gsi/gsi_files.cmake @@ -101,6 +101,7 @@ bkgvar_rewgt.f90 blacklist.f90 blendmod.f90 buddycheck_mod.f90 +cads.f90 calc_fov_conical.f90 calc_fov_crosstrk.f90 calctends.f90 diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90 index 43fb40e8ad..8158f35e11 100644 --- a/src/gsi/gsi_rfv3io_mod.f90 +++ b/src/gsi/gsi_rfv3io_mod.f90 @@ -559,7 +559,7 @@ subroutine gsi_rfv3io_get_ens_grid_specs(grid_spec,ierr) integer(i_kind) :: nio,nylen integer(i_kind),allocatable :: gfile_loc_layout(:) character(len=180) :: filename_layout - integer(i_kind) imiddle,jmiddle + integer(i_kind) imiddle,jmiddle,grid_ens_type_fv3_regional iret=nf90_open(trim(grid_spec),nf90_nowrite,gfile_grid_spec) @@ -635,33 +635,19 @@ subroutine gsi_rfv3io_get_ens_grid_specs(grid_spec,ierr) if(mype==0)write(6,*),'ny_layout_bens=',ny_layout_bens if(mype==0)write(6,*),'ny_layout_eens=',ny_layout_eens -! -! need to decide the grid orientation of the FV regional model -! -! grid_type_fv3_regional = 0 : decide grid orientation based on -! grid_lat/grid_lon -! 1 : input is E-W N-S grid -! 2 : input is W-E S-N grid -! - if(grid_type_fv3_regional == 0) then - imiddle=nxens/2 - jmiddle=nyens/2 - if( (grid_latt(imiddle,1) < grid_latt(imiddle,nyens)) .and. & - (grid_lont(1,jmiddle) < grid_lont(nxens,jmiddle)) ) then - grid_type_fv3_regional = 2 - else - grid_type_fv3_regional = 1 - endif + imiddle=nxens/2 + jmiddle=nyens/2 + if( (grid_latt(imiddle,1) < grid_latt(imiddle,nyens)) .and. & + (grid_lont(1,jmiddle) < grid_lont(nxens,jmiddle)) ) then + grid_ens_type_fv3_regional = 2 + else + grid_ens_type_fv3_regional = 1 endif ! check the grid type - if( grid_type_fv3_regional == 1 ) then - if(mype==0) write(6,*) 'FV3 regional input grid is E-W N-S grid' - grid_reverse_flag=.true. ! grid is revered comparing to usual map view - else if(grid_type_fv3_regional == 2) then - if(mype==0) write(6,*) 'FV3 regional input grid is W-E S-N grid' - grid_reverse_flag=.false. ! grid orientated just like we see on map view + if( grid_type_fv3_regional == grid_ens_type_fv3_regional ) then + if(mype==0) write(6,*) 'Ensemble has the same orientation as the control, Cool!' else - write(6,*) 'Error: FV3 regional input grid is unknown grid' + write(6,*) 'Warning! Ensemble has a different orientation as the control. This case needs further tests, Abort!' call stop2(678) endif ! diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index c24c485ce1..45d88887a3 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -93,8 +93,14 @@ module gsimod erradar_inflate,tdrerr_inflate,use_poq7,qc_satwnds,& init_qcvars,vadfile,noiqc,c_varqc,gps_jacqc,qc_noirjaco3,qc_noirjaco3_pole,& buddycheck_t,buddydiag_save,njqc,vqc,nvqc,hub_norm,vadwnd_l2rw_qc, & - pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check + pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check, & + cris_cads, iasi_cads, airs_cads use qcmod, only: troflg,lat_c,nrand + use cads, only: M__Sensor,N__Num_Bands,N__GradChkInterval,N__Band_Size,N__Bands,N__Window_Width, & + N__Window_Bounds,R__BT_Threshold,R__Grad_Threshold,R__Window_Grad_Threshold, L__Do_Quick_Exit, & + L__Do_CrossBand, N__BandToUse,L__Do_Imager_Cloud_Detection, N__Num_Imager_Chans, & + N__Num_Imager_Clusters,N__Imager_Chans,R__Stddev_Threshold,R__Coverage_Threshold, & + R__FG_Departure_Threshold, CADS_Setup_Cloud use pcpinfo, only: npredp,diag_pcp,dtphys,deltim,init_pcp use jfunc, only: iout_iter,iguess,miter,factqmin,factqmax,superfact,limitqobs, & factql,factqi,factqr,factqs,factqg, & @@ -507,6 +513,9 @@ module gsimod ! 2. fv3_regional = .true. ! 3. fv3_cmaq_regional = .true. ! 4. berror_fv3_cmaq_regional = .true. +! 09-02-2022 Jung Added namelist entries to call a new IR cloud detection routine +! the original cloud detection routine is the default. To use the new +! cloud detection routine, set the flags to .true. ! 09-15-2022 yokota - add scale/variable/time-dependent localization ! 2023-07-30 Zhao - added namelist options for analysis of significant wave height ! (aka howv in GSI code): corp_howv, hwllp_howv @@ -1051,6 +1060,13 @@ module gsimod ! wind observations. ! vad_near_analtime - assimilate newvadwnd obs around analysis time only +! +! Flags to use the new IR cloud detection routine. Flag must be set to true to use the new routine. The default +! (no flag or .false.) will use the default. +! airs_cads: use the clod and aerosool detection software for the AIRS instrument +! cris_cads: use the cloud and aerosol detection software for CrIS instruments +! iasi_cads: use the cloud and aerosol detection software for IASI instruments +! namelist/obsqc/dfact,dfact1,erradar_inflate,tdrerr_inflate,oberrflg,& vadfile,noiqc,c_varqc,blacklst,use_poq7,hilbert_curve,tcp_refps,tcp_width,& @@ -1061,7 +1077,7 @@ module gsimod q_doe_a_136,q_doe_a_137,q_doe_b_136,q_doe_b_137, & t_doe_a_136,t_doe_a_137,t_doe_b_136,t_doe_b_137, & uv_doe_a_236,uv_doe_a_237,uv_doe_a_213,uv_doe_b_236,uv_doe_b_237,uv_doe_b_213, & - vad_near_analtime + vad_near_analtime,airs_cads,cris_cads,iasi_cads ! OBS_INPUT (controls input data): ! dmesh(max(dthin))- thinning mesh for each group @@ -1663,6 +1679,40 @@ module gsimod ! fac_tsl - index to apply thermal skin layer or not: 0 = no; 1 = yes. namelist/nst/nst_gsi,nstinfo,zsea1,zsea2,fac_dtl,fac_tsl +! Initialize the Cloud and Aerosol Detection Software (CADS) +! +! M__Sensor Unique ID for sensor +! N__Num_Bands Number of channel bands +! N__Band_Size(:) Number of channels in each band +! N__Bands(:,:) Channel lists +! N__Window_Width(:) Smoothing filter window widths per band +! N__Window_Bounds(:,:) Channels in the spectral window gradient check +! N__GradChkInterval(:) Window width used in gradient calculation +! R__BT_Threshold(:) BT threshold for cloud contamination +! R__Grad_Threshold(:) Gradient threshold for cloud contamination +! R__Window_Grad_Threshold(:) Threshold for window gradient check in QE +! L__Do_Quick_Exit On/off switch for the Quick Exit scenario +! L__Do_CrossBand On/off switch for the cross-band method +! N__BandToUse(:) Band number assignment for each channel +! L__Do_Imager_Cloud_Detection On/off switch for the imager cloud detection +! N__Num_Imager_Chans No. of imager channels +! N__Num_Imager_Clusters No. of clusters to be expected +! N__Imager_Chans(:) List of imager channels +! R__Stddev_Threshold(:) St. Dev. threshold, one for each imager channel +! R__Coverage_Threshold Threshold for fractional coverage of a cluster +! R__FG_Departure_Threshold Threshold for imager FG departure + + NAMELIST / Cloud_Detect_Coeffs / M__Sensor, N__Num_Bands, & + N__Band_Size, N__Bands, N__Window_Width, N__Window_Bounds, & + N__GradChkInterval, R__BT_Threshold, R__Grad_Threshold, & + R__Window_Grad_Threshold, L__Do_Quick_Exit, & + L__Do_CrossBand, N__BandToUse, & + L__Do_Imager_Cloud_Detection, N__Num_Imager_Chans, & + N__Num_Imager_Clusters, N__Imager_Chans, & + R__Stddev_Threshold, R__Coverage_Threshold, & + R__FG_Departure_Threshold + + !EOC !--------------------------------------------------------------------------- @@ -1749,6 +1799,7 @@ subroutine gsimain_initialize call set_fgrid2agrid call gsi_nstcoupler_init_nml call init_radaruse_directDA + call CADS_Setup_Cloud if(mype==0) write(6,*)' at 0 in gsimod, use_gfs_stratosphere,nems_nmmb_regional = ', & use_gfs_stratosphere,nems_nmmb_regional diff --git a/src/gsi/hybrid_ensemble_isotropic.F90 b/src/gsi/hybrid_ensemble_isotropic.F90 index 63af28cc5e..ef6b53119c 100644 --- a/src/gsi/hybrid_ensemble_isotropic.F90 +++ b/src/gsi/hybrid_ensemble_isotropic.F90 @@ -1393,10 +1393,6 @@ subroutine load_ensemble endif ! regional_ensemble_option = 5: ensembles are fv3 regional. - if (l_both_fv3sar_gfs_ens) then ! first read in gfs ensembles for regional - call get_gefs_for_regional -!clthink do we need mpi_bar here? - endif call fv3_regional_enspert%get_fv3_regional_ensperts(en_perts,nelen,ps_bar) @@ -4042,7 +4038,7 @@ subroutine hybens_grid_setup logical,allocatable::vector(:) real(r_kind) eps,r_e real(r_kind) rlon_a(nlon),rlat_a(nlat),rlon_e(nlon),rlat_e(nlat) - character(:),allocatable:: fv3_spec_grid_filename + character(:),allocatable:: fv3_ens_spec_grid_filename integer :: ierr nord_e2a=4 ! soon, move this to hybrid_ensemble_parameters @@ -4130,9 +4126,9 @@ subroutine hybens_grid_setup else if(dual_res) then call get_region_dx_dy_ens(region_dx_ens,region_dy_ens) - if(regional_ensemble_option) then - fv3_spec_grid_filename="fv3_ens_grid_spec" - call gsi_rfv3io_get_ens_grid_specs(fv3_spec_grid_filename,ierr) + if(regional_ensemble_option == 5) then + fv3_ens_spec_grid_filename="fv3_ens_grid_spec" + call gsi_rfv3io_get_ens_grid_specs(fv3_ens_spec_grid_filename,ierr) endif else region_dx_ens=region_dx diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 7146ceff3e..9804965573 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -115,6 +115,9 @@ module qcmod ! def vadfile - local name of bufr file containing vad winds (used by read_radar) ! def use_poq7 - if true, accept sbuv/2 obs with profile ozone quality flag 7 ! def cao_check - if true, turn on cold-air-outbreak screening +! def airs_cads - if true, use the cloud and aerosol detection routine for Aqua/AIRS instrument +! def cris_cads - if true, use the cloud and aerosol detection routine for CrIS instruments +! def iasi_cads - if true, use the cloud and aerosol detection routine for IASI instruments ! ! following used for nonlinear qc: ! @@ -152,7 +155,7 @@ module qcmod use constants, only: r0_01,r0_02,r0_03,r0_04,r0_05,r10,r60,r100,h300,r400,r1000,r2000,r2400,r4000 use constants, only: deg2rad,rad2deg,t0c,one_tenth,rearth_equator use obsmod, only: rmiss_single - use radinfo, only: iuse_rad,passive_bc + use radinfo, only: iuse_rad,passive_bc,nuchan use radinfo, only: tzr_qc use radiance_mod, only: rad_obs_type implicit none @@ -183,6 +186,7 @@ module qcmod public :: qc_gmi public :: qc_amsr2 public :: qc_saphir + ! set passed variables to public public :: npres_print,nlnqc_iter,varqc_iter,pbot,ptop,c_varqc,njqc,vqc,nvqc,hub_norm public :: use_poq7,noiqc,vadfile,dfact1,dfact,erradar_inflate,gps_jacqc @@ -200,6 +204,7 @@ module qcmod public :: troflg public :: lat_c public :: nrand + public :: airs_cads, cris_cads, iasi_cads logical nlnqc_iter,njqc,vqc,nvqc,hub_norm logical noiqc @@ -215,6 +220,7 @@ module qcmod logical vadwnd_l2rw_qc logical troflg logical cao_check + logical airs_cads, cris_cads, iasi_cads character(10):: vadfile integer(i_kind) npres_print @@ -455,6 +461,10 @@ subroutine init_qcvars lat_c=21.0_r_kind nrand=13 + airs_cads = .false. + cris_cads = .false. + iasi_cads = .false. + return end subroutine init_qcvars @@ -2065,10 +2075,11 @@ subroutine qc_saphir(nchanl,sfchgt,luse,sea, & return end subroutine qc_saphir -subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & - cris, hirs, zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, & - wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, & - id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole) +subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr,airs, & + cris,iasi,hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, & + wavenumber,ptau5,prsltmp,tvp,temp,wmix,chan_level,emissivity_k,ts,tsim, & + id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole,cluster_fraction, & + cluster_bt, chan_stdev, model_bt) ! id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole,radmod) ! all-sky !$$$ subprogram documentation block @@ -2108,6 +2119,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & ! tzbgr - Tz over water ! tsavg5 - surface skin temperature ! tbc - simulated - observed BT with bias correction +! tsim - simulated BT ! tb_obs - observed Brightness temperatures ! tnoise - channel noise array ! wavenumber - array of channel wavenumbers @@ -2133,6 +2145,10 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & ! cld - cloud fraction ! cldp - cloud pressure ! zero_irjaco3_pole - logical to control use of ozone jacobians near poles +! cluster_fraction - size of imager derived cluster to determine clear cloudy profiles, used by CADS +! cluster_bt - imager brightness temperature of each cluster, used by CADS +! chan_stdev - standard deviation of cluster mean temperatures, used by CADS +! model_bt _ brightness temperature derived from the model's clear profile. used by CADS ! ! attributes: ! language: f90 @@ -2142,11 +2158,13 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & use kinds, only: r_kind, i_kind use radinfo, only: iomg_det, itopo_det, isst_det + use crtm_planck_functions, only: crtm_planck_radiance + use cads, only: cloud_aerosol_detection implicit none ! Declare passed variables - logical, intent(in ) :: sea,land,ice,snow,luse,goessndr, cris, hirs + logical, intent(in ) :: sea,land,ice,snow,luse,goessndr,airs,cris,hirs,iasi logical, intent(inout) :: zero_irjaco3_pole integer(i_kind), intent(in ) :: nsig,nchanl,ndat,is integer(i_kind),dimension(nchanl), intent(in ) :: ich @@ -2157,10 +2175,14 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & real(r_kind), intent( out) :: cld,cldp real(r_kind),dimension(40,ndat), intent(inout) :: aivals real(r_kind),dimension(nchanl), intent(in ) :: tbc,emissivity_k,ts,wavenumber,tb_obs,tbcnob - real(r_kind),dimension(nchanl), intent(in ) :: tnoise + real(r_kind),dimension(nchanl), intent(in ) :: chan_level + real(r_kind),dimension(nchanl), intent(in ) :: tnoise,tsim real(r_kind),dimension(nsig,nchanl),intent(in ) :: ptau5,temp,wmix real(r_kind),dimension(nsig), intent(in ) :: prsltmp,tvp real(r_kind),dimension(nchanl), intent(inout) :: errf,varinv,varinv_use + real(r_kind),dimension(7), intent(in ) :: cluster_fraction + real(r_kind),dimension(2,7), intent(in ) :: cluster_bt + real(r_kind),dimension(2), intent(in ) :: chan_stdev, model_bt ! Declare local parameters @@ -2168,21 +2190,29 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & real(r_kind) :: demisf,dtempf,efact,dtbf,term,cenlatx,sfchgtfact - real(r_kind) :: sum,sum2,sum3,cloudp,tmp,dts,delta - real(r_kind),dimension(nchanl) :: dtb - integer(i_kind) :: i,j,k,kk,lcloud,m + real(r_kind) :: sum1,sum2,sum3,tmp,dts,delta + integer(i_kind) :: i,j,lcloud,m,isurface_chan integer(i_kind), dimension(nchanl) :: irday real(r_kind) :: dtz,ts_ave,xindx,tzchks real(r_kind),parameter:: tbmax = 550._r_kind real(r_kind),parameter:: tbmin = 50._r_kind +! for cloud_aerosol_detect + integer(i_kind) :: I_Sensor_ID + integer(i_kind),dimension(nchanl) :: chan_array, i_flag_cloud + integer(i_kind),dimension(2) :: imager_chans + integer(i_kind) :: boundary_layer_pres, tropopause_height + integer(i_kind) :: ichan_10_micron, ichan_12_micron + real(r_kind),dimension(nchanl) :: tb_bc + real(r_kind) :: cloud_temperature, radiance_chan, radiance_model, radiance_cloud + real(r_kind) :: tb_obs_10, tb_obs_12, tb_obs_diff ! Reduce weight given to obs for shortwave ir if ! solar zenith angle tiny_r_kind irday = 1 if (pangs <= 89.0_r_kind .and. frac_sea > zero) then ! QC2 in statsrad - if(luse)aivals(9,is) = aivals(9,is) + one + if(luse) aivals(9,is) = aivals(9,is) + one do i=1,nchanl if(wavenumber(i) > r2000)then if(wavenumber(i) > r2400)then @@ -2225,7 +2255,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & ! If GOES and lza > 60. do not use if( goessndr .and. zasat*rad2deg > r60) then ! QC5 in statsrad - if(luse)aivals(12,is) = aivals(12,is) + one + if(luse) aivals(12,is) = aivals(12,is) + one do i=1,nchanl varinv(i) = zero varinv_use(i)=zero @@ -2237,7 +2267,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & sfchgtfact=one if (zsges > r2000) then ! QC1 in statsrad - if(luse)aivals(8,is) = aivals(8,is) + one + if(luse) aivals(8,is) = aivals(8,is) + one sfchgtfact = (r2000/zsges)**4 endif @@ -2265,114 +2295,196 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & cld=zero cldp=r10*prsltmp(1) - do k=1,nsig - if(prsltmp(k) > trop5)then - do i=1,nchanl - dtb(i)=(tvp(k)-tsavg5)*ts(i) - end do - do kk=1,k-1 - do i=1,nchanl - dtb(i)=dtb(i)+(tvp(k)-tvp(kk))*temp(kk,i) - end do - end do - sum=zero - sum2=zero - do i=1,nchanl - if(varinv_use(i) > tiny_r_kind)then - sum=sum+tbc(i)*dtb(i)*varinv_use(i) - sum2=sum2+dtb(i)*dtb(i)*varinv_use(i) - end if - end do - if (abs(sum2) < tiny_r_kind) sum2 = sign(tiny_r_kind,sum2) - cloudp=min(max(sum/sum2,zero),one) - sum=zero - do i=1,nchanl - if(varinv_use(i) > tiny_r_kind)then - tmp=tbc(i)-cloudp*dtb(i) - sum=sum+tmp*tmp*varinv_use(i) - end if - end do - if(sum < sum3)then - sum3=sum - lcloud=k - cld=cloudp - cldp=r10*prsltmp(k) - end if - end if +! Cloud and aerosol detection routines (ECMWF) + if (cris .and. cris_cads) then + I_Sensor_ID = 27 + chan_array = nuchan(ich) ! channel numbers + tb_bc = tbc + tsim ! observation BT with bias correction + boundary_layer_pres = nint(0.8_r_kind*prsltmp(1)) ! boundary layer set to be 80% of surface pressure + tropopause_height = nint(trop5) + imager_chans = (/15,16/) ! imager channel numbers (from satinfo) + isurface_chan = 501 ! surface channel + ichan_10_micron = 458 ! ~10.7 micron channel for low level cloud test + ichan_12_micron = 295 ! ~12.0 micron channel for low level cloud test + + call cloud_aerosol_detection( I_Sensor_ID, nchanl, chan_array, & + tropopause_height, boundary_layer_pres, tb_bc, tsim, chan_level, imager_chans, cluster_fraction, & + cluster_bt, chan_stdev, model_bt, i_flag_cloud, cldp ) + + elseif ( iasi .and. iasi_cads ) then + I_Sensor_ID = 16 + chan_array = nuchan(ich) ! channel numbers + tb_bc = tbc + tsim ! observation BT with bias correction + boundary_layer_pres = nint(0.8_r_kind*prsltmp(1)) ! boundary layer set to be 80% of surface pressure + tropopause_height = nint(trop5) + imager_chans = (/2,3/) ! imager channel numbers (from satinfo) + isurface_chan = 1271 ! surface channel + ichan_10_micron = 1173 ! ~10.7 micron channel for low level cloud test + ichan_12_micron = 756 ! ~12.0 micron channel for low level cloud test + + call cloud_aerosol_detection( I_Sensor_ID, nchanl, chan_array, & + tropopause_height, boundary_layer_pres, tb_bc, tsim, chan_level, imager_chans, cluster_fraction, & + cluster_bt, chan_stdev, model_bt, i_flag_cloud, cldp ) + + elseif ( airs .and. airs_cads ) then + I_Sensor_ID = 11 + chan_array = nuchan(ich) ! channel numbers + tb_bc = tbc + tsim ! observation BT with bias correction + boundary_layer_pres = nint(0.8_r_kind*prsltmp(1)) ! boundary layer set to be 80% of surface pressure + tropopause_height = nint(trop5) + isurface_chan = 914 ! surface channel + imager_chans = (/0,0/) ! imager channel numbers (from satinfo) + ichan_10_micron = 843 ! ~10.7 micron channel for low level cloud test + ichan_12_micron = 587 ! ~12.0 micron channel for low level cloud test + + call cloud_aerosol_detection( I_Sensor_ID, nchanl, chan_array, & + tropopause_height, boundary_layer_pres, tb_bc, tsim, chan_level, imager_chans, cluster_fraction, & + cluster_bt, chan_stdev, model_bt, i_flag_cloud, cldp ) - end do - if ( lcloud > 0 ) then ! If cloud detected, reject channels affected by it. + else + call emc_legacy_cloud_detect(nchanl,nsig,tsavg5,trop5,prsltmp,tvp,ts,tbc,temp,varinv_use,lcloud,cld,cldp) - do i=1,nchanl + endif ! end of which cloud test to use +! compute cloud stats +! If using CADS + if ((cris .and. cris_cads) .or. (iasi .and. iasi_cads) .or. (airs .and. airs_cads)) then + +! Reject channels affected by clouds + do i=1, nchanl + if ( i_flag_cloud(i) == 1) then +! QC4 in statsrad + if(luse) aivals(11,is) = aivals(11,is) + one + varinv(i) = zero + varinv_use(i) = zero + if(id_qc(i) == igood_qc) id_qc(i) = ifail_cloud_qc + endif + end do + +! Derive cloud amount for CADS + cld = zero + if ( cldp < prsltmp(1) ) then ! if cloud in this profile exists + cloud_layer: do i=2, nsig ! determine which layer the cloud exists. + if (prsltmp(i) < cldp) then + lcloud = i + do j=1, nchanl ! use surface channel to derive cloud amount + m = nuchan(ich(j)) + if ( m == isurface_chan ) then ! interpolate cloud top temperature + cloud_temperature = ((tvp(lcloud) -tvp(lcloud -1)/ log(prsltmp(lcloud) / prsltmp(lcloud - 1))) & + * log(cldp/prsltmp(lcloud-1))) + tvp(lcloud-1) + call crtm_planck_radiance(1,m,tb_bc(j),radiance_chan) ! observation radiance. same as tb_obs + bias correction + call crtm_planck_radiance(1,m,tsim(j),radiance_model) ! model derived radiance + call crtm_planck_radiance(1,m,cloud_temperature,radiance_cloud) ! cloud top temperature radiance + cld = (radiance_chan - radiance_model) / (radiance_cloud - radiance_model) + cld = min(max(cld,zero),one) + cldp = cldp * r10 + exit cloud_layer ! cloud layer foound and cloud amount computed + endif ! surface channel found + end do !surface_chan + endif ! cloud found (prsltmp(i) < cldp) + end do cloud_layer + +! If clear, do a 10.7 - 12 micron test for low level clouds + else ! lcloud = 0 + do i=1, nchanl + if ( nuchan(ich(i)) == ichan_10_micron ) tb_obs_10 = tb_obs(i) + if ( nuchan(ich(i)) == ichan_12_micron ) tb_obs_12 = tb_obs(i) + end do + if ( tb_obs_10 > zero .and. tb_obs_12 > zero ) then + tb_obs_diff = tb_obs_10 - tb_obs_12 + if ( tb_obs_diff > 2.20_r_kind ) then ! Assume a cloud exists + cldp = prsltmp(1) * r10 ! Assume near surface cloud + cld = one ! Assume overcast cloud + lcloud = 1 + endif + endif + endif + +! If more than 2% of the transmittance comes from the cloud layer, reject the channel (0.02 is a tunable parameter). +! or CADS flagged a channel to have cloud. + if ( lcloud > 0 ) then + do i=1, nchanl + if ( ptau5(lcloud,i) > 0.02_r_kind ) then + if(luse) aivals(11,is) = aivals(11,is) + one ! QC4 in statsrad + varinv(i) = zero + varinv_use(i) = zero + if(id_qc(i) == igood_qc) id_qc(i) = ifail_cloud_qc + end if + end do + endif + +! default compute cloud stats, emc_legacy_cloud_detect + else + if ( lcloud > 0 ) then + + do i=1,nchanl ! reject channels with iuse_rad(j)=-1 when they are peaking below the cloud j=ich(i) if (passive_bc .and. iuse_rad(j)==-1) then - if (lcloud .ge. kmax(i)) then - if(luse)aivals(11,is) = aivals(11,is) + one - varinv(i) = zero - if(id_qc(i) == igood_qc)id_qc(i)=ifail_cloud_qc - cycle - end if + if (lcloud .ge. kmax(i)) then + if(luse)aivals(11,is) = aivals(11,is) + one + varinv(i) = zero + varinv_use(i) = zero + if(id_qc(i) == igood_qc)id_qc(i)=ifail_cloud_qc + cycle + end if end if ! If more than 2% of the transmittance comes from the cloud layer, ! reject the channel (0.02 is a tunable parameter) if ( ptau5(lcloud,i) > 0.02_r_kind) then -! QC4 in statsrad - if(luse)aivals(11,is) = aivals(11,is) + one - varinv(i) = zero - if(id_qc(i) == igood_qc)id_qc(i)=ifail_cloud_qc +! QC4 in statsrad + if(luse) aivals(11,is) = aivals(11,is) + one + varinv(i) = zero + varinv_use(i) = zero + if(id_qc(i) == igood_qc) id_qc(i) = ifail_cloud_qc end if - end do - -! If no clouds check surface temperature/emissivity + end do - else ! If no cloud was detected, do surface temp/emiss checks - sum=zero - sum2=zero - do i=1,nchanl + else ! surface consistency and sensitivity chacks. ( if lcoud = 0 ) + sum1=zero + sum2=zero + do i=1,nchanl if ( varinv_use(i) > tiny_r_kind .and. ts(i) > 0.0001_r_kind) then - sum=sum+tbc(i)*ts(i)*varinv_use(i) - sum2=sum2+ts(i)*ts(i)*varinv_use(i) + sum1 = sum1 +tbc(i)*ts(i)*varinv_use(i) + sum2 = sum2+ts(i)*ts(i)*varinv_use(i) endif - end do - if (abs(sum2) < tiny_r_kind) sum2 = sign(tiny_r_kind,sum2) - dts=abs(sum/sum2) - if(abs(dts) > one)then + end do + if (abs(sum2) < tiny_r_kind) sum2 = sign(tiny_r_kind,sum2) + dts=abs(sum1/sum2) + if(abs(dts) > one)then if(.not. sea)then - dts=min(dtempf,dts) + dts=min(dtempf,dts) else - dts=min(three,dts) + dts=min(three,dts) end if do i=1,nchanl - delta=max(r0_05*tnoise(i),r0_02) - if(abs(dts*ts(i)) > delta)then -! QC3 in statsrad - if(luse .and. varinv(i) > zero) aivals(10,is) = aivals(10,is) + one - varinv(i) = zero - if(id_qc(i) == igood_qc)id_qc(i)=ifail_sfcir_qc - end if - end do - end if - endif + delta=max(r0_05*tnoise(i),r0_02) + if(abs(dts*ts(i)) > delta)then +! QC3 in statsrad + if(luse .and. varinv(i) > zero) aivals(10,is) = aivals(10,is) + one + varinv(i) = zero + if(id_qc(i) == igood_qc)id_qc(i)=ifail_sfcir_qc + endif + enddo + endif + endif -! ! Temporary additional check for CrIS to reduce influence of land points on window channels (particularly important for bias correction) -! - if (cris .and. .not. sea) then - do i=1,nchanl - if (ts(i) > 0.2_r_kind) then + if (cris .and. .not. sea) then + do i=1,nchanl + if (ts(i) > 0.2_r_kind) then ! QC3 in statsrad - if(luse .and. varinv(i) > zero) aivals(10,is) = aivals(10,is) + one - varinv(i) = zero - if(id_qc(i) == igood_qc)id_qc(i)=ifail_sfcir_qc - end if - end do - end if - + if(luse .and. varinv(i) > zero) & + aivals(10,is) = aivals(10,is) + one + varinv(i) = zero + if(id_qc(i) == igood_qc) id_qc(i) = ifail_sfcir_qc + end if + end do + end if + endif ! derive cloud stats ! ! Apply Tz retrieval ! @@ -2402,7 +2514,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & cenlatx=abs(cenlat)*r0_04 if (cenlatx < one) then - if(luse)aivals(6,is) = aivals(6,is) + one + if(luse) aivals(6,is) = aivals(6,is) + one efact = half*(cenlatx+one) do i=1,nchanl if(varinv(i) > tiny_r_kind) errf(i)=efact*errf(i) @@ -2414,7 +2526,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & if(varinv(i) > tiny_r_kind)then dtbf = demisf*abs(emissivity_k(i))+dtempf*abs(ts(i)) term = dtbf*dtbf - if(term > tiny_r_kind)varinv(i)=varinv(i)/(one+varinv(i)*term) + if(term > tiny_r_kind) varinv(i) = varinv(i)/(one+varinv(i)*term) end if end do @@ -2497,12 +2609,113 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & endif !! if (hirs) !---mkim + return +end subroutine qc_irsnd +subroutine emc_legacy_cloud_detect(nchanl,nsig,tsavg5,trop5,prsltmp,tvp,ts,tbc,temp,varinv_use,lcloud,cld,cldp) + +!$$$ subprogram documentation block +! . . . +! subprogram: emc_legacy_cloud_detect determine clear/cloudy profiles from hirs,goessndr,airs,iasi,cris instruments +! +! prgmmr: derber ??? org: np23 date: ??? +! +! abstract: determine if a profile is clear/cloudy. If cloudy, determine model layer of the lcoud. +! This subroutine is designed for infrared sounders. +! +! program history log: +! 2022-06-20 jung moved into a subroutine +! +! input argument list: +! nchanl - number of channels per obs +! nsig - number of model layers +! tsavg5 - surface skin temperature +! trop5 - tropopause pressure +! prsltmp - array of layer pressure in vertical (surface to toa) +! tvp - array of temperatures in vertical (surface to toa) +! ts - skin temperature sensitivity +! tbc - simulated - observed BT with bias correction +! temp - temperature sensitivity array +! varinv_use - observation weight used (modified obs var error inverse) +! +! output argument list: +! lcloud - model layer of cloud +! cld - derived cloud amount +! cldp - model layer pressure (hPa) of cloud +! +! attributes: +! language: f90 +! machine: ibm RS/6000 SP +! +!$$$ end documentation block + +use kinds, only: i_kind, r_kind +use constants, only: tiny_r_kind, zero, r10 +implicit none + +integer(i_kind), intent(in ) :: nchanl, nsig +integer(i_kind), intent( out) :: lcloud +real(r_kind), intent(in ) :: tsavg5, trop5 +real(r_kind), intent( out) :: cld, cldp +real(r_kind), dimension(nchanl), intent(in ) :: tbc, ts, varinv_use +real(r_kind), dimension(nsig,nchanl), intent(in ) :: temp +real(r_kind), dimension(nsig), intent(in ) :: tvp, prsltmp + +integer(i_kind) :: i, k, kk + +real(r_kind) :: sum,sum2,sum3,cloudp,tmp +real(r_kind),dimension(nchanl) :: dtb + + sum3=zero + do i=1,nchanl + sum3=sum3+tbc(i)*tbc(i)*varinv_use(i) + end do + sum3=0.75_r_kind*sum3 + lcloud=0 + cld=zero + cldp=r10*prsltmp(1) + + do k=1,nsig + if(prsltmp(k) > trop5)then + do i=1,nchanl + dtb(i)=(tvp(k)-tsavg5)*ts(i) + end do + do kk=1,k-1 + do i=1,nchanl + dtb(i)=dtb(i)+(tvp(k)-tvp(kk))*temp(kk,i) + end do + end do + sum=zero + sum2=zero + do i=1,nchanl + if(varinv_use(i) > tiny_r_kind)then + sum=sum+tbc(i)*dtb(i)*varinv_use(i) + sum2=sum2+dtb(i)*dtb(i)*varinv_use(i) + end if + end do + if (abs(sum2) < tiny_r_kind) sum2 = sign(tiny_r_kind,sum2) + cloudp=min(max(sum/sum2,zero),one) + sum=zero + do i=1,nchanl + if(varinv_use(i) > tiny_r_kind)then + tmp=tbc(i)-cloudp*dtb(i) + sum=sum+tmp*tmp*varinv_use(i) + end if + end do + if(sum < sum3)then + sum3=sum + lcloud=k + cld=cloudp + cldp=r10*prsltmp(k) + end if + end if + + end do + +end subroutine emc_legacy_cloud_detect - return -end subroutine qc_irsnd subroutine qc_avhrr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & zsges,cenlat,frac_sea,pangs,trop5,tzbgr,tsavg5,tbc,tb_obs,tnoise, & @@ -2593,7 +2806,6 @@ subroutine qc_avhrr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & real(r_kind),parameter:: oneover400=1.0_r_kind/400.0_r_kind - real(r_kind) :: demisf,dtempf,efact,dtbf,term,cenlatx,sfchgtfact real(r_kind) :: sum1,sum2,sum3,cloudp,tmp,dts real(r_kind),dimension(nchanl,nsig) :: dtb @@ -4304,7 +4516,7 @@ subroutine qc_geocsr(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse, & integer(i_kind), intent(in ) :: nchanl,ndat,nsig,is integer(i_kind),dimension(nchanl),intent(in ) :: ich integer(i_kind),dimension(nchanl),intent(inout) :: id_qc - integer(i_kind),dimension(nchanl), intent(in ) :: kmax + integer(i_kind),dimension(nchanl),intent(in ) :: kmax real(r_kind), intent(in ) :: zsges real(r_kind), intent(in ) :: tzbgr real(r_kind),dimension(40,ndat), intent(inout) :: aivals diff --git a/src/gsi/read_airs.f90 b/src/gsi/read_airs.f90 index 8dd22ffd5e..c5392dad14 100644 --- a/src/gsi/read_airs.f90 +++ b/src/gsi/read_airs.f90 @@ -622,7 +622,6 @@ subroutine read_airs(mype,val_airs,ithin,isfcalc,rmesh,jsatid,gstime,& end do bufr_chans end if - ! Channel based quality control if(amsua)then diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index 9843f919d7..b8bf4ff92b 100644 --- a/src/gsi/read_cris.f90 +++ b/src/gsi/read_cris.f90 @@ -93,6 +93,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth,gsi_nstcoupler_deter use mpimod, only: npe use gsi_io, only: verbose + use qcmod, only: cris_cads ! use radiance_mod, only: rad_obs_type implicit none @@ -145,7 +146,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& character(len=40) :: infile2 integer(i_kind) :: kidsat, ksatid integer(i_kind) :: iret,ireadsb,ireadmg,irec,next, nrec_startx - integer(i_kind) :: bufr_nchan,maxinfo + integer(i_kind) :: bufr_nchan,maxinfo,dval_info integer(i_kind),allocatable,dimension(:)::nrec @@ -178,8 +179,8 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& logical :: outside,iuse,assim,valid,clear logical :: cris,quiet - integer(i_kind) :: ifov, ifor, iscn, instr, ioff, ilat, ilon, sensorindex - integer(i_kind) :: i, l, iskip, bad_line, llll + integer(i_kind) :: ifov, ifor, iscn, instr, ioff, ilat, ilon, sensorindex_cris + integer(i_kind) :: i, j, l, iskip, bad_line, llll integer(i_kind) :: nreal, isflg integer(i_kind) :: itx, k, nele, itt, n integer(i_kind):: idomsfc(1) @@ -187,7 +188,23 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind):: error_status, irecx,ierr integer(i_kind):: radedge_min, radedge_max integer(i_kind):: bufr_size - character(len=20),dimension(1):: sensorlist + character(len=20),allocatable,dimension(:) :: sensorlist + +! Imager cluster information for CADS + integer(i_kind) :: iexponent, sensorindex_imager, cads_info + integer(i_kind),dimension(7) :: imager_cluster_index + logical :: imager_coeff + logical,dimension(7) :: imager_cluster_flag + character(len=80) :: spc_filename + character(len=20) :: sensorlist_imager + real(r_kind),dimension(83,7) :: imager_info + real(r_kind),dimension(7) :: imager_cluster_size + real(r_kind),dimension(2) :: imager_mean, imager_std_dev, imager_conversion + real(r_kind) :: imager_cluster_tot + +! bufr error codes +! real(r_kind),dimension(7,3) :: error_codes + ! scan angle calculation geometry based on: ! C. Root 2014: JPSS Ground Project Code 474-00032 @@ -209,6 +226,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" integer(i_kind),parameter:: sfc_channel=501 !used in thinning routine if cloud informatino is not available + integer(i_kind),parameter:: band_2_start=714 !for CADS, if any of band 1 (chans 1 - 713) are missing, reject profile integer(i_kind),parameter:: ichan=-999 ! fov-based surface code is not channel specific for cris real(r_kind),parameter:: expansion=one ! exansion factor for fov-based surface code. ! use one for ir sensors. @@ -227,8 +245,12 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& maxinfo = 31 disterrmax=zero ntest=0 - if(dval_use) maxinfo = maxinfo + 2 - nreal = maxinfo + nstinfo + dval_info = 0 + if(dval_use) dval_info = 2 + cads_info = 0 + if(cris_cads) cads_info = 23 + nreal = maxinfo + cads_info + dval_info + nstinfo + ndata = 0 nodata = 0 cris= obstype == 'cris' .or. obstype == 'cris-fsr' @@ -301,46 +323,89 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& 'SAID YEAR MNTH DAYS HOUR MINU SECO CLATH CLONH SAZA BEARAZ SOZA SOLAZI' ! Load spectral coefficient structure - sensorlist(1)=sis quiet=.not. verbose + + imager_coeff = .false. +!TODO spc_filename = trim(crtm_coeffs_path)//'viirs-m_'//trim(jsatid)//'.SpcCoeff.bin' ! when viirs naming convention becomes standarized + if ( trim(jsatid) == 'npp' ) then + spc_filename = trim(crtm_coeffs_path)//'viirs-m_npp.SpcCoeff.bin' + sensorlist_imager = 'viirs-m_npp' + elseif ( trim(jsatid) == 'n20' ) then + spc_filename = trim(crtm_coeffs_path)//'viirs-m_n20.SpcCoeff.bin' + sensorlist_imager = 'viirs-m_n20' + inquire(file=trim(spc_filename), exist=imager_coeff) + if ( .not. imager_coeff ) spc_filename = trim(crtm_coeffs_path)//'viirs-m_j1.SpcCoeff.bin' + sensorlist_imager = 'viirs-m_j1' + elseif ( trim(jsatid) == 'n21' ) then + spc_filename = trim(crtm_coeffs_path)//'viirs-m_n21.SpcCoeff.bin' + sensorlist_imager = 'viirs-m_n21' + inquire(file=trim(spc_filename), exist=imager_coeff) + if ( .not. imager_coeff ) spc_filename = trim(crtm_coeffs_path)//'viirs-m_j2.SpcCoeff.bin' + sensorlist_imager = 'viirs-m_j2' + endif + inquire(file=trim(spc_filename), exist=imager_coeff) + if ( imager_coeff ) then + allocate( sensorlist(2)) + sensorlist(1) = sis +!TODO sensorlist(2) = 'viirs-m_'//trim(jsatid) !when viirs naming conventions becomes standardized + sensorlist(2) = trim(sensorlist_imager) + else + allocate( sensorlist(1)) + sensorlist(1) = sis + endif + if( crtm_coeffs_path /= "" ) then if(mype_sub==mype_root .and. print_verbose) write(6,*)'READ_CRIS: crtm_spccoeff_load() on path "'//trim(crtm_coeffs_path)//'"' error_status = crtm_spccoeff_load(sensorlist,& - File_Path = crtm_coeffs_path,quiet=quiet ) + File_Path = crtm_coeffs_path,quiet=quiet) else error_status = crtm_spccoeff_load(sensorlist,quiet=quiet) endif if (error_status /= success) then write(6,*)'READ_CRIS: ***ERROR*** crtm_spccoeff_load error_status=',error_status,& - ' TERMINATE PROGRAM EXECUTION' + ' TERMINATE PROGRAM EXECUTION' call stop2(71) endif +! find CRIS sensorindex. + sensorindex_cris = 0 + if ( sc(1)%sensor_id(1:4) == 'cris' )then + sensorindex_cris = 1 + else + write(6,*)'READ_CRIS: ***ERROR*** sensorindex_cris not set NO CRIS DATA USED' + write(6,*)'READ_CRIS: We are looking for ', sc(1)%sensor_id, ' TERMINATE PROGRAM EXECUTION' + call stop2(71) + end if + +! find imager sensorindex. + sensorindex_imager = 0 + if ( cris_cads .and. imager_coeff ) then + if ( sc(2)%sensor_id(1:4) == 'viir' )then + sensorindex_imager = 2 + else + write(6,*)'READ_CRIS: ***ERROR*** sensorindex_viirs not set NO VIIRS CLUSTER INFO USED BY CADS' + write(6,*)'READ_CRIS: We are looking for ', sc(2)%sensor_id, ' TERMINATE PROGRAM EXECUTION' + imager_coeff = .false. + end if + else + imager_coeff = .false. + end if + ! Find the channels being used (from satinfo file) in the spectral coef. structure. do i=subset_start,subset_end channel_number(i -subset_start +1) = nuchan(i) end do sc_index(:) = 0 satinfo_chan: do i=1,satinfo_nchan - spec_coef: do l=1,sc(1)%n_channels - if ( channel_number(i) == sc(1)%sensor_channel(l) ) then + spec_coef: do l=1,sc(sensorindex_cris)%n_channels + if ( channel_number(i) == sc(sensorindex_cris)%sensor_channel(l) ) then sc_index(i) = l exit spec_coef endif end do spec_coef end do satinfo_chan -! find CRIS sensorindex. - sensorindex = 0 - if ( sc(1)%sensor_id(1:4) == 'cris' )then - sensorindex = 1 - else - write(6,*)'READ_CRIS: ***ERROR*** sensorindex not set NO CRIS DATA USED' - write(6,*)'READ_CRIS: We are looking for ', sc(1)%sensor_id, ' TERMINATE PROGRAM EXECUTION' - call stop2(71) - end if - ! Calculate parameters needed for FOV-based surface calculation. if (isfcalc==1)then instr=17 ! CrIS is similar to AIRS for this purpose. @@ -687,7 +752,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& clear = .false. pred = zero -! Cloud information may be missing depending on how the VIIRS granules align +! Cloud information may be missing depending on how the imager granules align ! with the CrIS granules. ! Cloud Amount, TOCC is total cloud cover [%], HOCT is cloud height [m] call ufbint(lnbufr,cloud_properties,2,1,iret,'TOCC HOCT') @@ -699,7 +764,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& else pred1 = cloud_properties(2) *7.0_r_kind / r1000 ! Assume a lapse rate to convert hgt to delta TB. radiance = allchan(2,sfc_channel_index) * r1000 ! Conversion from W to mW - call crtm_planck_temperature(sensorindex,sfc_channel,radiance,temperature(sfc_channel_index)) ! radiance to BT calculation + call crtm_planck_temperature(sensorindex_cris,sfc_channel,radiance,temperature(sfc_channel_index)) ! radiance to BT calculation pred2 = tsavg *0.98_r_kind - temperature(sfc_channel_index) pred = max(pred1,pred2) ! use the largest of lapse rate (pred1) or sfc channel-surface difference (pred2) endif @@ -709,7 +774,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& ! the surface channel is fixed and set earlier in the code (501). radiance = allchan(2,sfc_channel_index) * r1000 ! Conversion from W to mW - call crtm_planck_temperature(sensorindex,sfc_channel,radiance,temperature(sfc_channel_index)) ! radiance to BT calculation + call crtm_planck_temperature(sensorindex_cris,sfc_channel,radiance,temperature(sfc_channel_index)) ! radiance to BT calculation if (temperature(sfc_channel_index) > tbmin .and. temperature(sfc_channel_index) < tbmax ) then if ( tsavg*0.98_r_kind <= temperature(sfc_channel_index)) then ! 0.98 is a crude estimate of the surface emissivity clear = .true. @@ -743,7 +808,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& ! now such spectra are rejected. if (( allchan(2,bufr_chan) > zero .and. allchan(2,bufr_chan) < 99999._r_kind)) then ! radiance bounds radiance = allchan(2,bufr_chan) * r1000 ! Conversion from W to mW - call crtm_planck_temperature(sensorindex,sc_chan,radiance,temperature(bufr_chan)) ! radiance to BT calculation + call crtm_planck_temperature(sensorindex_cris,sc_chan,radiance,temperature(bufr_chan)) ! radiance to BT calculation else ! error with channel number or radiance temperature(bufr_chan) = tbmin endif @@ -756,12 +821,12 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& bufr_chan = bufr_index(i) if(temperature(bufr_chan) <= tbmin .or. temperature(bufr_chan) >= tbmax ) then temperature(bufr_chan) = tbmin - if(iuse_rad(ioff+i) >= 0) iskip = iskip + 1 + if(iuse_rad(ioff+i) >= 0 .or. (cris_cads .and. sc_index(i) < band_2_start)) iskip = iskip + 1 endif end do skip_loop if(iskip > 0 .and. print_verbose)write(6,*) ' READ_CRIS : iskip > 0 ',iskip -! if( iskip >= 10 )cycle read_loop + if( iskip >= 10 .and. cris_cads ) cycle read_loop crit1=crit1 + ten*real(iskip,r_kind) @@ -772,9 +837,86 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& call finalcheck(one,crit1,itx,iuse) endif if(.not. iuse)cycle read_loop -! + +! Read the imager cluster information for the Cloud and Aerosol Detection Software. +! Only channels 15 and 16 are used. + + if ( cris_cads ) then + call ufbseq(lnbufr,imager_info,83,7,iret,'CRISCS') + if ( iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. & + imager_info(3,1) >= zero .and. imager_coeff ) then ! if imager cluster info exists + imager_mean = zero + imager_std_dev = zero + imager_cluster_tot = zero + imager_cluster_flag = .TRUE. + imager_cluster_size = imager_info(3,1:7) + imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:)) + imager_conversion(1) = one / (sc(sensorindex_imager)%wavenumber(4) **2) + imager_conversion(2) = one / (sc(sensorindex_imager)%wavenumber(5) **2) + +! Order clusters from largest (1) to smallest (7) + imager_cluster_sort: do i=1,7 + j = maxloc(imager_cluster_size,dim=1,mask=imager_cluster_flag) + imager_cluster_index(i) = j + imager_cluster_flag(j) = .FALSE. + end do imager_cluster_sort + +! Convert from radiance to brightness temperature for mean and standard devation used by CADS +! Imager cluster info added to data_all array. + + imager_cluster_info: do j=1,7 + i = imager_cluster_index(j) + + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction + imager_cluster_tot = imager_cluster_tot + imager_info(3,i) + + iexponent = -(nint(imager_info(75,i)) -11) ! channel 15 radiance for each cluster + imager_info(76,i) = imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) + + iexponent = -(nint(imager_info(77,i)) -11) ! channel 15 radiance std dev for each cluster. + imager_info(78,i) = imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) + + call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx)) + data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + + iexponent = -(nint(imager_info(80,i)) -11) ! channel 16 radiance for each cluster + imager_info(81,i) = imager_info(81,i) * imager_conversion(2) * (ten ** iexponent) + + iexponent = -(nint(imager_info(82,i))-5 ) ! channel 16 radiance std dev for each cluster. + iexponent = -(nint(imager_info(82,i)) -11) ! channel 16 radiance std dev for each cluster. + imager_info(83,i) = imager_info(83,i) * imager_conversion(2) * (ten ** iexponent) + + call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx)) + data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + + + end do imager_cluster_info + +! Compute cluster averages for each channel + + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:)) ! Channel 15 radiance cluster average + imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2 + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 15 radiance RMSE + call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) + call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1)) ! Channel 15 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 15 BT std dev + data_all(maxinfo+22,itx) = imager_std_dev(1) + + imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:)) ! Channel 16 radiance cluster average + imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2 + imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 16 radiance RMSE + call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) + call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2)) ! Channel 16 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 16 BT std dev + data_all(maxinfo+23,itx) = imager_std_dev(2) + + else ! Imager cluster information is missing. Set everything to zero + data_all(maxinfo+1 : maxinfo+25,itx) = zero + endif + endif ! cris_cads + ! interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr -! + if ( nst_gsi > 0 ) then tref = ts(0) dtw = zero @@ -818,15 +960,17 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& data_all(31,itx)= dlat_earth_deg ! earth relative latitude (degrees) if(dval_use) then - data_all(32,itx)= val_cris - data_all(33,itx)= itt + data_all(maxinfo+cads_info+1,itx)= val_cris + data_all(maxinfo+cads_info+2,itx)= itt +! data_all(32+cads_info,itx)= val_cris +! data_all(33+cads_info,itx)= itt end if if ( nst_gsi > 0 ) then - data_all(maxinfo+1,itx) = tref ! foundation temperature - data_all(maxinfo+2,itx) = dtw ! dt_warm at zob - data_all(maxinfo+3,itx) = dtc ! dt_cool at zob - data_all(maxinfo+4,itx) = tz_tr ! d(Tz)/d(Tr) + data_all(maxinfo+cads_info+dval_info+1,itx) = tref ! foundation temperature + data_all(maxinfo+cads_info+dval_info+2,itx) = dtw ! dt_warm at zob + data_all(maxinfo+cads_info+dval_info+3,itx) = dtc ! dt_cool at zob + data_all(maxinfo+cads_info+dval_info+4,itx) = tz_tr ! d(Tz)/d(Tr) endif ! Put satinfo defined channel temperatures into data array diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90 index 9ab1d5446f..edd7a9b50e 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -127,6 +127,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& use gsi_nstcouplermod, only: gsi_nstcoupler_skindepth, gsi_nstcoupler_deter use mpimod, only: npe use gsi_io, only: verbose + use qcmod, only: iasi_cads ! use radiance_mod, only: rad_obs_type implicit none @@ -208,11 +209,11 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& logical :: outside,iuse,assim,valid logical :: iasi,quiet,cloud_info - integer(i_kind) :: ifov, instr, iscn, ioff, sensorindex + integer(i_kind) :: ifov, instr, iscn, ioff, sensorindex_iasi integer(i_kind) :: i, j, l, iskip, ifovn, bad_line, ksatid, kidsat, llll integer(i_kind) :: nreal, isflg integer(i_kind) :: itx, k, nele, itt, n - integer(i_kind):: iexponent,maxinfo, bufr_nchan + integer(i_kind):: iexponent,maxinfo, bufr_nchan, dval_info integer(i_kind):: idomsfc(1) integer(i_kind):: ntest integer(i_kind):: error_status, irecx,ierr @@ -221,8 +222,18 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind) :: sfc_channel_index integer(i_kind),allocatable, dimension(:) :: channel_number, sc_index, bufr_index integer(i_kind),allocatable, dimension(:) :: bufr_chan_test - character(len=20),dimension(1):: sensorlist - + character(len=20),allocatable, dimension(:):: sensorlist + +! Imager clouser information for CADS + integer(i_kind) :: sensorindex_imager, cads_info + integer(i_kind),dimension(7) :: imager_cluster_index + logical :: imager_coeff + logical,dimension(7) :: imager_cluster_flag + character(len=80) :: spc_filename + real(r_kind),dimension(33,7) :: imager_info + real(r_kind),dimension(7) :: imager_cluster_size + real(r_kind),dimension(2) :: imager_mean, imager_std_dev + real(r_kind) :: imager_cluster_tot ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" @@ -248,8 +259,11 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& maxinfo = 31 disterrmax=zero ntest=0 - if(dval_use) maxinfo=maxinfo+2 - nreal = maxinfo + nstinfo + dval_info = 0 + if(dval_use) dval_info = 2 + cads_info = 0 + if(iasi_cads) cads_info = 23 + nreal = maxinfo + cads_info + dval_info + nstinfo ndata = 0 nodata = 0 @@ -315,7 +329,19 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& ! load spectral coefficient structure quiet=.not. verbose - sensorlist(1)=sis + + imager_coeff = .false. + spc_filename =trim(crtm_coeffs_path)//'avhrr3_'//trim(jsatid)//'.SpcCoeff.bin' + inquire(file=trim(spc_filename), exist=imager_coeff) + if ( imager_coeff ) then + allocate( sensorlist(2)) + sensorlist(1) = sis + sensorlist(2) = 'avhrr3_'//trim(jsatid) + else + allocate( sensorlist(1)) + sensorlist(1) = sis + endif + if( crtm_coeffs_path /= "" ) then if(mype_sub==mype_root .and. print_verbose) write(6,*)'READ_IASI: crtm_spccoeff_load() on path "'//trim(crtm_coeffs_path)//'"' error_status = crtm_spccoeff_load(sensorlist,& @@ -330,6 +356,31 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& call stop2(71) endif +! find IASI sensorindex + sensorindex_iasi = 0 + if ( sc(1)%sensor_id(1:4) == 'iasi' ) then + sensorindex_iasi = 1 + else + write(6,*)'READ_IASI: ***ERROR*** sensorindex_iasi not set NO IASI DATA USED' + write(6,*)'READ_IASI: We are looking for ', sc(1)%sensor_id, ' TERMINATE PROGRAM EXECUTION' + call stop2(71) + end if + +! find imager sensorindex + sensorindex_imager = 0 + if ( iasi_cads .and. imager_coeff ) then + if ( sc(2)%sensor_id(1:4) == 'avhr' ) then + sensorindex_imager = 2 + imager_coeff = .true. + else + write(6,*)'READ_IASI: ***ERROR*** sensorindex_imager is not set NO IASI DATA USED' + write(6,*)'READ_IASI: We are looking for ', sc(2)%sensor_id + imager_coeff = .false. + end if + else + imager_coeff = .false. + end if + ! Find the channels being used (from satinfo file) in the spectral coef. structure. do i=subset_start,subset_end channel_number(i -subset_start +1) = nuchan(i) @@ -337,23 +388,13 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& sc_index(:) = 0 satinfo_chan: do i=1,satinfo_nchan spec_coef: do l=1,sc(1)%n_channels - if ( channel_number(i) == sc(1)%sensor_channel(l) ) then + if ( channel_number(i) == sc(sensorindex_iasi)%sensor_channel(l) ) then sc_index(i) = l exit spec_coef endif end do spec_coef end do satinfo_chan -! find IASI sensorindex - sensorindex = 0 - if ( sc(1)%sensor_id(1:4) == 'iasi' ) then - sensorindex = 1 - else - write(6,*)'READ_IASI: sensorindex not set NO IASI DATA USED' - write(6,*)'READ_IASI: We are looking for ', sc(1)%sensor_id, ' TERMINATE PROGRAM EXECUTION' - call stop2(71) - end if - ! Calculate parameters needed for FOV-based surface calculation. if (isfcalc==1)then instr=18 @@ -725,7 +766,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& if (( allchan(2,bufr_chan) > zero .and. allchan(2,bufr_chan) < 99999._r_kind)) then ! radiance bounds radiance = allchan(2,bufr_chan)*scalef(bufr_chan) sc_chan = sc_index(i) - call crtm_planck_temperature(sensorindex,sc_chan,radiance,temperature(bufr_chan)) + call crtm_planck_temperature(sensorindex_iasi,sc_chan,radiance,temperature(bufr_chan)) else temperature(bufr_chan) = tbmin endif @@ -750,7 +791,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& ! crit1=crit1 + ten*real(iskip,r_kind) -! If the surface channel exists (~960.0 cm-1) and the AVHRR cloud information is missing, use an +! If the surface channel exists (~960.0 cm-1) and the imager cloud information is missing, use an ! estimate of the surface temperature to determine if the profile may be clear. if (.not. cloud_info) then pred = tsavg*0.98_r_kind - temperature(sfc_channel_index) @@ -766,6 +807,78 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& endif if(.not. iuse)cycle read_loop +! Read the imager cluster information for the Cloud and Aerosol Detection Software. +! Only channels 4 and 5 are used. + + if ( iasi_cads ) then + call ufbseq(lnbufr,imager_info,33,7,iret,'IASIL1CS') + if (iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. & + imager_info(3,1) >= zero .and. imager_coeff ) then ! if imager cluster info exists + imager_mean = zero + imager_std_dev = zero + imager_cluster_tot = zero + imager_cluster_flag = .TRUE. + imager_cluster_size = imager_info(3,1:7) + imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:)) + +! Order clusters from largest (1) to smallest (7) + imager_cluster_sort: do i=1,7 + j = maxloc(imager_cluster_size,dim=1,mask=imager_cluster_flag) + imager_cluster_index(i) = j + imager_cluster_flag(j) = .FALSE. + end do imager_cluster_sort + +! Convert from radiance to brightness temperature for mean and standard deviation used by CADS. +! Imager cluster info added to data_all array + + imager_cluster_info: do j=1,7 + i = imager_cluster_index(j) + + data_all(maxinfo+j,itx) = imager_cluster_size(i) ! Imager cluster fraction + imager_cluster_tot = imager_cluster_tot + imager_info(3,i) + + iexponent = -(nint(imager_info(25,i))-5 ) ! channel 4 radiance for each cluster. + imager_info(26,i) = imager_info(26,i) * (ten ** iexponent) + + iexponent = -(nint(imager_info(27,i))-5 ) ! channel 4 radiance std dev for each cluster. + imager_info(28,i) = imager_info(28,i) * (ten ** iexponent) + + call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx)) + data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero) + + iexponent = -(nint(imager_info(30,i))-5 ) ! channel 5 radiance for each cluster + imager_info(31,i) = imager_info(31,i) * (ten ** iexponent) + + iexponent = -(nint(imager_info(32,i))-5 ) ! channel 5 radiance std dev for each cluser. + imager_info(33,i) = imager_info(33,i) * (ten ** iexponent) + + call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx)) + data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero) + + end do imager_cluster_info + +! Compute cluster averages for each channel + + imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:)) ! Channel 4 radiance cluster average + imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(26,:)**2 + imager_info(28,:)**2)) - imager_mean(1)**2 + imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero)) ! Channel 4 radiance RMSE + call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1)) + call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1)) ! Channel 4 average BT + imager_std_dev(1) = imager_std_dev(1) - imager_mean(1) ! Channel 4 BT std dev + data_all(maxinfo+22,itx) = imager_std_dev(1) + + imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:)) ! Channel 5 radiance cluster average + imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(31,:)**2 + imager_info(33,:)**2)) - imager_mean(1)**2 + imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero)) ! Channel 5 radiance RMSE + call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2)) + call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2)) ! Channel 5 average BT + imager_std_dev(2) = imager_std_dev(2) - imager_mean(2) ! Channel 5 BT std dev + data_all(maxinfo+23,itx) = imager_std_dev(2) + + else ! Imager cluster information is missing. Set everything to zero + data_all(maxinfo+1 : maxinfo+25,itx) = zero + endif + endif ! iasi_cads = .true. ! ! interpolate NSST variables to Obs. location and get dtw, dtc, tz_tr ! @@ -813,15 +926,15 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& data_all(31,itx)= dlat_earth_deg ! earth relative latitude (degrees) if(dval_use)then - data_all(32,itx)= val_iasi - data_all(33,itx)= itt + data_all(maxinfo+cads_info+1,itx)= val_iasi + data_all(maxinfo+cads_info+2,itx)= itt end if if ( nst_gsi > 0 ) then - data_all(maxinfo+1,itx) = tref ! foundation temperature - data_all(maxinfo+2,itx) = dtw ! dt_warm at zob - data_all(maxinfo+3,itx) = dtc ! dt_cool at zob - data_all(maxinfo+4,itx) = tz_tr ! d(Tz)/d(Tr) + data_all(maxinfo+cads_info+dval_info+1,itx) = tref ! foundation temperature + data_all(maxinfo+cads_info+dval_info+2,itx) = dtw ! dt_warm at zob + data_all(maxinfo+cads_info+dval_info+3,itx) = dtc ! dt_cool at zob + data_all(maxinfo+cads_info+dval_info+4,itx) = tz_tr ! d(Tz)/d(Tr) endif ! Put satinfo defined channel temperatures into data array diff --git a/src/gsi/setupaod.f90 b/src/gsi/setupaod.f90 index 58707acd6a..e0964ee972 100644 --- a/src/gsi/setupaod.f90 +++ b/src/gsi/setupaod.f90 @@ -180,6 +180,7 @@ subroutine setupaod(obsLL,odiagLL,lunin,mype,nchanl,nreal,nobs,& real(r_kind) :: qcall, smask real(r_kind) :: styp, dbcf + real(r_kind),dimension(nchanl):: chan_level real(r_kind),dimension(nchanl):: emissivity,ts,emissivity_k real(r_kind),dimension(nchanl):: tsim real(r_kind),dimension(nsig,nchanl):: wmix,temp,ptau5 @@ -409,7 +410,7 @@ subroutine setupaod(obsLL,odiagLL,lunin,mype,nchanl,nreal,nobs,& call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,qsat,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & - tsim,emissivity,ptau5,ts,emissivity_k, & + tsim,emissivity,chan_level,ptau5,ts,emissivity_k, & temp,wmix,jacobian,error_status,layer_od=layer_od,jacobian_aero=jacobian_aero) ! interpolate aerosols at observation locations for diag files here if (aero_diagsave) then diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index 20ab63456e..1ff2474e20 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -297,12 +297,14 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& use radinfo, only: iland_det, isnow_det, iwater_det, imix_det, iice_det, & iomg_det, itopo_det, isst_det,iwndspeed_det, optconv use qcmod, only: setup_tzr_qc,ifail_scanedge_qc,ifail_outside_range + use qcmod, only: iasi_cads, cris_cads use state_vectors, only: svars3d, levels, svars2d, ns3d use oneobmod, only: lsingleradob,obchan,oblat,oblon,oneob_type use correlated_obsmod, only: corr_adjust_jacobian, idnames use radiance_mod, only: rad_obs_type,radiance_obstype_search,radiance_ex_obserr,radiance_ex_biascor use sparsearr, only: sparr2, new, writearray, size, fullarray use radiance_mod, only: radiance_ex_obserr_gmi,radiance_ex_biascor_gmi + use cads, only: cads_imager_calc implicit none @@ -400,6 +402,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& real(r_kind),dimension(nsig):: qvp,tvp,qs real(r_kind),dimension(nsig):: prsltmp real(r_kind),dimension(nsig+1):: prsitmp + real(r_kind),dimension(nchanl):: chan_level real(r_kind),dimension(nchanl):: weightmax real(r_kind),dimension(nchanl):: cld_rbc_idx,cld_rbc_idx2 real(r_kind),dimension(nchanl):: tcc @@ -440,6 +443,11 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& logical:: muse_ii +! variables added for CADS + real(r_kind),dimension(7,nobs) :: imager_cluster_fraction + real(r_kind),dimension(2,7,nobs) :: imager_cluster_bt + real(r_kind),dimension(2,nobs) :: imager_chan_stdev, imager_model_bt + ! Notations in use: for a single obs. or a single obs. type ! nchanl : a known channel count of a given type obs stream ! nchanl_diag : a subset of "iuse" @@ -591,6 +599,26 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& return endif +! Load data array for current satellite + read(lunin) data_s,luse,ioid + + if (nobskeep>0) then +! write(6,*)'setuprad: nobskeep',nobskeep + call stop2(275) + end if + + call dtime_setup() +! If using CADS setup arrays and calculate imager BTs + imager_cluster_fraction=zero + imager_cluster_bt=zero + imager_chan_stdev=zero + imager_model_bt=zero + if ((iasi_cads .and. iasi) .or. (cris_cads .and. cris)) then + + call cads_imager_calc(obstype,isis,nobs,nreal,nchanl,nsig,data_s,init_pass,mype, & + imager_cluster_fraction,imager_cluster_bt,imager_chan_stdev, imager_model_bt) + endif ! using cads + if ( mype == 0 .and. .not.l_may_be_passive) write(6,*)mype,'setuprad: passive obs',is,isis ! Logic to turn off print of reading coefficients if not first interation or not mype_diaghdr or not init_pass @@ -697,8 +725,6 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& endif endif - - ! Find number of channels written to diag file if(reduce_diag)then nchanl_diag=0 @@ -768,24 +794,9 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& if (netcdf_diag) call init_netcdf_diag_ endif -! Load data array for current satellite - read(lunin) data_s,luse,ioid - - if (nobskeep>0) then -! write(6,*)'setuprad: nobskeep',nobskeep - call stop2(275) - end if - - if (abi2km .and. regional) then - abi2km_bc = zero - abi2km_bc(2) = 233.5_r_kind - abi2km_bc(3) = 241.7_r_kind - abi2km_bc(4) = 250.5_r_kind - end if ! PROCESSING OF SATELLITE DATA ! Loop over data in this block - call dtime_setup() do n = 1,nobs ! Extract analysis relative observation time. dtime = data_s(itime,n) @@ -911,7 +922,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & - tsim,emissivity,ptau5,ts,emissivity_k, & + tsim,emissivity,chan_level,ptau5,ts,emissivity_k, & temp,wmix,jacobian,error_status,tsim_clr=tsim_clr,tcc=tcc, & tcwv=tcwv,hwp_ratio=hwp_ratio,stability=stability) if(gmi) then @@ -922,7 +933,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & - tsim2,emissivity2,ptau52,ts2,emissivity_k2, & + tsim2,emissivity2,chan_level,ptau52,ts2,emissivity_k2, & temp2,wmix2,jacobian2,error_status,tsim_clr=tsim_clr2,tcc=tcc,& tcwv=tcwv,hwp_ratio=hwp_ratio,stability=stability) ! merge @@ -946,7 +957,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & - tsim,emissivity,ptau5,ts,emissivity_k, & + tsim,emissivity,chan_level,ptau5,ts,emissivity_k, & temp,wmix,jacobian,error_status) if(gmi) then gmi_low_angles(1:3)=data_s(ilzen_ang:iscan_ang,n) @@ -956,7 +967,7 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& call call_crtm(obstype,dtime,data_s(:,n),nchanl,nreal,ich, & tvp,qvp,qs,clw_guess,ciw_guess,rain_guess,snow_guess,prsltmp,prsitmp, & trop5,tzbgr,dtsavg,sfc_speed, & - tsim2,emissivity2,ptau52,ts2,emissivity_k2, & + tsim2,emissivity2,chan_level,ptau52,ts2,emissivity_k2, & temp2,wmix2,jacobian2,error_status) ! merge emissivity(10:13) = emissivity2(10:13) @@ -1091,6 +1102,11 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& endif predbias=zero + abi2km_bc = zero + abi2km_bc(2) = 233.5_r_kind + abi2km_bc(3) = 241.7_r_kind + abi2km_bc(4) = 250.5_r_kind + !$omp parallel do schedule(dynamic,1) private(i,mm,j,k,tlap,node,bias) do i=1,nchanl mm=ich(i) @@ -1376,10 +1392,12 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& varinv_use(i) = zero end if end do - call qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse(n),goessndr, & - cris,hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, & - wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, & - id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole(n)) + + call qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse(n),goessndr,airs,cris,iasi, & + hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, & + wavenumber,ptau5,prsltmp,tvp,temp,wmix,chan_level,emissivity_k,ts,tsim, & + id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole(n), & + imager_cluster_fraction(:,n), imager_cluster_bt(:,:,n), imager_chan_stdev(:,n),imager_model_bt(:,n)) ! --------- MSU ------------------- ! QC MSU data @@ -2745,9 +2763,10 @@ subroutine contents_netcdf_diag_(odiags,idv,iob) deallocate(predbias_angord) endif end subroutine contents_netcdf_diag_ - subroutine final_binary_diag_ close(4) end subroutine final_binary_diag_ end subroutine setuprad + + end module rad_setup