! $Id: muri_abi_ocean_retrieval_mod.f90 2300 2017-08-29 10:33:06Z awalther $ module muri_abi_ocean_retrieval_mod use muri_definitions_mod, only: & muri_input_type & , muri_output_type use muri_abi_ocean_lut_mod, only: & ocean_lut use lib_array, only:interp1d use aw_lib_array, only: interp4d implicit none private public :: muri_abi_ocean_algorithm public :: sedimental contains ! ! ! subroutine muri_abi_ocean_algorithm (inp, out) implicit none type( muri_input_type), intent(in) :: inp type( muri_output_type),intent(out) :: out integer :: i_fm, i_cm integer :: n_opt integer, parameter :: N_FMR=11 integer, parameter :: N_BANDS = 6 real :: fmr real, allocatable:: refl_toa (:,:,:) integer :: i_cha,i_opt, i_fmr real, allocatable :: refl_reference (:) real :: aot_temp(N_FMR) real :: refl_corrsp(N_FMR,6) real :: temp(8) real :: aot_allbands(N_FMR,6) real :: n1 (N_FMR) real :: n2 real :: err(N_FMR,6) real :: err_sqrt (N_FMR) real :: aod_nleta_temp(4,5) real :: err_nleta_temp(4,5) real :: fmf_nleta_temp(4,5) integer :: idx(1), idx2(2) real :: val,val2 real :: aod_allbands(4,5,8) integer:: k,i,j integer :: CHANNEL_REFERENCE ! call inp%info CHANNEL_REFERENCE = 3 if ( inp% is_sedimental) CHANNEL_REFERENCE = 5 err_nleta_temp = huge(err_nleta_temp) err_sqrt = huge(err_sqrt) call ocean_lut % read_lut ( inp % sol,inp%sat,inp%azi,inp%ws, path = inp % path) call ocean_lut % sub_table (inp % sol,inp%sat,inp%azi,inp%ws) ! allocate with n_opt allocate ( refl_toa ( ocean_lut%N_OPT, N_FMR, N_BANDS)) allocate ( refl_reference ( ocean_lut% N_OPT)) n_opt = ocean_lut % n_opt !print*,'n_opt',ocean_lut % n_opt if (inp % rfl(CHANNEL_REFERENCE).GE.1.1*ocean_lut % refl_fine(CHANNEL_REFERENCE,1,1)) then ! - loop over fine and coarse mode do i_fm = 1, 4 do i_cm = 1, 5 ! loop over fine mode ratio do i_fmr = 1,N_FMR fmr = (i_fmr-1)/10. do i_opt = 1, N_OPT do i_cha = 1,n_BANDS refl_toa (i_opt,i_fmr,i_cha) = fmr * (ocean_lut % refl_fine(i_cha,i_opt,i_fm)) & & + (1 - fmr) * (ocean_lut % refl_coarse(i_cha,i_opt,i_cm)) end do end do !- reference channel is #3 refl_reference = refl_toa(:,i_fmr,CHANNEL_REFERENCE) aot_temp(i_fmr) = interp1d(refl_reference & , ocean_lut % aot_550nm & , inp % rfl(CHANNEL_REFERENCE) & , bounds_error = .false. & , FILL_VALUE = -999.) end do ! end do of (i_fmr) do i_cha = 1,6 do i_fmr =1,11 refl_corrsp(i_fmr,i_cha) = interp1d(ocean_lut % aot_550nm,refl_toa(:,i_fmr,i_cha) & , aot_temp(i_fmr),bounds_error = .false., FILL_VALUE = -999. ) end do end do do i_cha = 1,6 n1 = inp % rfl(i_cha) - refl_corrsp(:,i_cha) n2= inp % rfl(i_cha) - refl_toa(1,1,i_cha) + 0.01 if(n2.LT.0.01) n2=0.01 err(:,i_cha) = n1/n2 end do if ( inp% is_sedimental) then err_sqrt = sqrt ( ( err(:,3)**2 + err(:,5)**2.+ err(:,6)**2.)/3.) else !if (inp % rfl(CHANNEL_REFERENCE).GE.1.5*ocean_lut % refl_fine(CHANNEL_REFERENCE,1,1)) then err_sqrt = sqrt ( (err(:,2)**2. +err(:,3)**2.+ err(:,5)**2.+ err(:,6)**2.)/4) !else !err_sqrt = sqrt ( ( err(:,2)**2. + err(:,3)**2.+err(:,5)**2.)/3) !end if end if !print*,'error sqrt:', err_sqrt val= minval(err_sqrt) idx=minloc(err_sqrt) aod_nleta_temp(i_fm,i_cm) = aot_temp(idx(1)) err_nleta_temp(i_fm,i_cm) = val fmf_nleta_temp(i_fm,i_cm) = (idx(1) - 1 ) /10. ! print*,i_fm,i_cm,idx(1) !aod_allbands (i_fm,i_cm,:) = aot_allbands(idx(1),:) end do end do val2 = minval(err_nleta_temp) idx2 = minloc(err_nleta_temp) !- once find a good match between fwd and measurment give output out% aot = aod_nleta_temp(idx2(1),idx2(2)) out % fm_mode = idx2(1) out % cm_mode = idx2(2) ! print*,fmf_nleta_temp !stop out % fmf = fmf_nleta_temp(idx2(1),idx2(2)) !do i_cha =1,6 !out% aot_channel(i_cha) = aod_allbands(idx2(1),idx2(2),i_cha) !end do else ! if 0.86um is greather than 1.1 0.86um Rayleigh out% aot = -999. out % fm_mode = -1 out % cm_mode = -1 out % fmf = -1 end if end subroutine muri_abi_ocean_algorithm ! ! ! Sedimental test for ABI logical function sedimental (inp) type( muri_input_type), intent(in) :: inp integer :: i_cha real :: n= 0.0 real :: b real :: m real :: sumx = 0.0 real :: sumx2 = 0.0 real :: sumxy = 0.0 real :: sumy = 0.0 ! real :: sumy2 = 0.0 real :: x real ::y real :: b2_reference real :: delta_b2 real, parameter :: SWAV(4)= [0.47, 1.38, 1.61, 2.25] real :: Irfl(4) sedimental = .false. n=0.0 sumx=0.0 sumx2=0.0 sumxy=0.0 sumy=0.0 Irfl(1)=inp%rfl(1) Irfl(2)=inp%rfl(4) Irfl(3)=inp%rfl(5) Irfl(4)=inp%rfl(6) do i_cha = 1,4 n=n+1 sumx = sumx + ALOG(SWAV(i_cha)) sumx2 = sumx2 + ALOG(SWAV(i_cha)) * ALOG(SWAV(i_cha)) sumxy = sumxy + ALOG(SWAV(i_cha)) * ALOG(Irfl(i_cha)) sumy = sumy + ALOG(Irfl(i_cha)) ! sumy2 = sumy2 + ALOG(Irfl(i_cha)) * ALOG(Irfl(i_cha)) end do m = (n * sumxy - sumx * sumy )/ (n * sumx2 - sumx **2) b = (sumy * sumx2 - sumx * sumxy) / (n * sumx2 - sumx**2) ! b2 is 0.64 b2_reference = EXP(m*ALOG(0.64)+b) ! print*,'b2 refrence',b2_reference delta_b2=inp%rfl(2)-b2_reference ! print*,'delta b2',delta_b2 if(delta_b2.GT.0.05)then sedimental = .true. if ( inp % rfl(1).GE. 0.25 ) then sedimental = .false. end if end if ! print*,'sedimental',sedimental end function sedimental end module muri_abi_ocean_retrieval_mod