var/da/da_radar/da_radzicevar.inc

   1 subroutine da_radzicevar(qvp0,qra0,qsn0,qgr0,qnr0,qns0,qng0,tmk0,prs,dbz,                 &
   2                          in0r,in0s,in0g,rn0_r,rn0_s,rn0_g,                                &
   3                          rhos,rhog,dtmk,dqvp,dqra,dqsn,dqgr,dqnr,dqns,dqng,zmm,tlopt,     &
   4                          gropt,zmm_ref)
   5
   6 !===================================================================
   7 !      Following Jung et al., 2008
   8 !===================================================================
   9
  10   implicit none
  11
  12   integer,intent(in) :: in0r,in0s,in0g
  13   integer,intent(in) :: tlopt,gropt  ! tlopt=0 nonlinear,=1 tl linear; gropt=0 tl linear; >=1 adj
  14
  15   integer,parameter  :: dualpol_opt=0 ! completed but not yet been tested
  16
  17   real           ::  rn0_r,rn0_s,rn0_g    ! intercept parameter 1/m^4
  18   real           ::  qra0,qsn0,qgr0, &    ! these are read in background states
  19                      qnr0,qns0,qng0, &
  20                      tmk0,qvp0
  21
  22   real           :: qra,qvp,        &    ! clone the background
  23                     qsn,qgr,tmk,    &
  24                     prs,dbz,qnr,    &
  25                     qns,qng,ref,    &
  26                     zrs,zss,zhs,    &
  27                     zmm,zmm_ref,    &
  28  ! above: background, below: increment for tlopt==1 and gropt==0 but adj for tlopt==1 and gropt>=1
  29                     dqra,dqsn,dqgr, &
  30                     dqnr,dqns,dqng, &
  31                     dtmk,dqvp
  32
  33   real           :: rgas=287.04
  34   real           :: z_e
  35
  36   integer        :: i,j,k,ii,jj,kk
  37
  38 ! reflectivity operator constant
  39   real,parameter :: rdrwave = 107.0  ! unit mm S band
  40   real,parameter :: lambda = rdrwave
  41   real,parameter :: Kw2 = 0.93
  42   real           :: pi = 3.1415926
  43   real,parameter :: mm3todBZ=1.0E+9
  44   real           :: rhor = 1000    ! kg m^-3 rainwater density
  45   real           :: rhos           ! kg m^-3 snow density
  46   real           :: rhog           ! kg m^-3 graupel density
  47   real           :: rhoair         ! air density
  48
  49   real :: rhows,rhowg  !  wet snow, wet graupel
  50
  51 ! temporal mixing ratio
  52   real           :: prain  ! pure rain mixing ratio
  53   real           :: dsnow  ! dry snow mixing ratio
  54   real           :: dgr    ! dry graupel mixing ration
  55   real           :: wsnow  ! wet snow mixing ratio
  56   real           :: wgr    ! wet graupel mixing ratio
  57
  58 !     parameters for rain
  59   real           :: alpha_ra = 4.28e-4
  60   real           :: alpha_rb = 4.28e-4
  61   real           :: beta_ra = 3.04
  62   real           :: beta_rb = 2.77
  63   real           :: alphar = 0
  64   real           :: dr = 3
  65   real           :: cr = 3.1415926/6
  66   real           :: zrh,zrv
  67
  68   real           :: para10r  ! =(1.-(2.*betarx+1.0)/4.)
  69   real           :: para14r  ! mm3todBZ*(4*lambda**4*alpharx**2/(pi**4*Kw2))
  70                              ! *(pi*rhor/rhoa)**(p2/4.)**(n0r)**para10r
  71                              ! *gamma(-(2.*betarx+1.0)*1d0)
  72                              ! for rainwater
  73   real           :: ronv,gonv,sonv                ! intercept parameter
  74
  75 ! parameters for snow and graupel/hail
  76
  77   integer,parameter :: npara_alpharxa=7
  78   real           :: para_alpha_rxa(npara_alpharxa,2) ! second dimension: 1 for snow
  79                                                      ! 2 for hail/graupel
  80   real,save      :: para_alpha_rxb(npara_alpharxa,2)       ! precalculated coefficients in contribution equations of ice species
  81   real           :: phimean=0,sigma,ice_abc(3),ice_bac(3)  ! A, B, and C in contribution equations of ice species
  82   real           :: fw                                  ! water fraction
  83   real           :: pxabk_all(3)                        !
  84   real           :: pxkh,pxkv                           ! the sum of fwx term
  85   real           :: pxkh_tlr,pxkv_tlr,pxkh_tlx,pxkv_tlx ! particle derivative of pxkh/pxkv for rain (tlr) and ice species (tlx)
  86   real           :: zsh,zsv,zgh,zgv,zdsh,zdsv,zdgh,zdgv ! contribution from snow/graupel
  87                                                         ! (thelast character means horizontal/vertical)
  88
  89   real           :: para1sg   ! para1sg=mm3todBZ*gamma(7)*lambda**4/pi**4/Kw2*(pi*rhox/rhoa)**1.75
  90
  91   real           :: cs=3.1415926/6
  92   real           :: cice=440.
  93   real           :: ds=3.
  94   real           :: dg=3.
  95   real           :: alphas = 0
  96   real           :: alphag = 0
  97
  98   real           :: alpha_rdsa=0.194*10.**(-4)
  99   real           :: alpha_rdsb=0.191*10.**(-4)
 100   !---------------------------------------------
 101   !real           :: alpha_rdha=0.191*10.**(-3)  ! hail
 102   !real           :: alpha_rdhb=0.165*10.**(-3)  ! hail
 103   !--------------------------------------------
 104   real           :: alpha_rdha=0.105*10.**(-3)  ! graupel
 105   real           :: alpha_rdhb=0.092*10.**(-3)  ! graupel
 106
 107   real           :: zh_e,zv_e
 108
 109   ! for dry graupel
 110   REAL,PARAMETER :: sigmahd = 1.0472
 111   REAL,PARAMETER :: Ahd = 0.4308
 112   REAL,PARAMETER :: Bhd = 0.3192
 113   REAL,PARAMETER :: Chd = 0.1250
 114   REAL,PARAMETER :: Dhd = 0.3750
 115   REAL,PARAMETER :: Ckhd = 0.1116
 116   ! for dry snow
 117   REAL,PARAMETER :: sigmas = 0.3491
 118   REAL,PARAMETER :: Asd = 0.8140
 119   REAL,PARAMETER :: Bsd = 0.0303
 120   REAL,PARAMETER :: Csd = 0.0778
 121   REAL,PARAMETER :: Dsd = 0.4221
 122   REAL,PARAMETER :: Cksd = 0.7837
 123
 124   real           :: ice_abc_d(3),ice_bac_d(3)
 125   real           :: pdfrrs,pdfrrg,pdfsrs,pdfgrg    ! partial derivative of F (r&s,r&g) with respect to (r,s,g)
 126   real           :: pdfrhot,pdfrhoq                ! partial derivative of rho for t and qv
 127
 128   real           :: prain_coef,dsnow_coef,dgr_coef ! coefficient to determine the ratio of
 129                                                    !pure rain, dry swno/graupel in total qr and qs/qg, respectively.
 130   real           :: qthres
 131
 132   real           :: temr01,temr02,temr03    ! temporary variable for real
 133   integer        :: temi01,temi02,temi03    ! temporary variable for integer
 134
 135   real           :: para_pr,para_pdx_dq,para_pdx_df        ! for
 136
 137   !---------------------------------------------------------------------------------
 138
 139   data para_alpha_rxa(:,1)/0.194e-4,7.094e-4,2.135e-4,-5.225e-4,0,0,0/
 140   data para_alpha_rxb(:,1)/0.191e-4,6.916e-4,-2.841e-4,-1.160e-4,0,0,0/
 141   !----hail--------------
 142   !data para_alpha_rxa(:,2)/0.191e-3,2.39e-3,-12.57e-3,38.71e-3,-65.53e-3,56.16e-3,  &
 143   !                        -18.98e-3/
 144   !data para_alpha_rxb(:,2)/0.165e-3,1.72e-3,-9.92e-3,32.15e-3,-56.0e-3,48.83e-3,   &
 145   !                        -16.69e-3/
 146   !----graupel-----------
 147   data para_alpha_rxa(:,2)/1.05E-04,1.82E-03,-3.77E-03,-7.97E-04,1.63E-02,-2.20E-02, &
 148                            8.74E-03/
 149   data para_alpha_rxb(:,2)/9.25E-05,1.93E-03,-9.79E-03,2.92E-02,-4.82E-02,3.93E-02, &
 150                           -1.22E-02/
 151
 152
 153   qthres=rf_qthres
 154
 155   qra=qra0
 156   qsn=qsn0
 157   qgr=qgr0
 158   tmk=tmk0
 159   qvp=qvp0
 160 !  if(qra0<qthres) qra=qthres
 161 !  if(qsn0<qthres) qsn=qthres
 162 !  if(qgr0<qthres) qgr=qthres
 163
 164   if(tlopt>=1.and.gropt>=1) then
 165     dqra=0
 166     dqgr=0
 167     dqsn=0
 168     dqvp=0
 169     dtmk=0
 170     dqnr=0
 171     dqns=0
 172     dqng=0
 173   endif
 174
 175   zrh=0
 176   zsh=0
 177   zgh=0
 178   zdsh=0
 179   zdgh=0
 180   zrv=0
 181   zsv=0
 182   zgv=0
 183   zdsv=0
 184   zdgv=0
 185
 186 !  if(qra<qthres.and.qsn<qthres.and.qgr<qthres) then
 187 !     zrs=0
 188 !     zss=0
 189 !     zhs=0
 190 !     ref=0
 191 !     zmm=0
 192 !     dbz=0
 193 !     return
 194 !  endif
 195
 196   call da_radzicevar_rhoair_tl(tlopt,rhoair,prs,rgas,tmk,qvp,pdfrhot,pdfrhoq)
 197
 198 !-------------------------------------------------------
 199 ! Calculate variable intercept parameters if wanted
 200 !-------------------------------------------------------
 201   call da_radzicevar_prepare_interceptpara(in0s,in0g,in0r,rn0_s,rn0_g,rn0_r,sonv,gonv,ronv)
 202 !-------------------------------------------------------
 203 ! Calculate mixing ratios (pure rain, dry snow/graupel, wet snow/graupel)
 204 !-------------------------------------------------------
 205   call da_radzicevar_prepare_mixingratios(tlopt,prain_coef,dsnow_coef,dgr_coef, &
 206                             prain,dsnow,wsnow,dgr,wgr,            &
 207                             qra,qgr,qsn,qthres,                   &
 208                             pdfrrs,pdfrrg,pdfsrs,pdfgrg           &
 209                             )
 210
 211   call da_radzicevar_prepare_zmm_adj(tlopt,gropt,zmm,zmm_ref)
 212 ! ==================FOR RAIN=============================
 213   call da_radzicevar_parameter_zrx(para10r,para14r,para_pr,rhoair,rhor,prain,   &
 214                      beta_ra,alpha_ra,mm3todBZ,lambda,Kw2,pi,     &
 215                      ronv)
 216
 217   if(in0r.eq.1) then
 218     ! for two moment microphysics scheme, not yet completed
 219   else
 220     if(tlopt==0) zrh=para14r*(rhoair*prain)**(1.-para10r)
 221
 222   endif
 223
 224   if(dualpol_opt==1) then
 225     call da_radzicevar_parameter_zrx(para10r,para14r,para_pr,rhoair,rhor,prain,   &
 226                        beta_rb,alpha_rb,mm3todBZ,lambda,Kw2,pi,     &
 227                        ronv)
 228
 229     if(in0r.eq.1) then
 230      ! for two moment microphysics scheme, not yet completed
 231     else
 232       if(tlopt==0) zrv=para14r*(rhoair*prain)**(1.-para10r)
 233
 234     endif
 235   endif
 236   zrs=zrh
 237 ! ==================FOR snow=============================
 238   if(rf_noice==0) then
 239 ! -------------------for wet snow----------------------
 240
 241     call da_radzicevar_waterfraction(qra,qsn,fw)
 242     rhows=rhos*(1.-fw**2)+rhor*fw**2
 243
 244     call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhows,wsnow,cs,ds,alphas,   &
 245                        mm3todBZ,lambda,Kw2,pi,sonv)
 246
 247     call da_radzicevar_sigma_in_abc(qsn,fw,1,sigma)  ! for snow
 248     call da_radzicevar_calc_ice_abc(phimean,sigma,ice_abc)
 249
 250     ice_bac(1)=ice_abc(2)
 251     ice_bac(2)=ice_abc(1)
 252     ice_bac(3)=ice_abc(3)
 253     pxkh=0
 254     pxkv=0
 255     pxkh_tlr=0
 256     pxkv_tlr=0
 257     pxkh_tlx=0
 258     pxkv_tlx=0
 259     call da_radzicevar_cal_tl_fw4wetice(1,qsn,pxabk_all,para_alpha_rxa,para_alpha_rxb,  &
 260                                pxkh,pxkv,ice_abc,ice_bac,fw,pxkh_tlr,          &
 261                                pxkv_tlr,pxkh_tlx,pxkv_tlx,qra,tlopt,           &
 262                                npara_alpharxa)
 263
 264      if(in0s.eq.1) then
 265       ! for two moment microphysics scheme, not yet completed
 266      else
 267        if(tlopt==0) zsh= para1sg*sonv**(-0.75)*(rhoair*wsnow)**(1.75)*pxkh
 268
 269
 270      endif
 271
 272      if(dualpol_opt==1) then
 273        if(in0s.eq.1) then
 274          ! for two moment microphysics scheme, not yet completed
 275        else
 276          if(tlopt==0) zsv= para1sg*sonv**(-0.75)*(rhoair*wsnow)**(1.75)*pxkv
 277
 278        endif
 279      endif
 280
 281 !  -------------------for dry snow----------------------
 282      call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhos,dsnow,cs,ds,alphas,   &
 283                         mm3todBZ,lambda,Kw2,pi,sonv)
 284      ice_abc_d(1)=Asd
 285      ice_abc_d(2)=Bsd
 286      ice_abc_d(3)=Csd
 287      ice_bac_d(1)=Bsd
 288      ice_bac_d(2)=Asd
 289
 290      pxabk_all(1)=alpha_rdsa**2
 291      pxabk_all(2)=alpha_rdsb**2
 292      pxabk_all(3)=alpha_rdsb*alpha_rdsa
 293      call da_radzicevar_pkx(ice_abc_d,pxabk_all,pxkh)
 294      call da_radzicevar_pkx(ice_bac_d,pxabk_all,pxkv)
 295
 296      if(in0s.eq.1) then
 297        ! for two moment microphysics scheme, not yet completed
 298      else
 299        if(tlopt==0) zdsh=para1sg*sonv**(-0.75)*(rhoair*dsnow)**(1.75)*pxkh
 300
 301      endif
 302
 303      if(dualpol_opt==1) then
 304        if(in0s.eq.1) then
 305         ! for two moment microphysics scheme, not yet completed
 306        else
 307          if(tlopt==0) zdsv=para1sg*sonv**(-0.75)*(rhoair*dsnow)**(1.75)*pxkv
 308
 309        endif
 310      endif
 311      zss=zsh+zdsh
 312
 313 ! ==================FOR graupel==========================
 314 ! -------------------for wet graupel---------------------
 315
 316     call da_radzicevar_waterfraction(qra,qgr,fw)
 317     rhowg=rhog*(1.-fw**2)+rhor*fw**2
 318
 319     call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhowg,wgr,cice,dg,alphag,   &
 320                        mm3todBZ,lambda,Kw2,pi,gonv)
 321
 322     call da_radzicevar_sigma_in_abc(qgr,fw,2,sigma)  ! for snow
 323     call da_radzicevar_calc_ice_abc(phimean,sigma,ice_abc)
 324
 325     ice_bac(1)=ice_abc(2)
 326     ice_bac(2)=ice_abc(1)
 327     ice_bac(3)=ice_abc(3)
 328     pxkh=0
 329     pxkv=0
 330     pxkh_tlr=0
 331     pxkv_tlr=0
 332     pxkh_tlx=0
 333     pxkv_tlx=0
 334     call da_radzicevar_cal_tl_fw4wetice(2,qgr,pxabk_all,para_alpha_rxa,para_alpha_rxb,       &
 335                                pxkh,pxkv,ice_abc,ice_bac,fw,pxkh_tlr,          &
 336                                pxkv_tlr,pxkh_tlx,pxkv_tlx,qra,tlopt,           &
 337                                npara_alpharxa)
 338
 339     if(in0g.eq.1) then
 340        ! for two moment microphysics scheme, not yet completed
 341     else
 342       if(tlopt==0) zgh= para1sg*gonv**(-0.75)*(rhoair*wgr)**(1.75)*pxkh
 343
 344
 345     endif
 346
 347     if(dualpol_opt==1) then
 348        if(in0g.eq.1) then
 349        ! for two moment microphysics scheme, not yet completed
 350        else
 351          if(tlopt==0) zgv= para1sg*sonv**(-0.75)*(rhoair*wgr)**(1.75)*pxkv
 352
 353        endif
 354      endif
 355
 356 ! --------------------------for dry graupel----------------------
 357
 358      call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhog,dgr,cice,dg,alphag,   &
 359                         mm3todBZ,lambda,Kw2,pi,gonv)
 360
 361      ice_abc_d(1)=Ahd
 362      ice_abc_d(2)=Bhd
 363      ice_abc_d(3)=Chd
 364      ice_bac_d(1)=Bhd
 365      ice_bac_d(2)=Ahd
 366
 367      pxabk_all(1)=alpha_rdha**2
 368      pxabk_all(2)=alpha_rdhb**2
 369      pxabk_all(3)=alpha_rdhb*alpha_rdha
 370      call da_radzicevar_pkx(ice_abc_d,pxabk_all,pxkh)
 371      call da_radzicevar_pkx(ice_bac_d,pxabk_all,pxkv)
 372
 373      if(in0s.eq.1) then
 374       ! for two moment microphysics scheme, not yet completed
 375      else
 376        if(tlopt==0) zdgh=para1sg*gonv**(-0.75)*(rhoair*dgr)**(1.75)*pxkh
 377
 378      endif
 379
 380      if(dualpol_opt==1) then
 381        if(in0s.eq.1) then
 382         ! for two moment microphysics scheme, not yet completed
 383        else
 384          if(tlopt==0) zdgv=para1sg*gonv**(-0.75)*(rhoair*dgr)**(1.75)*pxkv
 385
 386        endif
 387      endif
 388
 389      zhs=zgh+zdgh
 390 ! ============================done=================================
 391   endif  !(if rf_noice == 0)
 392
 393   z_e =zrh+zsh+zgh+zdsh+zdgh
 394   zv_e=zrv+zsv+zgv+zdsv+zdgv
 395 !-------------------------------------
 396 ! Convert to dBZ
 397 !-------------------------------------
 398   if(tlopt==0) then
 399      dbz = 10. * log10(z_e)
 400   endif
 401 ! save z_e mm^6 mm^-3
 402   zmm=z_e
 403
 404 end subroutine da_radzicevar