var/da/da_radiance/mod_clddet_geoir.f90

   1   module mod_clddet_geoir
   2 !   use netcdf
   3 !  use mod_para
   4   use da_control, only: missing_r
   5   implicit none
   6   contains
   7
   8   subroutine qc_SDob(nlongitude,nlatitude,tbb,SDob)
   9    ! ------------QC: step 1 ----
  10    ! Kozo Okamoto, QJ, 2017 and Zhuge and Zou, JAMC, 2016
  11    ! 10.4wm(Okamoto, QJ, 2017),Chan13, 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14
  12    ! only 10.8 for AGRI, !!!cha12!!!
  13     implicit none
  14     integer, intent(in):: nlongitude,nlatitude
  15     real, intent(in)    :: tbb(nlongitude,nlatitude)
  16     real, intent(inout) :: SDob(nlongitude,nlatitude)
  17 ! temp
  18     integer, parameter :: npixs_sdob = 3*3    ! pixles
  19     real :: temp_tb(npixs_sdob), temp_std
  20     integer :: i, j
  21
  22     SDob=missing_r
  23 !    call SD_ob(nlongitude,nlatitude,tbb, SDob)
  24      do i=2, nlongitude-1
  25        do j=2, nlatitude-1
  26          temp_tb(1) = tbb(i-1,j-1)
  27          temp_tb(2) = tbb(i,  j-1)
  28          temp_tb(3) = tbb(i+1,j-1)
  29          temp_tb(4) = tbb(i-1,j)
  30          temp_tb(5) = tbb(i,  j)
  31          temp_tb(6) = tbb(i+1,j)
  32          temp_tb(7) = tbb(i-1,j+1)
  33          temp_tb(8) = tbb(i,  j+1)
  34          temp_tb(9) = tbb(i+1,j+1)
  35          call find_std(nlongitude,nlatitude,npixs_sdob, temp_tb, temp_std, missing_r, missing_r)
  36          SDob(i,j) = temp_std
  37        end do
  38      end do
  39
  40    end subroutine
  41
  42   subroutine qc_RTCT(nlongitude,nlatitude,tbb,ter,rtct,rtct2)
  43    ! ------------QC: step 2 ----
  44    ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14
  45    ! Relative thermal Contrast test
  46    ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14, and terrain
  47    ! only 10.8 for AGRI, !!!cha12!!!
  48     implicit none
  49     integer, intent(in):: nlongitude,nlatitude
  50     real, intent(in)    :: tbb(nlongitude,nlatitude)
  51     real, intent(in)    :: ter(nlongitude,nlatitude)
  52     real, intent(out) :: rtct(nlongitude,nlatitude)
  53     real, intent(inout) :: rtct2(nlongitude,nlatitude)
  54 !    real, intent(inout) :: rtct2(nlongitude,nlatitude)
  55 ! temp
  56     integer, parameter :: npixs_rtct = 3*3    ! pixles
  57     real :: temp_ter(npixs_rtct), temp_tb(npixs_rtct), temp_std
  58     integer :: i, j, m, n
  59
  60     rtct=missing_r  ! all missing_r for variables are same
  61 !    call SD_ob(nlongitude,nlatitude,tbb, SDob)
  62      do i=2, nlongitude-1
  63        do j=2, nlatitude-1
  64          temp_ter(1) = ter(i-1,j-1)
  65          temp_ter(2) = ter(i,  j-1)
  66          temp_ter(3) = ter(i+1,j-1)
  67          temp_ter(4) = ter(i-1,j  )
  68          temp_ter(5) = ter(i,  j  )
  69          temp_ter(6) = ter(i+1,j  )
  70          temp_ter(7) = ter(i-1,j+1)
  71          temp_ter(8) = ter(i,  j+1)
  72          temp_ter(9) = ter(i+1,j+1)
  73          temp_tb(1) = tbb(i-1,j-1 )
  74          temp_tb(2) = tbb(i,  j-1 )
  75          temp_tb(3) = tbb(i+1,j-1 )
  76          temp_tb(4) = tbb(i-1,j   )
  77          temp_tb(5) = tbb(i,  j   )
  78          temp_tb(6) = tbb(i+1,j   )
  79          temp_tb(7) = tbb(i-1,j+1 )
  80          temp_tb(8) = tbb(i,  j+1 )
  81          temp_tb(9) = tbb(i+1,j+1 )
  82          m = COUNT( temp_ter /= missing_r)
  83          n = COUNT( temp_tb /= 0. )
  84          if ( n == npixs_rtct .and. m == npixs_rtct) then
  85            call find_std(nlongitude,nlatitude,npixs_rtct, temp_ter, temp_std, missing_r, missing_r)
  86            rtct(i,j) = ( maxval(temp_tb)-tbb(i,j) ) - 3*temp_std*7*0.001
  87            rtct2(i,j) =  3*temp_std*7*0.001
  88          else
  89            rtct(i,j) = 1 !missing_r
  90          end if
  91        end do
  92      end do
  93
  94   end subroutine
  95
  96   subroutine qc_rtmt(nlongitude,nlatitude,tb12,tb13,rtmt)
  97    ! ------------QC: step 2 ----
  98    ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14
  99    ! Relative fourteen minus fifteen test for AHI
 100    ! Relative thirteen minus twelve for AGRI
 101    implicit none
 102     integer, intent(in):: nlongitude,nlatitude
 103    real, intent(in)  :: tb12(nlongitude,nlatitude)
 104    real, intent(in)  :: tb13(nlongitude,nlatitude)
 105    real, intent(out) :: rtmt(nlongitude,nlatitude)
 106    integer, parameter :: npixs_rtmt = 11*11  ! pixles
 107    real :: rtmt_tb12(npixs_rtmt),rtmt_tb13(npixs_rtmt)
 108    real :: max_tb12,max_tb13
 109    integer :: nstar, nend
 110    integer :: i, j, ij, m
 111           RTMT = missing_r
 112      do i=6, nlongitude-5   ! 11*11
 113        do j=6, nlatitude-5   ! 11*11
 114               ij=1
 115           do m=-5, 5
 116               nstar = (ij-1)*11+1
 117               nend  = ij*11
 118               rtmt_tb12(nstar:nend) = tb12(i-5:i+5,j+m)
 119               rtmt_tb13(nstar:nend) = tb13(i-5:i+5,j+m)
 120               ij = ij + 1
 121           end do
 122           call calc_rtmt(npixs_rtmt,rtmt_tb12,rtmt_tb13,max_tb12,max_tb13)
 123           if (max_tb12/=missing_r .and. max_tb13/=missing_r) then
 124              rtmt(i,j)=abs( (tb12(i,j)-tb13(i,j))- &
 125                     (max_tb12-max_tb13) )
 126           else
 127              rtmt(i,j)=missing_r
 128           end if
 129        end do
 130      end do
 131   end subroutine
 132
 133   subroutine qc_cwvt(nlongitude,nlatitude,tb10,tb12,cwvt1)
 134    ! ------------QC: step 3 ----
 135    ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14
 136    ! Cirrus water vapor test for AHI using Ch10, Ch14 5*5pixel
 137    ! Ch10, Ch12 for AGRI
 138    implicit none
 139     integer, intent(in):: nlongitude,nlatitude
 140    real, intent(in) :: tb10(nlongitude,nlatitude)
 141    real, intent(in) :: tb12(nlongitude,nlatitude)
 142    real, intent(out) :: cwvt1(nlongitude,nlatitude)
 143 !   real, intent(out) :: cwvt2(nlongitude,nlatitude)
 144    integer, parameter :: npixs_rtmt = 5*5  ! pixles
 145    real :: cwvt_tb10(npixs_rtmt),cwvt_tb12(npixs_rtmt)
 146    real :: cwvtmp
 147    integer :: nstar, nend
 148    integer :: i, j, ij, m
 149           cwvt1 = missing_r
 150           cwvtmp = 0 !missing_r
 151      do i=3, nlongitude-2   ! 11*11
 152        do j=3, nlatitude-2   ! 11*11
 153               ij=1
 154           do m=-2, 2
 155               nstar = (ij-1)*5+1
 156               nend  = ij*5
 157               cwvt_tb10(nstar:nend) = tb10(i-2:i+2,j+m)
 158               cwvt_tb12(nstar:nend) = tb12(i-2:i+2,j+m)
 159             ij = ij + 1
 160           end do
 161           call calc_correlation (npixs_rtmt,cwvt_tb10,cwvt_tb12,cwvtmp)
 162                   cwvt1(i,j)=cwvtmp
 163 !          call calc_correlation2(npixs_rtmt,cwvt_tb10,cwvt_tb12,cwvt2(i,j))
 164        end do
 165       end do
 166   end subroutine
 167
 168       subroutine qc_tit(nlongitude,nlatitude,tbb,tbb15m,tit,missing_r)
 169    ! ------------QC: step 4 ----
 170    ! 11.2(Zhuge and Zou, JAMC, 2016) for AHI Chan14
 171    ! Cirrus water vapor test for AHI using Ch14, Ch14min before
 172    ! Ch12, Ch12 for AGRI
 173    implicit none
 174     integer, intent(in):: nlongitude,nlatitude
 175    real, intent(in) :: tbb(nlongitude,nlatitude)
 176    real, intent(in) :: tbb15m(nlongitude,nlatitude)
 177    real, intent(in) :: missing_r
 178    real, intent(out) :: tit(nlongitude,nlatitude)
 179    integer :: nstar, nend
 180    integer :: i, j, ij, m
 181
 182    tit=tbb15m-tbb
 183
 184   end subroutine
 185
 186
 187   subroutine find_std(nlongitude,nlatitude,n, arr, std_dev, miss1, miss2)
 188     implicit none
 189     integer, intent(in):: nlongitude,nlatitude
 190     integer, intent(in) :: n
 191     real, intent(in) :: arr(n)
 192     real, intent(in) :: miss1, miss2
 193     real, intent(out):: std_dev
 194     real :: variance, avg
 195     real :: arr_temp(n)
 196     integer :: i, m
 197
 198     avg = 0.
 199     m=0
 200     do i=1, n
 201       if(arr(i)/=miss1 .and. arr(i)/=miss2)then
 202         avg = avg + arr(i)
 203         arr_temp(m+1)=arr(i)
 204         m=m+1
 205       end if
 206     end do
 207
 208     if (m>=6) then
 209       avg=avg/m
 210       variance=0.
 211       do i=1,m
 212         variance=variance+(arr_temp(i)-avg)**2
 213       end do
 214       variance=variance/m
 215       std_dev=sqrt(variance)
 216     else
 217       std_dev=missing_r
 218     end if
 219   end subroutine
 220
 221   subroutine calc_rtmt(npixs_rtmt,temp_tb12,temp_tb13,max_tb12,max_tb13)
 222     implicit none
 223     integer, intent(in) :: npixs_rtmt
 224     real, intent(in) :: temp_tb12(npixs_rtmt),temp_tb13(npixs_rtmt)
 225     real, intent(out):: max_tb12, max_tb13
 226     real :: variance, avg
 227     real :: arr_temp12(npixs_rtmt),arr_temp13(npixs_rtmt)
 228     integer :: i, m
 229
 230     m=0
 231     do i=1, npixs_rtmt
 232       if( temp_tb12(i) /=missing_r .and. temp_tb12(i) /=missing_r .and. &
 233           temp_tb13(i) /=missing_r .and. temp_tb13(i) /=missing_r )then
 234           arr_temp12(m+1)=temp_tb12(i)
 235           arr_temp13(m+1)=temp_tb13(i)
 236         m=m+1
 237       end if
 238     end do
 239     ! 80 per 121 would be ok
 240     if (m>=80) then
 241       max_tb12=maxval(arr_temp12(1:m))
 242       max_tb13=maxval(arr_temp13(1:m))
 243     else
 244       max_tb12=missing_r
 245       max_tb13=missing_r
 246     end if
 247   end subroutine
 248
 249   subroutine calc_correlation( n, arr1, arr2, r)
 250     implicit none
 251
 252     integer, intent(in) :: n
 253     real, intent(in)    :: arr1(n)
 254     real, intent(in)    :: arr2(n)
 255     real, intent(out) :: r
 256 ! below variables can be output: m, t(significant)
 257     integer, parameter :: lag = 0
 258     real    :: t
 259     integer :: m
 260 ! temp variables
 261     real    :: arr1_tem(n), arr2_tem(n)
 262 !    real, allocatable :: x(:)
 263 !    real, allocatable :: y(:)
 264 !    real, allocatable :: xdev(:)
 265 !    real, allocatable :: ydev(:)
 266 !    real, allocatable :: xdevydev(:)
 267 !    real, allocatable :: xdevxdev(:)
 268 !    real, allocatable :: ydevydev(:)
 269     real :: x(3*n)
 270     real :: y(3*n)
 271     real :: xdev(3*n)
 272     real :: ydev(3*n)
 273     real :: xdevydev(3*n)
 274     real :: xdevxdev(3*n)
 275     real :: ydevydev(3*n)
 276     real :: xmn
 277     real :: ymn
 278     real :: COVxy
 279     real :: Sx
 280     real :: Sy
 281     integer :: i
 282
 283     arr1_tem=arr1
 284     arr2_tem=arr2
 285 !    allocate( x(1:3*N) )
 286 !    allocate( y(1:3*N) )
 287
 288     x(:) = missing_r
 289     y(:) = missing_r
 290     x(N+1:2*N) = arr1_tem(:)
 291     y(N+1:2*N) = arr2_tem(:)
 292
 293     y = EOSHIFT( y, SHIFT = -lag, BOUNDARY = missing_r )
 294
 295     WHERE ( x == missing_r ) y = missing_r
 296     WHERE ( y == missing_r ) x = missing_r
 297
 298     m = COUNT( x /= missing_r )
 299     r = missing_r
 300   if (m >= 16) then
 301
 302     xmn = sum( x, dim = 1, mask = x /= missing_r ) / real(m)
 303     ymn = sum( y, dim = 1, mask = y /= missing_r ) / real(m)
 304 !    allocate ( xdev(1:3*N) )
 305 !    allocate ( ydev(1:3*N) )
 306     xdev(:) = x(:) - xmn
 307     where ( x == missing_r ) xdev(:) = missing_r
 308     ydev(:) = y(:) - ymn
 309     where ( y == missing_r ) ydev(:) = missing_r
 310 !    allocate ( xdevydev(1:3*N) )
 311     xdevydev(:) = xdev(:) * ydev(:)
 312     where ( x == missing_r ) xdevydev(:) = missing_r
 313
 314     COVxy = sum( xdevydev, dim = 1, mask = xdevydev /= missing_r ) / real(m)
 315
 316 !    allocate ( xdevxdev(1:3*N) )
 317     xdevxdev(:) = xdev(:) * xdev(:)
 318     where ( x == missing_r ) xdevxdev(:) = missing_r
 319     Sx = sqrt( sum( xdevxdev, dim = 1, mask = xdevxdev /= missing_r ) / real(m) )
 320
 321 !    allocate ( ydevydev(1:3*N) )
 322     ydevydev(:) = ydev(:) * ydev(:)
 323     where ( y == missing_r ) ydevydev(:) = missing_r
 324     Sy = SQRT( SUM( ydevydev, DIM = 1, MASK = ydevydev /= missing_r ) / REAL(m) )
 325     if ( Sx /= missing_r .and. Sy/= missing_r .and. (Sx * Sy/=0) ) r = COVxy / ( Sx * Sy )
 326   else
 327     r = missing_r
 328   end if
 329 !   deallocate( x, y )
 330
 331   end subroutine
 332
 333
 334
 335
 336
 337
 338 end module