! $Id$ ! module cloud_mask_mod implicit none integer , parameter :: et_class_FIRST = 0 integer , parameter:: et_class_T110 = 0 integer , parameter :: et_class_TMAX_T = 1 integer , parameter :: et_class_T_STD = 2 integer , parameter :: et_class_E_TROP = 3 integer , parameter :: et_class_FMFT = 4 integer , parameter :: et_class_BTD_110_067 = 5 integer , parameter :: et_class_BTD_110_067_COV = 6 integer , parameter :: et_class_BTD_110_085 = 7 integer , parameter :: et_class_E_037 = 8 integer , parameter :: et_class_E_037_DAY = 9 integer , parameter :: et_class_E_037_NGT = 10 integer , parameter :: et_class_BTD_037_110_NGT = 11 integer , parameter :: et_class_R_006_DAY = 12 integer , parameter :: et_class_R_STD = 13 integer , parameter :: et_class_R_006_MIN_3x3_DAY = 14 integer , parameter :: et_class_R_RATIO_DAY = 15 integer , parameter :: et_class_R_013_DAY = 16 integer , parameter :: et_class_R_016_DAY = 17 integer , parameter :: et_class_LAST = 17 type bayes_coef_type character (len =120) :: cvs_version character ( len =1000) :: file logical :: is_read = .false. integer :: n_class integer :: n_bounds real, dimension(:) , allocatable :: prior_no real, dimension(:) , allocatable :: prior_yes real, dimension(:) , allocatable :: Optimal_Posterior_Prob logical, dimension(:) , allocatable :: do_Skip_Sfc_Type real, dimension(:,:) , allocatable :: Classifier_Bounds_Min real, dimension(:,:) , allocatable :: Classifier_Bounds_Max real, dimension(:,:) , allocatable :: Delta_Classifier_Bounds real, dimension(:,:,:), allocatable :: class_cond_no real, dimension(:,:,:), allocatable :: class_cond_yes real, dimension(:,:,:), allocatable :: cond_no real, dimension(:,:,:), allocatable :: cond_yes character (len=20), dimension(:) , allocatable :: Classifier_Value_Name integer , dimension(:) , allocatable :: Classifier_Value_Name_enum integer , dimension(:), allocatable :: flag_idx contains procedure :: alloc => alloc_bayes_coef procedure :: dealloc => dealloc_bayes_coef end type bayes_coef_type type ( bayes_coef_type) , private , save :: bayes_coef contains ! ! ! subroutine cloud_mask ( conf , sat , sfc , imp , cld_mask ) use sat_mod use config_mod use sfc_data_mod use imp_mod use data_write_tools implicit none type (conf_user_opt_type ) , intent ( in) :: conf type (sat_type ) , intent ( in), target :: sat type (sfc_main_type ) , intent ( in) :: sfc type ( imp_type ) , intent ( in) :: imp type ( cr_data_type ) :: cld_mask integer , dimension(:,:), allocatable :: sfc_type_number integer :: elem_idx, line_idx , class_idx, sfc_idx real :: Classifier_Value real, dimension(:), allocatable :: Cond_yes real, dimension(:), allocatable :: Cond_no integer :: bin_idx logical :: is_on_test logical :: is_mountain = .true. logical :: has_cold_btd real :: post = -999. real , allocatable :: erg (:,:) character ( len =10) , allocatable :: names(:) !integer,dimension(:,:)::bayes_sfc_type ! - read in classifer coeffiecients bayes_coef % file = trim(conf%ancil_path)//'clavrx_ancil_data/naive_bayes_mask/' & //trim(conf % mask %bayesian_mask_classifier ) if ( .not. bayes_coef % is_read) call read_bayes_coeff ( ) allocate ( names (bayes_coef % n_class) ) names = bayes_coef %Classifier_Value_Name ! - determine sfc type allocate (sfc_type_number (sat% num_pix , sat%num_elem_this_seg)) allocate (erg (sat% num_pix , sat%num_elem_this_seg)) sfc_type_number = bayes_sfc_type ( sat % geo % lat , sat % geo % lat , sfc ) line_loop: do line_idx = 1, sat%num_elem_this_seg elem_loop: do elem_idx= 1,sat% num_pix allocate ( Cond_Yes ( bayes_coef % n_class ) ) allocate ( Cond_No ( bayes_coef % n_class ) ) sfc_idx = sfc_type_number(Elem_Idx,Line_Idx) is_mountain = sfc % z % data (Elem_Idx,Line_Idx) > 2000.0 & & .and. sfc_idx /= 6 has_cold_btd = .false. class_loop: do class_idx= 1 , bayes_coef % n_class ! - init is_on_test = .false. cond_yes ( class_idx) = 1.0 cond_no ( class_idx) = 1.0 ! - compute cond_yes cond_no for each test class select case ( bayes_coef % Classifier_Value_Name_enum (class_idx)) case( et_class_T110 ) if ( sat % data(31) % bt (Elem_Idx,Line_Idx) <= 0 ) cycle class_loop Classifier_Value = sat % data(31) % bt (Elem_Idx,Line_Idx) if ( sfc % land_class % data(Elem_Idx,Line_Idx) == ET_land_class % DEEP_OCEAN ) then Classifier_Value = imp % bt11_lrc (Elem_Idx,Line_Idx) end if is_on_test = .true. case( et_class_TMAX_T) if ( is_mountain ) cycle class_loop Classifier_Value = imp % bt_stats(31) % max(Elem_Idx,Line_Idx) - sat % data(31) % bt (Elem_Idx,Line_Idx) is_on_test = .true. case (et_class_T_STD) if ( is_mountain ) cycle class_loop Classifier_Value = imp % bt_stats(31) % std(Elem_Idx,Line_Idx) is_on_test = .true. case (et_class_E_TROP) ! - needs rtm computation of emis at tropopause(?) !Classifier_Value = Emiss_11um_Tropo_Rtm(Elem_Idx,Line_Idx) is_on_test = .false. case (et_class_FMFT) ! - needs rtm computation of bt_11_clear and bt_12_clear is_on_test = .false. case (et_class_BTD_110_067) if ( has_cold_btd ) cycle class_loop if ( sat % chan_on(27) .eqv. .false. ) cycle class_loop if ( .not. sat % chan_on(31) ) cycle class_loop Classifier_Value = sat % data(31) % bt (Elem_Idx,Line_Idx) - sat % data(27) % bt (Elem_Idx,Line_Idx) is_on_test = .true. case ( et_class_BTD_110_067_COV) is_on_test = .false. case( et_class_BTD_110_085) if ( .not. sat % chan_on(29) ) cycle class_loop if ( .not. sat % chan_on(31) ) cycle class_loop if ( has_cold_btd ) cycle class_loop Classifier_Value = sat % data(31) % bt (Elem_Idx,Line_Idx) - sat % data(29) % bt (Elem_Idx,Line_Idx) is_on_test = .true. case( et_class_E_037 ) if ( .not. sat % chan_on(20) ) cycle class_loop ! - solar contamin ! - oceanic_glint if ( sat % data(20) % bt (Elem_Idx,Line_Idx) <= 0 ) cycle class_loop if ( has_cold_btd ) cycle class_loop Classifier_Value = imp % emis_ch20(Elem_Idx,Line_Idx) is_on_test = .true. case( et_class_E_037_DAY) is_on_test = .false. case( et_class_E_037_NGT) is_on_test = .false. case( et_class_BTD_037_110_NGT) is_on_test = .false. case( et_class_R_006_DAY) is_on_test = .false. case( et_class_R_STD ) is_on_test = .false. case( et_class_R_006_MIN_3x3_DAY) is_on_test = .false. case( et_class_R_RATIO_DAY) is_on_test = .false. case( et_class_R_013_DAY) if ( .not. sat % chan_on(26) ) cycle class_loop Classifier_Value = sat % data(26) % ref (Elem_Idx,Line_Idx) is_on_test = .true. case ( et_class_R_016_Day) if ( .not. sat % chan_on(6) ) cycle class_loop Classifier_Value = sat % data(6) % ref (Elem_Idx,Line_Idx) is_on_test=.true. case default print*,'wrong class ', bayes_coef % Classifier_Value_Name_enum (class_idx) , bayes_coef % Classifier_Value_Name (class_idx) end select if ( is_on_test) then bin_idx = int (( classifier_value & & - bayes_coef % Classifier_Bounds_Min(Class_Idx,Sfc_Idx)) / & & bayes_coef % Delta_Classifier_Bounds(Class_Idx,Sfc_Idx)) + 1 bin_idx = max(1,min(bayes_coef % N_bounds-1,Bin_Idx)) Cond_yes (class_idx) = bayes_coef % class_cond_yes ( bin_idx, class_idx , sfc_idx ) Cond_no (class_idx) = bayes_coef % class_cond_no ( bin_idx, class_idx , sfc_idx ) end if end do class_loop ! - compute Posterior_Cld_Probability ! Posterior_Cld_Probability(Elem_Idx,Line_Idx) = & erg (Elem_Idx,Line_Idx) = & (bayes_coef % Prior_Yes(Sfc_Idx)*product(Cond_Yes)) / & (bayes_coef % Prior_Yes(Sfc_Idx)*product(Cond_Yes) + & bayes_coef % Prior_No(Sfc_Idx)*product(Cond_No)) deallocate ( Cond_Yes ) deallocate ( Cond_No ) end do elem_loop end do line_loop deallocate ( sfc_type_number ) call cld_mask % set_property ( data = erg ) end subroutine cloud_mask ! ++++++++++++ ! Reads the coeffeicients from file ! ++++++++++++++++ subroutine read_bayes_coeff ( ) use file_tools , only : getlun implicit none integer :: lun integer :: ios real :: time_diff_max integer :: n_class , n_bounds , n_sfc_bayes character (120) :: first_line character ( 72) :: header integer :: i_class integer :: i_sfc ! loop variable integer :: i_sfc_file ! number of surface in file integer :: int_dummy lun = getlun() open ( unit=lun, file = trim(bayes_coef % file) , & action ="read", form="formatted",status="old",iostat=ios) if (ios/= 0) then print*, 'no bayesain cloud mask ',trim(bayes_coef % file) return end if ! - first two lines are header meta-data read ( unit=lun, fmt="(a120)") first_line ! - in case we have an additional cvs id version line - very risky! ! if ( index (first_line, '$Id:') > 0 ) then bayes_coef % cvs_version = first_line read (unit=lun, fmt=*) header else bayes_coef % cvs_version= 'unknown' header = first_line end if read(unit=lun,fmt=*) read(unit=lun,fmt=*) time_diff_max, n_class, n_bounds, n_sfc_bayes bayes_coef % n_class = n_class bayes_coef % n_bounds = n_bounds call bayes_coef % alloc ( n_class, n_bounds, n_sfc_bayes ) do i_sfc = 1 , n_sfc_bayes read (unit=lun,fmt=*,iostat=ios) i_sfc_file, Header read(unit=lun,fmt=*,iostat=ios) bayes_coef % Prior_Yes(i_sfc), bayes_coef % Prior_No(i_sfc) read(unit=lun,fmt=*,iostat=ios) bayes_coef % Optimal_Posterior_Prob(i_Sfc) do i_class = 1 , n_class if ( i_sfc == 1) then read(unit=lun,fmt=*,iostat=ios) & int_dummy & , int_dummy & , int_dummy & , bayes_coef %Classifier_Value_Name(i_class) end if if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'T_11' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_T110 if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'T_max-T' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_TMAX_T if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'T_std' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_T_STD if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Emiss_tropo' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_E_TROP if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'FMFT' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_FMFT if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Btd_11_67' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_BTD_110_067 if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Bt_11_67_Covar' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_BTD_110_067_COV if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Btd_11_85' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_BTD_110_085 if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Emiss_375' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_E_037 if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Emiss_375_Day' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_E_037_DAY if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Emiss_375_Night' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_E_037_NGT if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Btd_375_11_Night' ) & & bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_BTD_037_110_NGT if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Ref_063_Day' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_R_006_DAY if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Ref_std' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_R_STD if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Ref_063_Min_3x3_Day' ) & & bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_R_006_MIN_3x3_DAY if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Ref_Ratio_Day' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_R_RATIO_DAY if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Ref_138_Day' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_R_013_DAY if ( trim(bayes_coef %Classifier_Value_Name(i_class)) == 'Ref_160_Day' ) bayes_coef %Classifier_Value_Name_enum(i_class) = et_class_R_016_DAY read(unit=lun,fmt=*,iostat=ios) & bayes_coef % Classifier_Bounds_Min(i_class,i_sfc) & , bayes_coef % Classifier_Bounds_Max(i_class,i_sfc) & , bayes_coef % Delta_Classifier_Bounds(i_class,i_sfc) read(unit=lun,fmt=*,iostat=ios) bayes_coef %Class_Cond_Yes(:,i_class,i_sfc) read(unit=lun,fmt=*,iostat=ios) bayes_coef %Class_Cond_No(:,i_class,i_sfc) end do end do bayes_coef % is_read = .true. end subroutine read_bayes_coeff ! ! ! 7 sfc bayesian surface types ! function bayes_sfc_type ( lat , lon, sfc ) use sfc_data_mod , only : sfc_main_type & , ET_land_class , ET_snow_class real , dimension ( : , : ) :: lat ,lon type ( sfc_main_type ) :: sfc integer ,dimension(size(lat,1),size(lat,2)):: bayes_sfc_type bayes_sfc_type = 0 ! # deep ocean where ( sfc % land_class % data == ET_land_class % DEEP_OCEAN) bayes_sfc_type = 1 end where ! # 2 - shallow ocean/water where ( sfc % land_class % data == ET_land_class % MODERATE_OCEAN .or. & sfc % land_class % data == ET_land_class % DEEP_INLAND_WATER .or. & sfc % land_class % data == ET_land_class % SHALLOW_INLAND_WATER .or. & sfc % land_class % data == ET_land_class % SHALLOW_OCEAN ) bayes_sfc_type = 2 end where ! # 3 unfrozen land where ( sfc % land_class % data == ET_land_class % LAND .or. & sfc % land_class % data == ET_land_class % COASTLINE .or. & sfc % land_class % data == ET_land_class % EPHEMERAL_WATER .or. & sfc % coast_mask % data == 1 ) bayes_sfc_type = 3 end where ! -TODO Sst_Back_Uni_Temp ! -- # 4 snow covered land where ( (lat > - 60.) .and. (sfc % snow_class % data == ET_snow_class % SNOW)) bayes_sfc_type = 4 end where ! -- # 5 Arctic where ( (lat >= 0.) .and. (sfc % snow_class % data == ET_snow_class % SEA_ICE )) bayes_sfc_type = 5 end where ! -- # 6 Antarctic where ( (lat <= -60.) .and. (sfc % snow_class % data == ET_snow_class % SNOW )) bayes_sfc_type = 6 end where where ( (lat <= -60.) .and. (sfc % snow_class % data == ET_snow_class % SEA_ICE )) bayes_sfc_type = 6 end where ! -- # 6 greenland where ( lat >= 60.0 & & .and. lon > -75. & & .and. lon < -10.0 & & .and. ( sfc % land_class % data == ET_land_class % LAND & .or. sfc % land_class % data == ET_land_class % COASTLINE ) & & .and. sfc % snow_class % data == ET_snow_class % SNOW ) bayes_sfc_type = 6 end where !-TODO ! -- # 7 Desert !-TODO end function bayes_sfc_type ! ! ! subroutine alloc_bayes_coef ( this, n_class , n_bounds , n_sfc_bayes ) class ( bayes_coef_type ) :: this integer :: n_class , n_bounds , n_sfc_bayes allocate(this%Prior_Yes(N_sfc_bayes)) allocate(this%Prior_No(N_sfc_bayes)) allocate(this%Optimal_Posterior_Prob(N_sfc_bayes)) allocate(this%do_Skip_Sfc_Type(N_sfc_bayes)) allocate(this%Classifier_Bounds_Min(N_class,N_sfc_bayes)) allocate(this%Classifier_Bounds_Max(N_class,N_sfc_bayes)) allocate(this%Delta_Classifier_Bounds(N_class,N_sfc_bayes)) allocate(this%Class_Cond_Yes(N_bounds-1,N_class,N_sfc_bayes)) allocate(this%Class_Cond_No(N_bounds-1,N_class,N_sfc_bayes)) allocate(this%Classifier_Value_Name(N_class)) allocate(this%Classifier_Value_Name_enum(N_class)) allocate(this%Flag_Idx(N_class)) end subroutine alloc_bayes_coef ! ! ! subroutine dealloc_bayes_coef ( this ) class ( bayes_coef_type ) :: this deallocate(this%Prior_Yes) deallocate(this%Prior_No) deallocate(this%Optimal_Posterior_Prob) deallocate(this%do_Skip_Sfc_Type) deallocate(this%Classifier_Bounds_Min) deallocate(this%Classifier_Bounds_Max) deallocate(this%Delta_Classifier_Bounds) deallocate(this%Class_Cond_Yes) deallocate(this%Class_Cond_No) deallocate(this%Classifier_Value_Name) deallocate(this%Classifier_Value_Name_enum) deallocate(this%Flag_Idx) end subroutine dealloc_bayes_coef !------------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------------- real elemental function split_window_test( t11_clear, t12_clear, t11, t12) real, intent(in):: t11_clear real, intent(in):: t12_clear real, intent(in):: t11 real, intent(in):: t12 split_window_test = (t11_clear - t12_clear) * (t11 - 260.0) / (t11_clear - 260.0) if (t11_clear <=265.0) split_window_test = 0.0 split_window_test = (t11 - t12) - split_window_test end function split_window_test end module cloud_mask_mod