!-------------------------------------------------------------------------------------- ! test Fortran implentation of new naive bayes cloud mask/phase using level2 file !-------------------------------------------------------------------------------------- program TEST_NEW_ECM_LEVEL2 use NB_CLOUD_MASK_SERVICES use READ_LEVEL2 use NB_CLOUD_MASK_LUT_MODULE use READ_PRIOR_MOD use CALC_PROB_MASK_PHASE, only: & GET_PROB_MASK_PHASE use NB_CLOUD_MASK_NETCDF_READ_MODULE, only: & OPEN_NETCDF_WRITE & , DEFINE_SDS_DIMS_NETCDF & , WRITE_SDS_NETCDF_2D_R4 & , WRITE_SDS_NETCDF_2D_I4 & , WRITE_SDS_NETCDF_3D_R4 & , CLOSE_NETCDF implicit none type(Classifier), dimension(:), allocatable:: Lut type(Level2) :: L2 real, dimension(:), allocatable :: Clear_Cond_Ratio real, dimension(:), allocatable :: Water_Cond_Ratio real, dimension(:), allocatable :: Ice_Cond_Ratio real, dimension(:), allocatable :: Obs_Prob real, dimension(:,:), allocatable :: Post_Prob_Cloud real, dimension(:,:,:), allocatable :: Post_Prob_Cloud_By_Class real, dimension(:,:), allocatable :: Post_Prob_Clear real, dimension(:,:), allocatable :: Post_Prob_Water real, dimension(:,:), allocatable :: Post_Prob_Ice real, dimension(:,:), allocatable :: Pix_Prior_Probability real :: Post_Prob_Clear_By_Class real, parameter :: MISSING = -999.0 real, parameter :: Prob_Min_Thresh = 0.001 real, parameter :: R_Max = 1000.0 real :: x, y, z real :: R_Clear real :: R_Water real :: R_Ice real :: Value_Dim real :: Prior_Prob_Clear, Prior_Prob_Water, Prior_Prob_Ice real :: Z1, Z2, Z3, Z4 real :: Post_Sum character(len=1020) :: Lut_File_Name_Full character(len=1020) :: Level2_Path character(len=1020) :: Level2_File character(len=1020) :: Level2_File_Name_Full character(len=1020) :: Prior_File_Name_Full character(len=1020) :: Output_Path character(len=1020) :: Output_File character(len=1020) :: Output_File_Name_Full character(len=30) :: Dim_Name integer, parameter :: Nchan = 45 integer :: Nclass integer :: Class_Idx integer :: Elem_Idx, Line_Idx integer :: i integer :: Sfc_Idx integer :: Ncid_Out integer :: Error integer :: Error_Io integer :: Is_Land integer :: Chan_Idx integer :: Indx_Str integer :: Number_Valid_Classes integer, dimension(3) :: Sds_Out_Dims integer, dimension(:,:), allocatable :: Cloud_Mask integer, dimension(:,:), allocatable :: Cloud_Mask_Binary integer, dimension(:,:), allocatable :: TUT integer, dimension(:,:), allocatable :: RUT integer, dimension(:), allocatable :: Chan_On integer, dimension(:), allocatable :: Chan_Wvl integer, dimension(:), allocatable :: Class_Use_Flag integer :: TUT_ROUTINE !function integer :: RUT_ROUTINE !function ! --- set pathes and file names !Lut_File_Name_Full = '/home/dbotambekov/idlscripts/idl_cloud_mask/fortran/test.nc' !Lut_File_Name_Full ='/home/dbotambekov/idlscripts/idl_cloud_mask/fortran/modis_1d_test.nc' !Lut_File_Name_Full ='/home/dbotambekov/idlscripts/idl_cloud_mask/fortran/try1.nc' Lut_File_Name_Full = '/apollo/cloud/personal/dbotambekov/ecm2/combined_tables/test_try10.nc' !Level2_File_Name_Full = '/apollo/cloud/scratch/Satellite_Output/VIIRS/global/ecm_prov/2018/level2/258/clavrx_npp_d20180915_t0000437_e0002079_b35662.level2.hdf' !Level2_Path = '/apollo/cloud/personal/dbotambekov/ecm2/level2/' !Level2_File = 'clavrx_MYD02SSH.A2015231.0000.006.2015231165148.level2.hdf' !Level2_File = 'clavrx_MYD02SSH.A2015231.0840.006.2015231214604.level2.hdf' !Level2_File = 'clavrx_MYD02SSH.A2015231.1510.006.2015232165058.level2.hdf' Level2_Path = '/apollo/cloud/personal/dbotambekov/ecm2/level2_test_goes16/' !Level2_File = 'clavrx_OR_ABI-L1b-RadF-M6C01_G16_s20192340000210.level2.hdf' !Level2_Path = '/apollo/cloud/scratch/Satellite_Output/NASA-SNPP_VIIRS/global/for_andy/2015/level2/231/' !Level2_File = 'clavrx_snpp_viirs.A2015231.0000.001.2017314133533.uwssec_B00019729.level2.hdf' !Level2_File_Name_Full = trim(Level2_Path) // trim(Level2_File) Prior_File_Name_Full = '/apollo/cloud/Ancil_Data/clavrx_ancil_data/static/luts/nb_cloud_mask/nb_cloud_mask_modis_prior.nc' Output_Path = '/apollo/cloud/personal/dbotambekov/ecm2/output_fortran_test_try10/' !Output_File_Name_Full = '/home/dbotambekov/idlscripts/idl_cloud_mask/fortran/out.nc' ! --- read input files open(unit=40,file ='input_list',status="old",action="read",iostat=Error_Io) ! - BEGIN LOOP OVER FILES File_Loop: do read(unit=40,fmt="(a)",iostat=Error_Io) Level2_File print *,'Running = ',trim(Level2_File) Level2_File_Name_Full = trim(Level2_Path) // trim(Level2_File) ! --- create output file name Indx_Str = index(trim(Level2_File),'level2.hdf') - 1 Output_File = trim(Level2_File(1:Indx_Str)) // 'fortran_ECM.nc' print *,'Output = ',trim(Output_File) Output_File_Name_Full = trim(Output_Path) // trim(Output_File) ! IF (Error_Io /= 0) EXIT !enddo File_Loop !close(unit=40) !stop !------------------------------------------------------------------------------ ! set channel on/off and channel mapping ! Chan_On = channel on / off binary flags ! Chan_Wvl = channel nominal wavelength in nm ! This is set up using CLAVR-x channel numbers !------------------------------------------------------------------------------ allocate(Chan_On(Nchan)) allocate(Chan_Wvl(Nchan)) Chan_On = 0 Chan_Wvl = 0 Chan_On(1) = 1 Chan_On(2) = 1 Chan_On(26) = 1 Chan_On(6) = 1 Chan_On(20) = 1 Chan_On(27) = 1 Chan_On(28) = 1 Chan_On(29) = 1 Chan_On(30) = 1 Chan_On(31) = 1 Chan_On(32) = 1 Chan_On(33) = 1 Chan_On(37) = 1 Chan_On(38) = 1 Chan_On(44) = 1 Chan_Wvl(1) = 650 Chan_Wvl(2) = 860 Chan_Wvl(26) = 1380 Chan_Wvl(6) = 1600 Chan_Wvl(20) = 3750 Chan_Wvl(27) = 6700 Chan_Wvl(28) = 7300 Chan_Wvl(29) = 8500 Chan_Wvl(30) = 9700 Chan_Wvl(31) = 11000 Chan_Wvl(32) = 12000 Chan_Wvl(33) = 13300 Chan_Wvl(37) = 6200 Chan_Wvl(38) = 10400 Chan_Wvl(44) = 700 !DNB !print *,'start, call NB_CLOUD_MASK_LUT_READ' !CALL system('date') !------------------------------------------------------------------- ! --- read lut !------------------------------------------------------------------- call NB_CLOUD_MASK_LUT_READ(trim(Lut_File_Name_Full), Lut, Nclass) !------------------------------------------------------------------- ! turn off tables based on chan_on !------------------------------------------------------------------- if (.not. allocated(Class_Use_Flag)) allocate(Class_Use_Flag(Nclass)) Class_Use_Flag = 1 do Class_Idx = 1, Nclass do Chan_Idx = 1, Nchan if (Chan_Idx > Lut(Class_Idx)%Nchan_Used) cycle if (Chan_On(Chan_Idx) == 0) then if(Chan_Wvl(Chan_Idx) == Lut(Class_Idx)%Wvl(Chan_Idx) .and. Chan_Wvl(Chan_Idx) /= 0) then Class_Use_Flag(Class_Idx) = 0 endif endif enddo !print *, Classifier_Names(Class_Idx), Class_Use_Flag(Class_Idx) enddo Number_Valid_Classes = sum(Class_Use_Flag) if (Number_Valid_Classes == 0) then print *, 'no active tests, Number_Valid_Classes =',Number_Valid_Classes stop endif !print *,'call READ_LEVEL2_FILE_DATA' !CALL system('date') !------------------------------------------------------------------- ! --- read level2 file !------------------------------------------------------------------- call READ_LEVEL2_FILE_DATA(trim(Level2_File_Name_Full), L2) !print *,'call COMPUTE_PRIOR' !CALL system('date') !------------------------------------------------------------------- ! --- compute prior !------------------------------------------------------------------- call COMPUTE_PRIOR(trim(Prior_File_Name_Full), L2%Longitude, L2%Latitude, L2%Month, Pix_Prior_Probability) !------------------------------------------------------------------- ! --- start calculating ECM !------------------------------------------------------------------- ! - allocate needed variables allocate (Clear_Cond_Ratio(Nclass)) allocate (Water_Cond_Ratio(Nclass)) allocate (Ice_Cond_Ratio(Nclass)) allocate (Obs_Prob(Nclass)) allocate (Post_Prob_Cloud(L2%Nx,L2%Ny)) allocate (Post_Prob_Cloud_By_Class(L2%Nx,L2%Ny,Nclass)) allocate (Post_Prob_Clear(L2%Nx,L2%Ny)) allocate (Post_Prob_Water(L2%Nx,L2%Ny)) allocate (Post_Prob_Ice(L2%Nx,L2%Ny)) allocate (Cloud_Mask_Binary(L2%Nx,L2%Ny)) allocate (Cloud_Mask(L2%Nx,L2%Ny)) allocate (TUT(L2%Nx,L2%Ny)) allocate (RUT(L2%Nx,L2%Ny)) ! - set to missing Post_Prob_Cloud = MISSING Post_Prob_Cloud_By_Class = MISSING Post_Prob_Clear = MISSING Post_Prob_Water = MISSING Post_Prob_Ice = MISSING Cloud_Mask = MISSING Cloud_Mask_Binary = MISSING TUT = MISSING RUT = MISSING do Class_Idx = 1, Nclass print *,Class_Idx-1,' ',Classifier_Names(Class_Idx) enddo !print *,'loop over pixels start' !CALL system('date') ! - loop over pixels do Elem_Idx = 1, L2%Nx do Line_Idx = 1, L2%Ny ! --- skip pixel if bad if (L2%Bad_Pixel_Mask(Elem_Idx,Line_Idx) == 1) cycle ! - reset table probabilities to missing Clear_Cond_Ratio = MISSING Water_Cond_Ratio = MISSING Ice_Cond_Ratio = MISSING Obs_Prob = MISSING R_Clear = 1.0 R_Water = 1.0 R_Ice = 1.0 ! - loop over classifiers do Class_Idx = 1, Nclass ! - set up input x = MISSING y = MISSING z = MISSING ! - loop over class rank do i = 1, Lut(Class_Idx)%Rank ! - select classifier name select case(i) case(1) Dim_Name = Lut(Class_Idx)%Class_Xname case(2) Dim_Name = Lut(Class_Idx)%Class_Yname case(3) Dim_Name = Lut(Class_Idx)%Class_Zname end select ! - select value based on the dimension name select case(trim(Dim_Name)) case('etropo10') Value_Dim = L2%Emiss_Tropo_10(Elem_Idx,Line_Idx) case('etropo11') Value_Dim = L2%Emiss_Tropo_11(Elem_Idx,Line_Idx) case('topa') Value_Dim = L2%Topa(Elem_Idx,Line_Idx) case('zopa') Value_Dim = L2%Zopa(Elem_Idx,Line_Idx) case('bt10') Value_Dim = L2%Bt10(Elem_Idx,Line_Idx) case('bt10std') Value_Dim = L2%Bt10Std(Elem_Idx,Line_Idx) case('logbt10std') Value_Dim = L2%Logbt10std(Elem_Idx,Line_Idx) case('bt11') Value_Dim = L2%Bt11(Elem_Idx,Line_Idx) case('bt11minsub') Value_Dim = L2%Bt11minsub(Elem_Idx,Line_Idx) case('bt11maxsub') Value_Dim = L2%Bt11maxsub(Elem_Idx,Line_Idx) case('dbt11max') Value_Dim = L2%Dbt11max(Elem_Idx,Line_Idx) case('dbt11maxsub') Value_Dim = L2%Dbt11maxsub(Elem_Idx,Line_Idx) case('bt11std') Value_Dim = L2%Bt11Std(Elem_Idx,Line_Idx) case('btdclr11') Value_Dim = L2%Bt11_Clr(Elem_Idx,Line_Idx) - L2%Bt11(Elem_Idx,Line_Idx) case('btdclr12') Value_Dim = L2%Bt12_Clr(Elem_Idx,Line_Idx) - L2%Bt12(Elem_Idx,Line_Idx) case('logbt11std') Value_Dim = L2%Logbt11std(Elem_Idx,Line_Idx) case('btd3810') Value_Dim = L2%Btd3810(Elem_Idx,Line_Idx) case('btd3811') Value_Dim = L2%Btd3811(Elem_Idx,Line_Idx) case('emiss3810') Value_Dim = L2%Emiss3810(Elem_Idx,Line_Idx) case('emiss3811') Value_Dim = L2%Emiss3811(Elem_Idx,Line_Idx) case('btd8511') Value_Dim = L2%Btd8511(Elem_Idx,Line_Idx) case('btd1112') Value_Dim = L2%Btd1112(Elem_Idx,Line_Idx) case('fmft') Value_Dim = L2%Fmft(Elem_Idx,Line_Idx) case('logcod065') Value_Dim = L2%Logcod065(Elem_Idx,Line_Idx) case('logcod138') Value_Dim = L2%Logcod138(Elem_Idx,Line_Idx) case('logcod160') Value_Dim = L2%Logcod160(Elem_Idx,Line_Idx) case('refl047') Value_Dim = L2%Ref047(Elem_Idx,Line_Idx) case('refl065') Value_Dim = L2%Ref065(Elem_Idx,Line_Idx) case('refl065clr') Value_Dim = L2%Ref065_Clr(Elem_Idx,Line_Idx) case('drefl065') Value_Dim = L2%Rgct(Elem_Idx,Line_Idx) case('drefl065clr') Value_Dim = L2%Dref065_Clr(Elem_Idx,Line_Idx) case('drefl065min') Value_Dim = L2%Dref065_Min(Elem_Idx,Line_Idx) case('drefl065minsub') Value_Dim = L2%Dref065_Min_Sub(Elem_Idx,Line_Idx) case('refratio') Value_Dim = L2%Refrat086065(Elem_Idx,Line_Idx) case('refrat086065') Value_Dim = L2%Refrat086065(Elem_Idx,Line_Idx) case('refrat138065') Value_Dim = L2%Refrat138065(Elem_Idx,Line_Idx) case('refl138') Value_Dim = L2%Ref138(Elem_Idx,Line_Idx) case('refl160') Value_Dim = L2%Ref160(Elem_Idx,Line_Idx) case('rgct') Value_Dim = L2%Rgct(Elem_Idx,Line_Idx) case('logrefl065std') Value_Dim = L2%Logref065Std(Elem_Idx,Line_Idx) case('ndsi') Value_Dim = L2%Ndsi(Elem_Idx,Line_Idx) case('bt6711covar') Value_Dim = L2%Bt6711Covar(Elem_Idx,Line_Idx) case default print *,"unknown classifier dimension name = ",trim(Dim_Name) Class_Use_Flag(Class_Idx) = 0 cycle end select ! - select input select case(i) case(1) x = Value_Dim case(2) y = Value_Dim case(3) z = Value_Dim end select enddo ! loop over class rank ! - select surface type index Sfc_Idx = int(L2%Sfc_Type_Mask(Elem_Idx,Line_Idx)) ! - set prior Prior_Prob_Clear = 1.0 - Pix_Prior_Probability(Elem_Idx,Line_Idx) Prior_Prob_Water = 0.5 * Pix_Prior_Probability(Elem_Idx,Line_Idx) Prior_Prob_Ice = 0.5 * Pix_Prior_Probability(Elem_Idx,Line_Idx) ! - compute for this class call GET_PROB_MASK_PHASE(x, y, z, L2%Zen(Elem_Idx,Line_Idx), & L2%Solzen(Elem_Idx,Line_Idx), L2%Lunzen(Elem_Idx,Line_Idx), & L2%Solglintzen(Elem_Idx,Line_Idx), L2%Lunglintzen(Elem_Idx,Line_Idx), & L2%Solscatang(Elem_Idx,Line_Idx), L2%Tpw(Elem_Idx,Line_Idx), & L2%Tsfc(Elem_Idx,Line_Idx), L2%Zsfc(Elem_Idx,Line_Idx), & L2%Land_Class(Elem_Idx,Line_Idx), L2%Snow_Class(Elem_Idx,Line_Idx), & L2%Solglint_Mask(Elem_Idx,Line_Idx), L2%Coast_Mask(Elem_Idx,Line_Idx), & Sfc_Idx, Class_Idx, Lut, MISSING, Z1, Z2, Z3, Z4) Clear_Cond_Ratio(Class_Idx) = Z1 Water_Cond_Ratio(Class_Idx) = Z2 Ice_Cond_Ratio(Class_Idx) = Z3 Obs_Prob(Class_Idx) = Z4 !remove this? akh !-------------------------------------- ! question - we really do not need to save the Cond_Ratio's for each Class !-------------------------------------- !--- combine this class into the running product of all classes if (Clear_Cond_Ratio(Class_Idx) /= MISSING) & R_Clear = min(R_Max, R_Clear * Clear_Cond_Ratio(Class_Idx)) if (Water_Cond_Ratio(Class_Idx) /= MISSING) & R_Water = min(R_Max, R_Water * Water_Cond_Ratio(Class_Idx)) if (Ice_Cond_Ratio(Class_Idx) /= MISSING) & R_Ice = min(R_Max, R_Ice * Ice_Cond_Ratio(Class_Idx)) !--- store class value if (Clear_Cond_Ratio(Class_Idx) /= MISSING) then Post_Prob_Clear_By_Class = 1.0 / (1.0 + Clear_Cond_Ratio(Class_Idx) / & Prior_Prob_Clear - Clear_Cond_Ratio(Class_Idx)) ! - convert from clear to cloud prob Post_Prob_Cloud_By_Class(Elem_Idx,Line_Idx,Class_Idx) = 1.0 - Post_Prob_Clear_By_Class endif enddo ! loop over classifiers ! --- compute posterior probs for this pixel ! - turn ratios into probabilities via post_prob = 1.0 / ( 1.0 + R/Prior_Yes - R) Post_Prob_Clear(Elem_Idx,Line_Idx) = 1.0 / (1.0 + R_Clear/Prior_Prob_Clear - R_Clear) Post_Prob_Water(Elem_Idx,Line_Idx) = 1.0 / (1.0 + R_Water/Prior_Prob_Water - R_Water) Post_Prob_Ice(Elem_Idx,Line_Idx) = 1.0 / (1.0 + R_Ice/Prior_Prob_Ice - R_Ice) !--- normalize these probabilties so that they sum to 1 - won't happen if prior is not from lut Post_Sum = Post_Prob_clear(Elem_Idx,Line_Idx) + & Post_Prob_Water(Elem_Idx,Line_Idx) + Post_Prob_Ice(Elem_Idx,Line_Idx) Post_Prob_Clear(Elem_Idx,Line_Idx) = Post_Prob_Clear(Elem_Idx,Line_Idx) / Post_Sum Post_Prob_Water(Elem_Idx,Line_Idx) = Post_Prob_Water(Elem_Idx,Line_Idx) / Post_Sum Post_Prob_Ice(Elem_Idx,Line_Idx) = Post_Prob_Ice(Elem_Idx,Line_Idx) / Post_Sum ! --- make a cloud probability from the complement of the clear probability Post_Prob_Cloud(Elem_Idx,Line_Idx) = 1.0 - Post_Prob_Clear(Elem_Idx,Line_Idx) ! --- for some reason Sfc_Idx could be = 0, this is a filter to avoid seg !!!!! fault, hard code it to be = 1 if (Sfc_Idx == 0) Sfc_Idx = 1 ! --- make the binary mask if (Post_Prob_Cloud(Elem_Idx,Line_Idx) < Prob_Clear_Prob_Cloudy_Thresh(Sfc_Idx) .and. & Post_Prob_Cloud(Elem_Idx,Line_Idx) /= MISSING) Cloud_Mask_Binary(Elem_Idx,Line_Idx) = 0 if (Post_Prob_Cloud(Elem_Idx,Line_Idx) >= Prob_Clear_Prob_Cloudy_Thresh(Sfc_Idx) .and. & Post_Prob_Cloud(Elem_Idx,Line_Idx) /= MISSING) Cloud_Mask_Binary(Elem_Idx,Line_Idx) = 1 ! --- make the 4-level mask if (Post_Prob_Cloud(Elem_Idx,Line_Idx) >= 0.0 .and. & Post_Prob_Cloud(Elem_Idx,Line_Idx) < Conf_Clear_Prob_Clear_Thresh(Sfc_Idx)) & Cloud_Mask(Elem_Idx,Line_Idx) = 0 if (Post_Prob_Cloud(Elem_Idx,Line_Idx) >= Conf_Clear_Prob_Clear_Thresh(Sfc_Idx) .and. & Post_Prob_Cloud(Elem_Idx,Line_Idx) < Prob_Clear_Prob_Cloudy_Thresh(Sfc_Idx)) & Cloud_Mask(Elem_Idx,Line_Idx) = 1 if (Post_Prob_Cloud(Elem_Idx,Line_Idx) >= Prob_Clear_Prob_Cloudy_Thresh(Sfc_Idx) .and. & Post_Prob_Cloud(Elem_Idx,Line_Idx) < Prob_Cloudy_Conf_Cloudy_Thresh(Sfc_Idx)) & Cloud_Mask(Elem_Idx,Line_Idx) = 2 if (Post_Prob_Cloud(Elem_Idx,Line_Idx) >= Prob_Cloudy_Conf_Cloudy_Thresh(Sfc_Idx) .and. & Post_Prob_Cloud(Elem_Idx,Line_Idx) <= 1.0) & Cloud_Mask(Elem_Idx,Line_Idx) = 3 !--- Thermal Uniformity Test Filter (Needs to Use 10.4 when necessary) TUT(ELem_Idx,Line_Idx) = TUT_ROUTINE(L2%Coast_Mask(Elem_Idx,Line_Idx), & L2%Bt11std(Elem_Idx,Line_Idx), & L2%Zsfcstd(Elem_Idx,Line_Idx), & Tut_Clear_Prob_Clear_Thresh(Sfc_Idx)) !--- Reflectance Uniformity Test Filter if (L2%Solzen(Elem_Idx,Line_Idx) < Lut(1)%Rut_Solzen_Thresh) then Is_Land = 1 if (L2%sfc_type_mask(Elem_Idx,Line_Idx) == 0 .or. & L2%Sfc_Type_Mask(Elem_Idx,Line_Idx) == 1) Is_Land = 0 RUT(Elem_Idx,Line_Idx) = RUT_ROUTINE(L2%Coast_Mask(Elem_Idx,Line_Idx), & Is_Land, & L2%Ref065(Elem_Idx,Line_Idx), & L2%Ref065Std(Elem_Idx,Line_Idx), & Rut_Clear_Prob_Clear_Thresh(Sfc_Idx)) endif !--- apply uniformity if (Cloud_Mask(Elem_Idx,Line_Idx) == 0 .and. & (TUT(Elem_Idx,Line_Idx) == 1 .or. & RUT(Elem_Idx,Line_Idx) == 1)) Cloud_Mask(Elem_Idx,Line_Idx) = 1 enddo ! loop over pixels enddo ! loop over pixels !print *,'loop over pixels end, write output' !CALL system('date') ! --- write to the output call OPEN_NETCDF_WRITE(Output_File_Name_Full, Ncid_Out) call DEFINE_SDS_DIMS_NETCDF(Ncid_Out, L2%Nx, L2%Ny, Nclass, Sds_Out_Dims, Error) call WRITE_SDS_NETCDF_2D_R4(Ncid_Out,L2%Nx,L2%Ny,(/Sds_Out_Dims(1),Sds_Out_Dims(2)/),'cloud_probability',Post_Prob_Cloud, Error) call WRITE_SDS_NETCDF_2D_I4(Ncid_Out,L2%Nx,L2%Ny,(/Sds_Out_Dims(1),Sds_Out_Dims(2)/),'cloud_mask',Cloud_Mask, Error) call WRITE_SDS_NETCDF_2D_I4(Ncid_Out,L2%Nx,L2%Ny,(/Sds_Out_Dims(1),Sds_Out_Dims(2)/),'cloud_mask_binary',Cloud_Mask_Binary, Error) call WRITE_SDS_NETCDF_3D_R4(Ncid_Out,L2%Nx,L2%Ny,Nclass,Sds_Out_Dims,'cloud_probability_by_class',Post_Prob_Cloud_By_Class, Error) call WRITE_SDS_NETCDF_2D_I4(Ncid_Out,L2%Nx,L2%Ny,(/Sds_Out_Dims(1),Sds_Out_Dims(2)/),'tut',TUT, Error) call WRITE_SDS_NETCDF_2D_I4(Ncid_Out,L2%Nx,L2%Ny,(/Sds_Out_Dims(1),Sds_Out_Dims(2)/),'rut',RUT, Error) !print *,L2%(1:3,1:3) call CLOSE_NETCDF(Ncid_Out) !print *,'write output done' !CALL system('date') ! --- deallocate everything if (allocated(Class_Use_Flag)) deallocate (Class_Use_Flag) if (allocated(Clear_Cond_Ratio)) deallocate (Clear_Cond_Ratio) if (allocated(Water_Cond_Ratio)) deallocate (Water_Cond_Ratio) if (allocated(Ice_Cond_Ratio)) deallocate (Ice_Cond_Ratio) if (allocated(Obs_Prob)) deallocate (Obs_Prob) if (allocated(Post_Prob_Cloud)) deallocate (Post_Prob_Cloud) if (allocated(Post_Prob_Cloud_By_Class)) deallocate (Post_Prob_Cloud_By_Class) if (allocated(Post_Prob_Clear)) deallocate (Post_Prob_Clear) if (allocated(Post_Prob_Water)) deallocate (Post_Prob_Water) if (allocated(Post_Prob_Ice)) deallocate (Post_Prob_Ice) if (allocated(Cloud_Mask_Binary)) deallocate (Cloud_Mask_Binary) if (allocated(Cloud_Mask)) deallocate (Cloud_Mask) if (allocated(TUT)) deallocate (TUT) if (allocated(RUT)) deallocate (RUT) if (allocated(Lut)) deallocate(Lut) if (allocated(L2%Bad_Pixel_Mask)) deallocate(L2%Bad_Pixel_Mask) if (allocated(L2%Latitude)) deallocate(L2%Latitude) if (allocated(L2%Longitude)) deallocate(L2%Longitude) if (allocated(L2%Land_Class)) deallocate(L2%Land_Class) if (allocated(L2%Snow_Class)) deallocate(L2%Snow_Class) if (allocated(L2%Coast_Mask)) deallocate(L2%Coast_Mask) if (allocated(L2%Zsfc)) deallocate(L2%Zsfc) if (allocated(L2%Zsfcstd)) deallocate(L2%Zsfcstd) if (allocated(L2%Tsfc)) deallocate(L2%Tsfc) if (allocated(L2%Tpw)) deallocate(L2%Tpw) if (allocated(L2%Zen)) deallocate(L2%Zen) if (allocated(L2%Solzen)) deallocate(L2%Solzen) if (allocated(L2%Solglintzen)) deallocate(L2%Solglintzen) if (allocated(L2%Solglint_Mask)) deallocate(L2%Solglint_Mask) if (allocated(L2%Solscatang)) deallocate(L2%Solscatang) if (allocated(L2%Lunzen)) deallocate(L2%Lunzen) if (allocated(L2%Lunglintzen)) deallocate(L2%Lunglintzen) if (allocated(L2%Sfc_Type_Mask)) deallocate(L2%Sfc_Type_Mask) if (allocated(L2%Topa)) deallocate(L2%Topa) if (allocated(L2%Zopa)) deallocate(L2%Zopa) if (allocated(L2%Emiss_Tropo_11)) deallocate(L2%Emiss_Tropo_11) if (allocated(L2%Bt10)) deallocate(L2%Bt10) if (allocated(L2%Bt10Std)) deallocate(L2%Bt10Std) if (allocated(L2%Logbt10Std)) deallocate(L2%Logbt10Std) if (allocated(L2%Bt11)) deallocate(L2%Bt11) if (allocated(L2%Bt11Std)) deallocate(L2%Bt11Std) if (allocated(L2%Bt11Max)) deallocate(L2%Bt11Max) if (allocated(L2%Logbt11Std)) deallocate(L2%Logbt11Std) if (allocated(L2%Bt11_Clr)) deallocate(L2%Bt11_Clr) if (allocated(L2%Logbt11std)) deallocate(L2%Logbt11std) if (allocated(L2%Bt11minsub)) deallocate(L2%Bt11minsub) if (allocated(L2%Bt11maxsub)) deallocate(L2%Bt11maxsub) if (allocated(L2%Dbt11max)) deallocate(L2%Dbt11max) if (allocated(L2%Dbt11maxsub)) deallocate(L2%Dbt11maxsub) if (allocated(L2%Bt12)) deallocate(L2%Bt12) if (allocated(L2%Bt12_Clr)) deallocate(L2%Bt12_Clr) if (allocated(L2%Btd8511)) deallocate(L2%Btd8511) if (allocated(L2%Btd3810)) deallocate(L2%Btd3810) if (allocated(L2%Btd3811)) deallocate(L2%Btd3811) if (allocated(L2%Bt6711Covar)) deallocate(L2%Bt6711Covar) if (allocated(L2%Emiss3810)) deallocate(L2%Emiss3810) if (allocated(L2%Emiss3811)) deallocate(L2%Emiss3811) if (allocated(L2%Btd1112)) deallocate(L2%Btd1112) if (allocated(L2%Rgct)) deallocate(L2%Rgct) if (allocated(L2%Fmft)) deallocate(L2%Fmft) if (allocated(L2%Ndsi)) deallocate(L2%Ndsi) if (allocated(L2%Ref047)) deallocate(L2%Ref047) if (allocated(L2%Ref065)) deallocate(L2%Ref065) if (allocated(L2%Ref065_Min)) deallocate(L2%Ref065_Min) if (allocated(L2%Ref065_Clr)) deallocate(L2%Ref065_Clr) if (allocated(L2%Dref065_Clr)) deallocate(L2%Dref065_Clr) if (allocated(L2%Dref065_Min)) deallocate(L2%Dref065_Min) if (allocated(L2%Ref065_Max_Sub)) deallocate(L2%Ref065_Max_Sub) if (allocated(L2%Dref065_Min_Sub)) deallocate(L2%Dref065_Min_Sub) if (allocated(L2%Logref065Std)) deallocate(L2%Logref065Std) if (allocated(L2%Ref086)) deallocate(L2%Ref086) if (allocated(L2%Ref138)) deallocate(L2%Ref138) if (allocated(L2%Ref160)) deallocate(L2%Ref160) if (allocated(L2%Bt38)) deallocate(L2%Bt38) if (allocated(L2%Bt85)) deallocate(L2%Bt85) if (allocated(L2%Ref065Std)) deallocate(L2%Ref065Std) if (allocated(L2%Refrat086065)) deallocate(L2%Refrat086065) if (allocated(L2%Refrat138065)) deallocate(L2%Refrat138065) if (allocated(L2%Logcod065)) deallocate(L2%Logcod065) if (allocated(L2%Logcod138)) deallocate(L2%Logcod138) if (allocated(L2%Logcod160)) deallocate(L2%Logcod160) if (allocated(L2%Cmask)) deallocate(L2%Cmask) if (allocated(L2%Cmask_Aux)) deallocate(L2%Cmask_Aux) if (allocated(Classifier_Names)) deallocate(Classifier_Names) if (allocated(Conf_Clear_Prob_Clear_Thresh)) deallocate(Conf_Clear_Prob_Clear_Thresh) if (allocated(Prob_Clear_Prob_Cloudy_Thresh)) deallocate(Prob_Clear_Prob_Cloudy_Thresh) if (allocated(Prob_Cloudy_Conf_Cloudy_Thresh)) deallocate(Prob_Cloudy_Conf_Cloudy_Thresh) if (allocated(Rut_Clear_Prob_Clear_Thresh)) deallocate(Rut_Clear_Prob_Clear_Thresh) if (allocated(Tut_Clear_Prob_Clear_Thresh)) deallocate(Tut_Clear_Prob_Clear_Thresh) deallocate(Chan_On) deallocate(Chan_Wvl) !print *,'END' !CALL system('date') ! --- END OF FILE INPUT LOOP IF (Error_Io /= 0) EXIT enddo File_Loop close(unit=40) end program TEST_NEW_ECM_LEVEL2 !==================================================================== ! Function Name: TUT_Routine ! ! Function: ! Thermal Uniformity Routine (TUT) Test ! ! Description: ! Computes the Thermal Uniformity Routine (TUT) test (YES/NO), as described in ! the ABI cloud mask ATBD ! ! ! Calling Sequence: ! Test_Results(Num_Tests,Elem_Idx,Line_Idx) = TUT_Routine (& ! Is_Land, & ! Is_Coast, & ! BT_Chn14_Stddev_3x3, & ! Sfc_Hgt_Stddev_3x3) ! ! Inputs: ! Is pixel land (YES/NO) ! Is pixel coast (YES/NO) ! Standard Deviation of the 11 micron BT over a 3x3 box ! Standard Deviation of the surface height over a 3x3 box ! ! Outputs: ! Function returns pass (sym%YES) or fail (sym%NO) result of the test via Test_Result ! ! Dependencies: ! Cloud Mask threshold include file with various needed thresholds. ! ! Restrictions: None ! !==================================================================== integer function TUT_ROUTINE (Is_Coast, & Bt_Chn14_Stddev_3x3, & Sfc_Hgt_Stddev_3x3, & BT_CHN14_CLR_UNI_THRESH) real, intent(in) :: Is_Coast real, intent(in) :: Bt_Chn14_Stddev_3x3 real, intent(in) :: Sfc_Hgt_Stddev_3x3 real, intent(in) :: BT_CHN14_CLR_UNI_THRESH real :: Test_Threshold TUT_ROUTINE = 0 Test_Threshold = 0.0 IF (Is_Coast == 0) THEN ! !7K/km is the adiabatic lapse rate ! Test_Threshold = 3.0 * 7.0*Sfc_Hgt_Stddev_3x3/1000.0 IF (Bt_Chn14_Stddev_3x3 > BT_CHN14_CLR_UNI_THRESH + Test_Threshold) THEN TUT_ROUTINE = 1 ENDIF ENDIF ! print, 'in TUT ', Is_Coast, Sfc_Hgt_Stddev_3x3, BT_Chn14_Stddev_3x3, BT_CHN14_CLR_UNI_THRESH, Test_Threshold, Test_Result RETURN end function TUT_ROUTINE !==================================================================== ! Function Name: RUT_Routine !==================================================================== integer function RUT_ROUTINE (Is_Coast, & Is_Land, & Refl_Chn2_Clear, & Refl_Chn2_Stddev_3x3, & Rut_Thresh) real, intent(in) :: Is_Coast integer, intent(in) :: Is_Land real, intent(in) :: Refl_Chn2_Clear real, intent(in) :: Refl_Chn2_Stddev_3x3 real, intent(in) :: RUT_THRESH real :: Test_Threshold RUT_ROUTINE = 0 IF (Is_Coast == 0) THEN !--- compute threshold IF (Is_Land == 1) THEN Test_Threshold = MAX(0.5,Refl_Chn2_Clear * Rut_Thresh) ELSE Test_Threshold = RUT_THRESH ENDIF !--- apply test IF (Refl_Chn2_Stddev_3x3 > Test_Threshold) THEN RUT_ROUTINE = 1 ENDIF ENDIF RETURN end function RUT_ROUTINE