diff --git a/src/cloud_mask/ECM2/nbm_cloud_mask_lut_module.f90 b/src/cloud_mask/ECM2/nbm_cloud_mask_lut_module.f90 index 91153816f172ad6187635d56f2f6ce80c75c9733..45e7a8d8274ecda42701e296c4ab2ca391ca4663 100644 --- a/src/cloud_mask/ECM2/nbm_cloud_mask_lut_module.f90 +++ b/src/cloud_mask/ECM2/nbm_cloud_mask_lut_module.f90 @@ -142,134 +142,6 @@ end subroutine ECM2_CLOUD_MASK_COMPUTE_PRIOR !------------------------------------------------------------------- -subroutine CONVERT_R_TO_LOGLIKELIHOOD_1D(Table, Buffer_2d, water_prior, ice_prior, mode, MISSING) - real, dimension(:,:,:,:), intent(out) :: Table - real, dimension(:,:), intent(in) :: Buffer_2d - real, dimension(:), intent(in) :: water_prior, ice_prior - integer, intent(in) :: mode - real, intent(in) :: MISSING - real prior_prob, inv_prior_odds, inv_post_odds, post_prob, likelihood, r - integer :: i,j - Table(:,:,1,1) = MISSING - do i=1,size(Buffer_2d, 2) - prior_prob = -1.0 - if (mode .eq. 1) then - ! clear prior - if(water_prior(i) .ge. 0.0 .and. & - ice_prior(i) .ge. 0.0) then - prior_prob = 1.0 - water_prior(i) - ice_prior(i) - endif - else if (mode .eq. 2) then - if(water_prior(i) .ge. 0.0) then - prior_prob = water_prior(i) - endif - else - if(ice_prior(i) .ge. 0.0) then - prior_prob = ice_prior(i) - endif - endif - if(prior_prob .ge. 0.0) then - inv_prior_odds = (1.0-prior_prob)/prior_prob - do j=1,size(Buffer_2d, 1) - r = Buffer_2d(j,i) - if(r .ge. 0.0) then - inv_post_odds = r * inv_prior_odds - post_prob = 1.0 - inv_post_odds/(1.0+inv_post_odds) - likelihood = post_prob/prior_prob - Table(j,i,1,1) = log(likelihood) - endif - enddo - endif - enddo -end subroutine - -subroutine CONVERT_R_TO_LOGLIKELIHOOD_2D(Table, Buffer_3d, water_prior, ice_prior, mode, MISSING) - real, dimension(:,:,:,:), intent(out) :: Table - real, dimension(:,:,:), intent(in) :: Buffer_3d - real, dimension(:), intent(in) :: water_prior, ice_prior - integer, intent(in) :: mode - real, intent(in) :: MISSING - real prior_prob, inv_prior_odds, inv_post_odds, post_prob, likelihood, r - integer :: i,j,k - Table(:,:,:,1) = MISSING - do i=1,size(Buffer_3d, 3) - prior_prob = -1.0 - if (mode .eq. 1) then - ! clear prior - if(water_prior(i) .ge. 0.0 .and. & - ice_prior(i) .ge. 0.0) then - prior_prob = 1.0 - water_prior(i) - ice_prior(i) - endif - else if (mode .eq. 2) then - if(water_prior(i) .ge. 0.0) then - prior_prob = water_prior(i) - endif - else - if(ice_prior(i) .ge. 0.0) then - prior_prob = ice_prior(i) - endif - endif - if(prior_prob .ge. 0.0) then - inv_prior_odds = (1.0-prior_prob)/prior_prob - do j=1,size(Buffer_3d, 2) - do k=1,size(Buffer_3d, 1) - r = Buffer_3d(k,j,i) - if(r .ge. 0.0) then - inv_post_odds = r * inv_prior_odds - post_prob = 1.0 - inv_post_odds/(1.0+inv_post_odds) - likelihood = post_prob/prior_prob - Table(k,j,i,1) = log(likelihood) - endif - enddo - enddo - endif - enddo -end subroutine - -subroutine CONVERT_R_TO_LOGLIKELIHOOD_3D(Table, Buffer_4d, water_prior, ice_prior, mode, MISSING) - real, dimension(:,:,:,:), intent(out) :: Table - real, dimension(:,:,:,:), intent(in) :: Buffer_4d - real, dimension(:), intent(in) :: water_prior, ice_prior - integer, intent(in) :: mode - real, intent(in) :: MISSING - real prior_prob, inv_prior_odds, inv_post_odds, post_prob, likelihood, r - integer :: i,j,k,l - do i=1,size(Buffer_4d, 4) - prior_prob = -1.0 - if (mode .eq. 1) then - ! clear prior - if(water_prior(i) .ge. 0.0 .and. & - ice_prior(i) .ge. 0.0) then - prior_prob = 1.0 - water_prior(i) - ice_prior(i) - endif - else if (mode .eq. 2) then - if(water_prior(i) .ge. 0.0) then - prior_prob = water_prior(i) - endif - else - if(ice_prior(i) .ge. 0.0) then - prior_prob = ice_prior(i) - endif - endif - if(prior_prob .ge. 0.0) then - inv_prior_odds = (1.0-prior_prob)/prior_prob - do j=1,size(Buffer_4d, 3) - do k=1,size(Buffer_4d, 2) - do l=1,size(Buffer_4d, 1) - if(r .ge. 0.0) then - r = Buffer_4d(l,k,j,i) - inv_post_odds = r * inv_prior_odds - post_prob = 1.0 - inv_post_odds/(1.0+inv_post_odds) - likelihood = post_prob/prior_prob - Table(l,k,j,i) = log(likelihood) - endif - enddo - enddo - enddo - endif - enddo -end subroutine - subroutine ECM2_CLOUD_MASK_LUT_READ(Lut_File_Full_Path, N_Classifier) character(len=*), intent(in) :: Lut_File_Full_Path @@ -290,9 +162,6 @@ character(len=100) :: Var_Name logical :: File_Exist character(len=100) :: Temp_Name -real prior_prob, inv_prior_odds, inv_post_odds, post_prob, likelihood, r -integer :: i,j,k,l - ! --- initiate Is_Classifiers_Read = .FALSE. File_Exist = .FALSE. @@ -463,16 +332,13 @@ do Class_Idx = 1, N_Classifier Count_Read_2d(1) = Lut(Class_Idx)%Nbins_X Count_Read_2d(2) = Lut(Class_Idx)%N_Sfc call READ_NETCDF_2D_REAL(Group_Id, Start_Read_2d, (/1,1/), Count_Read_2d, 'clear_table', Buffer_2d) - call CONVERT_R_TO_LOGLIKELIHOOD_1D(Lut(Class_Idx)%Clear_Table, Buffer_2d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 1, MISSING) + Lut(Class_Idx)%Clear_Table(:,:,1,1) = Buffer_2d(:,:) call READ_NETCDF_2D_REAL(Group_Id, Start_Read_2d, (/1,1/), Count_Read_2d, 'water_table', Buffer_2d) - call CONVERT_R_TO_LOGLIKELIHOOD_1D(Lut(Class_Idx)%Water_Table, Buffer_2d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 2, MISSING) + Lut(Class_Idx)%Water_Table(:,:,1,1) = Buffer_2d(:,:) call READ_NETCDF_2D_REAL(Group_Id, Start_Read_2d, (/1,1/), Count_Read_2d, 'ice_table', Buffer_2d) - call CONVERT_R_TO_LOGLIKELIHOOD_1D(Lut(Class_Idx)%Ice_Table, Buffer_2d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 3, MISSING) + Lut(Class_Idx)%Ice_Table(:,:,1,1) = Buffer_2d(:,:) call READ_NETCDF_2D_REAL(Group_Id, Start_Read_2d, (/1,1/), Count_Read_2d, 'obs_table', Buffer_2d) Lut(Class_Idx)%Obs_Table(:,:,1,1) = Buffer_2d(:,:) @@ -494,16 +360,13 @@ do Class_Idx = 1, N_Classifier Count_Read_3d(2) = Lut(Class_Idx)%Nbins_Y Count_Read_3d(3) = Lut(Class_Idx)%N_Sfc call READ_NETCDF_3D_REAL(Group_Id, Start_Read_3d, (/1,1,1/), Count_Read_3d, 'clear_table', Buffer_3d) - call CONVERT_R_TO_LOGLIKELIHOOD_2D(Lut(Class_Idx)%Clear_Table, Buffer_3d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 1, MISSING) + Lut(Class_Idx)%Clear_Table(:,:,:,1) = Buffer_3d(:,:,:) call READ_NETCDF_3D_REAL(Group_Id, Start_Read_3d, (/1,1,1/), Count_Read_3d, 'water_table', Buffer_3d) - call CONVERT_R_TO_LOGLIKELIHOOD_2D(Lut(Class_Idx)%Water_Table, Buffer_3d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 2, MISSING) + Lut(Class_Idx)%Water_Table(:,:,:,1) = Buffer_3d(:,:,:) call READ_NETCDF_3D_REAL(Group_Id, Start_Read_3d, (/1,1,1/), Count_Read_3d, 'ice_table', Buffer_3d) - call CONVERT_R_TO_LOGLIKELIHOOD_2D(Lut(Class_Idx)%Ice_Table, Buffer_3d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 3, MISSING) + Lut(Class_Idx)%Ice_Table(:,:,:,1) = Buffer_3d(:,:,:) call READ_NETCDF_3D_REAL(Group_Id, Start_Read_3d, (/1,1,1/), Count_Read_3d, 'obs_table', Buffer_3d) Lut(Class_Idx)%Obs_Table(:,:,:,1) = Buffer_3d(:,:,:) @@ -527,17 +390,11 @@ do Class_Idx = 1, N_Classifier Count_Read_4d(3) = Lut(Class_Idx)%Nbins_Z Count_Read_4d(4) = Lut(Class_Idx)%N_Sfc call READ_NETCDF_4D_REAL(Group_Id, Start_Read_4d, (/1,1,1,1/), Count_Read_4d, 'clear_table', Buffer_4d) - call CONVERT_R_TO_LOGLIKELIHOOD_3D(Lut(Class_Idx)%Clear_Table, Buffer_4d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 1, MISSING) - + Lut(Class_Idx)%Clear_Table(:,:,:,:) = Buffer_4d(:,:,:,:) call READ_NETCDF_4D_REAL(Group_Id, Start_Read_4d, (/1,1,1,1/), Count_Read_4d, 'water_table', Buffer_4d) - call CONVERT_R_TO_LOGLIKELIHOOD_3D(Lut(Class_Idx)%Water_Table, Buffer_4d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 2, MISSING) - + Lut(Class_Idx)%Water_Table(:,:,:,:) = Buffer_4d(:,:,:,:) call READ_NETCDF_4D_REAL(Group_Id, Start_Read_4d, (/1,1,1,1/), Count_Read_4d, 'ice_table', Buffer_4d) - call CONVERT_R_TO_LOGLIKELIHOOD_3D(Lut(Class_Idx)%Ice_Table, Buffer_4d, & - Lut(Class_Idx)%Water_Fraction, Lut(Class_Idx)%Ice_Fraction, 3, MISSING) - + Lut(Class_Idx)%Ice_Table(:,:,:,:) = Buffer_4d(:,:,:,:) call READ_NETCDF_4D_REAL(Group_Id, Start_Read_4d, (/1,1,1,1/), Count_Read_4d, 'obs_table', Buffer_4d) Lut(Class_Idx)%Obs_Table(:,:,:,:) = Buffer_4d(:,:,:,:) diff --git a/src/cloud_mask/ECM2/nbm_cloud_mask_module.f90 b/src/cloud_mask/ECM2/nbm_cloud_mask_module.f90 index c17228fdc21f86965bf4707f4764bb830b4f75f3..281bb39661fed9bdbda421fa7ade1b4153d06356 100644 --- a/src/cloud_mask/ECM2/nbm_cloud_mask_module.f90 +++ b/src/cloud_mask/ECM2/nbm_cloud_mask_module.f90 @@ -118,15 +118,16 @@ MODULE ECM2_CLOUD_MASK_MODULE INTEGER, save:: Diag_Warning_Flag = 0 + 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 :: Posterior_Cld_Probability_By_Class REAL :: Post_Prob_Clear REAL :: Prior_Yes_Temp REAL :: Post_Prob_Clear_By_Class - REAL :: Post_Prob_Water_By_Class - REAL :: Post_Prob_Ice_By_Class REAL, PARAMETER :: Prob_Min_Thresh = 0.001 - REAL :: R_Max + REAL, PARAMETER :: R_Max = 1000.0 REAL :: Value_Dim REAL :: X_Dim, Y_Dim, Z_Dim REAL :: R_Clear @@ -273,10 +274,16 @@ MODULE ECM2_CLOUD_MASK_MODULE ! - ALLOCATE needed variables 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(Posterior_Cld_Probability_By_Class)) DEALLOCATE (Posterior_Cld_Probability_By_Class) ALLOCATE (Class_Use_Flag(Nclass)) + ALLOCATE (Clear_Cond_Ratio(Nclass)) + ALLOCATE (Water_Cond_Ratio(Nclass)) + ALLOCATE (Ice_Cond_Ratio(Nclass)) ALLOCATE (Obs_Prob(Nclass)) ALLOCATE (Posterior_Cld_Probability_By_Class(Nclass)) @@ -288,6 +295,9 @@ MODULE ECM2_CLOUD_MASK_MODULE Output%Posterior_Ice_Probability = MISSING_VALUE_REAL4 Output%Cld_Mask_Bayes = MISSING_VALUE_INT1 Output%Cld_Mask_Binary = MISSING_VALUE_INT1 + Clear_Cond_Ratio = MISSING_VALUE_REAL4 + Water_Cond_Ratio = MISSING_VALUE_REAL4 + Ice_Cond_Ratio = MISSING_VALUE_REAL4 Obs_Prob = MISSING_VALUE_REAL4 R_Clear = 1.0 @@ -787,28 +797,35 @@ MODULE ECM2_CLOUD_MASK_MODULE Input%Coastal_Mask, Output%Sfc_Idx, Class_Idx, Lut, & MISSING_VALUE_INT1, MISSING_VALUE_REAL4, 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 + !-------------------------------------- + ! 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 ( Z1 /= MISSING_VALUE_REAL4 & - .and. Z2 /= MISSING_VALUE_REAL4 & - .and. Z3 /= MISSING_VALUE_REAL4 ) then - R_Clear = R_Clear + Z1 - R_Water = R_Water + Z2 - R_Ice = R_Ice + Z3 - - Post_Prob_Clear_By_Class = exp(Z1)*Prior_Prob_Clear - Post_Prob_Water_By_Class = exp(Z2)*Prior_Prob_Water - Post_Prob_Ice_By_Class = exp(Z3)*Prior_Prob_Ice - - !--- normalize these probabilties so that they sum to 1 - Post_Sum = Post_Prob_Clear_By_Class + Post_Prob_Water_By_Class + Post_Prob_Ice_By_Class - Post_Prob_Clear_By_Class = Post_Prob_Clear_By_Class / Post_Sum - - ! --- make a cloud probability from the complement of the clear probability - Posterior_Cld_Probability_By_Class(Class_Idx) = 1.0 - Post_Prob_Clear_By_Class - ENDIF + IF (Clear_Cond_Ratio(Class_Idx) /= MISSING_VALUE_REAL4) & + R_Clear = min(R_Max, R_Clear * Clear_Cond_Ratio(Class_Idx)) + IF (Water_Cond_Ratio(Class_Idx) /= MISSING_VALUE_REAL4) & + R_Water = min(R_Max, R_Water * Water_Cond_Ratio(Class_Idx)) + IF (Ice_Cond_Ratio(Class_Idx) /= MISSING_VALUE_REAL4) & + R_Ice = min(R_Max, R_Ice * Ice_Cond_Ratio(Class_Idx)) + + !--- store class value + IF (Clear_Cond_Ratio(Class_Idx) /= MISSING_VALUE_REAL4) 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 + Posterior_Cld_Probability_By_Class(Class_Idx) = 1.0 - Post_Prob_Clear_By_Class + + ENDIF + + ENDDO ! loop over classifiers @@ -816,23 +833,23 @@ 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) ! NOTE ! if prior comes from LUT, then it comes from the last classifier in the LUT ! this could be bad, since filtering may cause distributional shift ! e.g. the last classifier is night-only +Post_Prob_Clear = 1.0 / (1.0 + R_Clear/Prior_Prob_Clear - R_Clear) +Output%Posterior_Water_Probability = 1.0 / (1.0 + R_Water/Prior_Prob_Water - R_Water) +Output%Posterior_Ice_Probability = 1.0 / (1.0 + R_Ice/Prior_Prob_Ice - R_Ice) -! Avoid overflow, make the max be zero -R_Max = max(R_Clear, R_Water, R_Ice) -R_Clear = R_Clear - R_Max -R_Water = R_Water - R_Max -R_Ice = R_Ice - R_Max - -Post_Prob_Clear = exp(R_Clear)*Prior_Prob_Clear -Output%Posterior_Water_Probability = exp(R_Water)*Prior_Prob_Water -Output%Posterior_Ice_Probability = exp(R_Ice)*Prior_Prob_Ice +!--- constrain cold clouds to be ice +if (Input%Topa < 240.0 .and. Input%Topa /= MISSING_VALUE_REAL4) then + Output%Posterior_Ice_Probability = Output%Posterior_Ice_Probability + Output%Posterior_Water_Probability + Output%Posterior_Water_Probability = 0.0 +endif -!--- normalize these probabilties so that they sum to 1 -Post_Sum = Post_Prob_Clear + Output%Posterior_Water_Probability + Output%Posterior_Ice_Probability +!--- normalize these probabilties so that they sum to 1 - won't happen IF prior is not from lut +Post_Sum = Post_Prob_clear + Output%Posterior_Water_Probability + Output%Posterior_Ice_Probability Post_Prob_Clear = Post_Prob_Clear / Post_Sum Output%Posterior_Water_Probability = Output%Posterior_Water_Probability / Post_Sum Output%Posterior_Ice_Probability = Output%Posterior_Ice_Probability / Post_Sum @@ -840,12 +857,6 @@ Output%Posterior_Ice_Probability = Output%Posterior_Ice_Probability / Post_Sum ! --- make a cloud probability from the complement of the clear probability Output%Posterior_Cld_Probability = 1.0 - Post_Prob_Clear -!--- constrain cold clouds to be ice -if (Input%Topa < 240.0 .and. Input%Topa /= MISSING_VALUE_REAL4) then - Output%Posterior_Ice_Probability = Output%Posterior_Ice_Probability + Output%Posterior_Water_Probability - Output%Posterior_Water_Probability = 0.0 -endif - ! --- make the binary mask IF (Output%Posterior_Cld_Probability < Mask_Thresh%Prob_Clear_Prob_Cloudy_Thresh(Output%Sfc_Idx) .and. & Output%Posterior_Cld_Probability /= MISSING_VALUE_REAL4) Output%Cld_Mask_Binary = symbol%CLEAR_BINARY @@ -952,6 +963,9 @@ ENDIF !IF Invalid_Data_Mask ! --- DEALLOCATE variables 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(Posterior_Cld_Probability_By_Class)) DEALLOCATE (Posterior_Cld_Probability_By_Class) IF (ALLOCATED(Chan_On)) DEALLOCATE (Chan_On)