var/da/da_radiance/da_cloud_detect.inc

   1 subroutine da_cloud_detect(isensor,nchannels,ndim,kts,kte,n,iv)
   2
   3 !** *CLOUD_DETECT* - CLOUD FLAGGING FOR IR Channels
   4
   5 !    AUTHOR: THOMAS AULIGNE      DATE : 01/08/2005
   6 !
   7 !    PURPOSE.
   8 !    -------
   9 !    FLAG THE PRESENCE OF CLOUD CONTAMINATION IN IR CHANNELS
  10 !
  11 !**  INTERFACE.
  12 !    ---------
  13 !    WHERE nchannels    : Number of channels
  14 !          kts          : model level corresponding to 100hPa (top of initial cloud search)
  15 !          kte          : model level corresponding to surface (lower extent of cloud)
  16 !          rad_obs      : Potentially cloudy observations
  17 !          rad_clr      : Clear radiance from Model
  18 !          rad_ovc      : Model overcast radiance estimates
  19 !          cloud_flag   : Cloud flag by channel; 1=clear, -1=cloudy
  20 !
  21 !**  EXTERNALS
  22 !    ---------
  23 !    N2QN1  - Minimization algorithm (double-precision constrained version of M1QN3)
  24 !**  METHODS
  25 !    CLOUD DETECTION SCHEME MMR FROM AULIGNÉ T. MWR (2014).OR. PF FROM XU ET AL., GMD (2016)
  26 !    MODIFICATIONS
  27 !    -------------
  28 !    BY DONGMEI XU 201904
  29 !    PURPOSE, ADD CLOUD DETECTION METHOD BASED ON PARTICLE FILTER
  30 !    METHOD, Xu et al., 2016: A method for retrieving clouds with satellite infrared radiances using the particle filter. Geosci. Model Dev., 9, 3919–3932
  31 !
  32 !**  -----------------------------------------
  33
  34 IMPLICIT NONE
  35
  36 !* 0.1 Global arrays
  37
  38 INTEGER,INTENT(IN)    :: isensor             ! sensor index.
  39 INTEGER,INTENT(IN)    :: nchannels           ! number of channels
  40 INTEGER,INTENT(IN)    :: ndim                ! model levels between surface (lower extent of cloud) and 100hPa (top of cloud search)
  41 INTEGER,INTENT(IN)    :: kts                 ! model level corresponding to 100hPa (top of initial cloud search)
  42 INTEGER,INTENT(IN)    :: kte                 ! model level corresponding to surface (lower extent of cloud)
  43 INTEGER,INTENT(IN)    :: n                   ! pixel index
  44 type (iv_type), intent(inout)  :: iv         ! O-B structure.
  45
  46 INTEGER,PARAMETER     :: NITER             = 100
  47 INTEGER,PARAMETER     :: NBAND             = 1
  48 LOGICAL,PARAMETER     :: LPRECON           = .false.
  49 INTEGER,PARAMETER     :: NEIGNVEC          = 4
  50 INTEGER,PARAMETER     :: AIRS_Max_Channels = 2378
  51 INTEGER,PARAMETER     :: IASI_Max_Channels = 8079
  52 !! local declarations
  53 INTEGER               :: ichan(nchannels)    ! AIRS and IASI channel IDs
  54 REAL                  :: rad_obs(nchannels)  ! Observed radiance
  55 REAL                  :: rad_clr(nchannels)  ! Model clear radiance estimates
  56 REAL                  :: rad_ovc(nchannels,ndim-1) ! RT overcast radiance estimates
  57 double precision      :: px(ndim) !neignvec)        ! Cloud fractions
  58 REAL                  :: rad_cld(nchannels)
  59 INTEGER               :: ich,ilev,jlev,i,j,JBAND
  60 double precision      :: ZF, ZF_CLR
  61 double precision      :: ZG(ndim)
  62 double precision      :: binf(ndim), bsup(ndim)
  63 REAL                  :: AMAT(nchannels,ndim)
  64 INTEGER               :: NCHAN,k
  65 LOGICAL               :: LMATCH
  66 INTEGER               :: Band_Size(5)
  67 INTEGER               :: Bands(IASI_Max_Channels,5)
  68 integer               :: cldtoplevel
  69
  70 ! Hessian evaluation
  71  REAL                 :: hessian(ndim,ndim), eignvec(ndim,ndim), eignval(ndim)
  72
  73 !! local declarations for N2QN1 !!
  74 INTEGER               :: NRZ, impres, io, IMODE, NSIM, nit, izs(2)
  75 double precision      :: ZDF1, ZDXMIN, ZEPSG
  76 double precision ,ALLOCATABLE :: ZRZ(:)
  77 real, allocatable     :: RZS(:)
  78 INTEGER, ALLOCATABLE  :: IZ(:)
  79 DOUBLE PRECISION, ALLOCATABLE :: DZS(:)
  80 INTEGER               :: gn
  81 REAL :: ZHOOK_HANDLE
  82
  83
  84 logical   :: iasi,airs, modis,imager,sounder,cris,giirs,ahi
  85 double precision, allocatable            :: ppx(:,:),wx(:),jo(:)
  86 integer                   :: p1,ii,jj
  87 double precision      :: tmp
  88 ! Initializations
  89
  90
  91 iasi  = trim(rttov_inst_name(rtminit_sensor(isensor))) == 'iasi'
  92 airs  = trim(rttov_inst_name(rtminit_sensor(isensor))) == 'airs'
  93 imager  = trim(rttov_inst_name(rtminit_sensor(isensor))) == 'imager'
  94 ahi  = trim(rttov_inst_name(rtminit_sensor(isensor))) == 'ahi'
  95 Band_Size(:)   = 0
  96 Bands(:,:)     = 0
  97
  98 if ( iasi ) then
  99       Band_Size(1:5) = (/ 193, 15, 116, 4, 15 /)
 100
 101       Bands(1:Band_Size(1),1) = &
 102 &      (/    16,   38,   49,   51,   55,   57,   59,   61,   63,   66, &
 103 &            70,   72,   74,   79,   81,   83,   85,   87,   89,   92, &
 104 &            95,   97,   99,  101,  104,  106,  109,  111,  113,  116, &
 105 &           119,  122,  125,  128,  131,  133,  135,  138,  141,  144, &
 106 &           146,  148,  151,  154,  157,  159,  161,  163,  165,  167, &
 107 &           170,  173,  176,  178,  179,  180,  183,  185,  187,  189, &
 108 &           191,  193,  195,  197,  199,  201,  203,  205,  207,  210, &
 109 &           212,  214,  217,  219,  222,  224,  226,  228,  230,  232, &
 110 &           234,  236,  239,  241,  242,  243,  246,  249,  252,  254, &
 111 &           256,  258,  260,  262,  265,  267,  269,  271,  272,  273, &
 112 &           275,  278,  280,  282,  284,  286,  288,  290,  292,  294, &
 113 &           296,  299,  301,  303,  306,  308,  310,  312,  314,  316, &
 114 &           318,  320,  323,  325,  327,  329,  331,  333,  335,  337, &
 115 &           339,  341,  343,  345,  347,  350,  352,  354,  356,  358, &
 116 &           360,  362,  364,  366,  369,  371,  373,  375,  377,  379, &
 117 &           381,  383,  386,  389,  398,  401,  404,  407,  410,  414, &
 118 &           416,  426,  428,  432,  434,  439,  445,  457,  515,  546, &
 119 &           552,  559,  566,  571,  573,  646,  662,  668,  756,  867, &
 120 &           921, 1027, 1090, 1133, 1191, 1194, 1271, 1805, 1884, 1946, &
 121 &          1991, 2094, 2239 /)
 122
 123
 124       Bands(1:Band_Size(2),2) = &
 125 &      (/ 1479, 1509, 1513, 1521, 1536, 1574, 1579, 1585, 1587, 1626, &
 126 &         1639, 1643, 1652, 1658, 1671  /)
 127
 128       Bands(1:Band_Size(3),3) = &
 129 &      (/ 2119, 2213, 2271, 2321, 2398, 2701, 2741, 2819, 2889, 2907, 2910, &
 130 &         2919, 2939, 2944, 2948, 2951, 2958, 2977, 2985, 2988, 2991, &
 131 &         2993, 3002, 3008, 3014, 3027, 3029, 3036, 3047, 3049, 3053, &
 132 &         3058, 3064, 3069, 3087, 3093, 3098, 3105, 3107, 3110, 3127, &
 133 &         3136, 3151, 3160, 3165, 3168, 3175, 3178, 3207, 3228, 3244, &
 134 &         3248, 3252, 3256, 3263, 3281, 3303, 3309, 3312, 3322, 3375, &
 135 &         3378, 3411, 3438, 3440, 3442, 3444, 3446, 3448, 3450, 3452, &
 136 &         3454, 3458, 3467, 3476, 3484, 3491, 3497, 3499, 3504, 3506, &
 137 &         3509, 3518, 3527, 3555, 3575, 3577, 3580, 3582, 3586, 3589, &
 138 &         3599, 3653, 3658, 3661, 4032, 5368, 5371, 5379, 5381, 5383, &
 139 &         5397, 5399, 5401, 5403, 5405, 5455, 5480, 5483, 5485, 5492, &
 140 &         5502, 5507, 5509, 5517, 5558  /)                                           !&    1812, 1826, 1843  /)
 141
 142       Bands(1:Band_Size(4),4) = &
 143 &      (/   5988, 5992, 5994, 6003  /)                              !&    1921, 1923, 1924, 1928, 1937  /)
 144
 145       Bands(1:Band_Size(5),5) = &
 146 &      (/  6982, 6985, 6987, 6989, 6991, 6993, 6995, 6997, 7267, 7269, &
 147 &          7424, 7426, 7428, 7885, 8007 /)
 148
 149 else if (airs) then
 150       Band_Size(1:5) = (/86, 0, 0, 16, 0 /)
 151
 152       Bands(1:Band_Size(1),1) = &
 153 &    (/                                                 &              !&      1,   6,   7,  10,  11,  15,  16,  17,  20,  21, &
 154 &                                                       &              !&     22,  24,  27,  28,  30,  36,  39,  40,  42,  51, &
 155 &                                                       &              !&     52,  54,  55,  56,  59,  62,  63,  68,  69,  71, &
 156 &                                                       &              !&     72,  73,  74,  75,  76,  77,  78,  79,  80,  82, &
 157 &                     92,  93,  98,  99, 101, 104, 105, &              !&     83,  84,  86,  92,  93,  98,  99, 101, 104, 105, &
 158 &     108, 110, 111, 113, 116, 117, 123, 124, 128, 129, &
 159 &     138, 139, 144, 145, 150, 151, 156, 157, 159, 162, &
 160 &     165, 168, 169, 170, 172, 173, 174, 175, 177, 179, &
 161 &     180, 182, 185, 186, 190, 192,      198, 201, 204, &              !&     180, 182, 185, 186, 190, 192, 193, 198, 201, 204, &
 162 &     207, 210,      215, 216,      221,      226, 227, &              !&     207, 210, 213, 215, 216, 218, 221, 224, 226, 227, &
 163 &     232,                     252, 253, 256, 257, 261, &              !&     232, 239, 248, 250, 251, 252, 253, 256, 257, 261, &
 164 &     262, 267, 272, 295, 299,      305,           310, &              !&     262, 267, 272, 295, 299, 300, 305, 308, 309, 310, &
 165 &          321, 325, 333, 338, 355, 362, 375, 453, 475, &              !&     318, 321, 325, 333, 338, 355, 362, 375, 453, 475, &
 166 &     484, 497, 528, 587, 672, 787, 791, 843, 870, 914, &
 167 &     950 /)
 168
 169       Bands(1:Band_Size(4),4) = &
 170 &    (/ 1852, 1865, 1866,       1868, 1869, 1872, 1873,       1876, &  !&    1852, 1865, 1866, 1867, 1868, 1869, 1872, 1873, 1875, 1876,
 171 &             1881, 1882, 1883,                   1911, 1917, 1918, &  !&    1877, 1881, 1882, 1883, 1884, 1897, 1901, 1911, 1917, 1918, &
 172 &                   1924, 1928        /)                               !&    1921, 1923, 1924, 1928, 1937  /)
 173
 174 else if (imager) then
 175       Band_Size(1) = 2
 176       Bands(1:Band_Size(1),1) = &
 177 &   (/4,5/)
 178 else if (ahi) then
 179       Band_Size(1) = 4
 180       Bands(1:Band_Size(1),1) = &
 181 &   (/7,8,9,10/)
 182 end if
 183
 184     allocate(ppx(ndim*11,ndim))
 185     allocate(wx(ndim*11))
 186     allocate(jo(ndim*11))
 187     wx=1.0
 188     jo=0.0
 189     px(1:ndim-1) = 1.0/ndim
 190     px(ndim)     = 1.0 - SUM(px(1:ndim-1))
 191     ichan        = iv%instid(isensor)%ichan(1:nchannels)
 192     rad_clr      = iv%instid(isensor)%rad_xb(1:nchannels,n)              !iv%instid(isensor)%tb_xb(1:nchan,n)
 193     rad_obs      = iv%instid(isensor)%rad_obs(1:nchannels,n)             !iv%instid(isensor)%tb_inv(1:nchan,n) + rad_clr
 194     rad_ovc      = iv%instid(isensor)%rad_ovc(1:nchannels,kts+1:kte,n)
 195
 196     nchan        = 0
 197     AMAT(:,:)    = 0.0
 198     px(1:ndim-1) = 0.0
 199     px(ndim)     = 1.0
 200     ZF_CLR       = 0.0
 201     nit          = niter
 202
 203 !do ich=1,nchannels
 204 !   CALL CRTM_Planck_Radiance(11,ichan(ich),tb_obs(ich),rad_obs(ich))
 205 !   CALL CRTM_Planck_Radiance(11,ichan(ich),tb_clr(ich),rad_clr(ich))
 206 !end do
 207
 208     !--------------------!
 209     !   Loop over band   !
 210     !--------------------!
 211     BAND_LOOP: DO JBAND = 1, NBAND
 212       DO i = 1, Band_Size(JBAND)
 213         LMATCH = .FALSE.
 214         DO ich=1,nchannels
 215           IF (ichan(ich)/= Bands(i,JBAND)) CYCLE
 216           IF ((rad_obs(ich)<=0.0).OR.(rad_obs(ich)>1000.0)) CYCLE
 217           IF ((rad_clr(ich)<=0.0).OR.(rad_clr(ich)>1000.0)) CYCLE
 218           IF (ANY(rad_ovc(ich,1:NDIM-1)<=0.0)) CYCLE
 219           IF (ANY(rad_ovc(ich,1:NDIM-1)>1000.0)) CYCLE
 220
 221           LMATCH               = .TRUE.    !! Found match for channel
 222           nchan                = nchan +1
 223           AMAT(nchan,1:ndim-1) = rad_ovc(ich,1:NDIM-1) / rad_obs(ich)
 224           AMAT(nchan,ndim)     = rad_clr(ich) / rad_obs(ich)
 225           ZF_CLR               = ZF_CLR + 0.5*(AMAT(nchan,ndim)-1.0)**2
 226
 227         if (use_clddet==2)  then
 228           p1=0
 229           do ii=0,10
 230             do jj=1,ndim-1
 231               p1=p1+1
 232               ppx(p1,1:ndim)=0 !initialization
 233               ppx(p1,jj)=real(ii)/10
 234               ppx(p1,ndim)=1-ppx(p1,jj) !initialization
 235               jo(p1)=iv%instid(isensor)%tb_error(ich,n)*(1-SUM(ppx(p1,1:ndim)*amat(nchan,1:ndim)))**2
 236            end do
 237           end do
 238           ! step 2 calculate the weight
 239           do ii=1,p1
 240              ! jo(ii)=jo(ii)-mv
 241               tmp=exp(-jo(ii))/sum(exp(-jo(1:p1))) ! normalize the weight
 242               wx(ii)=tmp*wx(ii) !1/jo  !exp(-jo)
 243           end do
 244           wx(:)=wx(:)/sum(wx(1:p1))
 245           px=0
 246           do k = 1,ndim
 247             px(k)= SUM(wx(1:p1)*ppx(1:p1,k))
 248           end do
 249           px(1:ndim) = px(1:ndim) / SUM(px(1:ndim))   ! Re-normalization
 250         end if
 251         ENDDO
 252         IF (.NOT. LMATCH) then
 253            if (print_detail_rad) then
 254               write(unit=message(1),fmt='(A,2I8)') 'CLOUD_DETECT: No match for channel:',i,Bands(i,JBAND)
 255               call da_message(message(1:1))
 256            endif
 257         ENDIF
 258       ENDDO
 259     ENDDO BAND_LOOP                      ! Loop over band
 260
 261 if (use_clddet==1) then
 262     !--------------------!
 263     ! Hessian evaluation !
 264     !--------------------!
 265     IF (LPRECON) THEN
 266       hessian(:,:)= 0.0
 267       DO ilev=1, NDIM
 268         DO jlev=ilev, NDIM
 269           DO J=1,NCHAN
 270               hessian(ilev,jlev) = hessian(ilev,jlev)  + &
 271                                   (AMAT(J,ilev)-AMAT(J,NDIM)) * &
 272                                   (AMAT(J,jlev)-AMAT(J,NDIM))
 273           ENDDO
 274           hessian(jlev,ilev) = hessian(ilev,jlev)
 275         ENDDO
 276       ENDDO
 277     ENDIF
 278
 279      !-----------------!
 280      ! n2qn1 minimizer !
 281      !-----------------!
 282       impres = 2
 283       io     = 66
 284       NSIM   = NITER+5
 285       ZDXMIN = 1.e-6
 286       ZEPSG  = 1.e-3 !e-9
 287       IMODE  = 1
 288       NRZ    = NDIM*(NDIM+9)/2 ! N2QN1
 289       ALLOCATE(IZ(2*NDIM +1))
 290       ALLOCATE(ZRZ(NRZ))
 291       ALLOCATE(DZS(NCHAN*NDIM))
 292       allocate(rzs(ndim*neignvec))
 293       binf   = -1000.0
 294       bsup   = 1000.0
 295       izs(1) = nchan
 296       izs(2) = neignvec
 297       rzs    = 0.0
 298       ZRZ    = 0.0
 299       dzs(1:NCHAN*NDIM)=RESHAPE(AMAT(1:NCHAN,1:NDIM),(/NCHAN*NDIM/))
 300
 301       IF (LPRECON) THEN
 302         IMODE = 2
 303         i     = 0
 304         DO ilev=1, NDIM
 305           DO jlev=ilev, NDIM
 306             i = i + 1
 307             ZRZ(i) = hessian(jlev,ilev)
 308           ENDDO
 309         ENDDO
 310       ENDIF
 311 !      rzs(1:ndim*neignvec)                  = RESHAPE(eignvec(1:ndim,1:neignvec),(/ndim*neignvec/))
 312 !      rzs(ndim*neignvec+1:(ndim+1)*neignvec)= eignval(1:neignvec)
 313
 314       call da_cloud_sim(0,NDIM,px,ZF,ZG,izs,RZS,DZS)
 315       ZDF1      = 1.e-1*ZF
 316
 317
 318 !      call da_error(__FILE__,__LINE__, &
 319 !             (/"inria_n2qn1 is not implemented here, please contact the author of this subroutine."/))
 320 !    call inria_n2qn1(da_cloud_sim,NDIM,px,ZF,ZG,(/(ZDXMIN,jlev=1,NDIM)/),ZDF1, &
 321 !                ZEPSG,impres,io,IMODE,nit,NSIM,binf,bsup,IZ,ZRZ,izs,RZS,DZS)
 322
 323       if (allocated(iz))  deallocate(iz)
 324       if (allocated(zrz)) deallocate(zrz)
 325       if (allocated(dzs)) deallocate(dzs)
 326       if (allocated(rzs)) deallocate(rzs)
 327 end if !mmr
 328       deallocate(ppx)
 329       deallocate(wx)
 330       deallocate(jo)
 331       !-----------------!
 332       ! Cloudy radiance !
 333       !-----------------!
 334       DO ich=1,nchannels
 335         rad_cld(ich) = SUM(px(1:ndim-1) * rad_ovc(ich,1:ndim-1)) + px(ndim) * rad_clr(ich)
 336
 337         if (ABS(rad_cld(ich)-rad_clr(ich)) < 0.01*rad_clr(ich)) then
 338            iv%instid(isensor)%cloud_flag(ich,n) = qc_good
 339         else
 340            iv%instid(isensor)%cloud_flag(ich,n) = qc_bad
 341         end if
 342       ENDDO
 343
 344     ! Dump cloud top pressure
 345     do ilev = kte, kts+2, -1
 346       if (px(ilev-kts+1) > 0.01) cldtoplevel = ilev
 347     end do
 348
 349     if (rtm_option == rtm_option_rttov) then
 350 #ifdef RTTOV
 351        iv%instid(isensor)%clwp(n) = coefs(isensor)%coef%ref_prfl_p(cldtoplevel)
 352 #endif
 353     elseif (rtm_option == rtm_option_crtm) then
 354 #ifdef CRTM
 355        iv%instid(isensor)%clwp(n) = iv%instid(isensor)%pm(cldtoplevel,n)
 356 #endif
 357     end if
 358
 359 end subroutine da_cloud_detect