var/da/da_radar/da_radzicevar_tl.inc

   1 subroutine da_radzicevar_tl(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   real           :: tc
 137
 138   !---------------------------------------------------------------------------------
 139
 140   data para_alpha_rxa(:,1)/0.194e-4,7.094e-4,2.135e-4,-5.225e-4,0,0,0/
 141   data para_alpha_rxb(:,1)/0.191e-4,6.916e-4,-2.841e-4,-1.160e-4,0,0,0/
 142   !----hail--------------
 143   !data para_alpha_rxa(:,2)/0.191e-3,2.39e-3,-12.57e-3,38.71e-3,-65.53e-3,56.16e-3,  &
 144   !                        -18.98e-3/
 145   !data para_alpha_rxb(:,2)/0.165e-3,1.72e-3,-9.92e-3,32.15e-3,-56.0e-3,48.83e-3,   &
 146   !                        -16.69e-3/
 147   !----graupel-----------
 148   data para_alpha_rxa(:,2)/1.05E-04,1.82E-03,-3.77E-03,-7.97E-04,1.63E-02,-2.20E-02, &
 149                            8.74E-03/
 150   data para_alpha_rxb(:,2)/9.25E-05,1.93E-03,-9.79E-03,2.92E-02,-4.82E-02,3.93E-02, &
 151                           -1.22E-02/
 152
 153
 154   qthres=rf_qthres
 155
 156   qra=qra0
 157   qsn=qsn0
 158   qgr=qgr0
 159   tmk=tmk0
 160   qvp=qvp0
 161 !  if(qra0<qthres) qra=qthres
 162 !  if(qsn0<qthres) qsn=qthres
 163 !  if(qgr0<qthres) qgr=qthres
 164
 165   tc=tmk-273.15
 166   if(tc<=0.0) then
 167      qra=max(qthres*(1+exp(-(tc-5.0)**2/6.25)),qra)
 168      qsn=max(qthres*2,qsn)
 169      qgr=max(qthres*2,qgr)
 170   else
 171      qra=max(qthres*2,qra)
 172      qsn=max(qthres*(1+exp(-(tc-5.0)**2/6.25)),qsn)
 173      qgr=max(qthres*(1+exp(-(tc-5.0)**2/6.25)),qgr)
 174   endif
 175
 176   if(tlopt>=1.and.gropt>=1) then
 177     dqra=0
 178     dqgr=0
 179     dqsn=0
 180     dqvp=0
 181     dtmk=0
 182     dqnr=0
 183     dqns=0
 184     dqng=0
 185   endif
 186
 187   zrh=0
 188   zsh=0
 189   zgh=0
 190   zdsh=0
 191   zdgh=0
 192   zrv=0
 193   zsv=0
 194   zgv=0
 195   zdsv=0
 196   zdgv=0
 197
 198 !  if(qra<qthres.and.qsn<qthres.and.qgr<qthres) then
 199 !     zrs=0
 200 !     zss=0
 201 !     zhs=0
 202 !     ref=0
 203 !     zmm=0
 204 !     dbz=0
 205 !     return
 206 !  endif
 207
 208   call da_radzicevar_rhoair_tl(tlopt,rhoair,prs,rgas,tmk,qvp,pdfrhot,pdfrhoq)
 209
 210 !-------------------------------------------------------
 211 ! Calculate variable intercept parameters if wanted
 212 !-------------------------------------------------------
 213   call da_radzicevar_prepare_interceptpara(in0s,in0g,in0r,rn0_s,rn0_g,rn0_r,sonv,gonv,ronv)
 214 !-------------------------------------------------------
 215 ! Calculate mixing ratios (pure rain, dry snow/graupel, wet snow/graupel)
 216 !-------------------------------------------------------
 217   call da_radzicevar_prepare_mixingratios(tlopt,prain_coef,dsnow_coef,dgr_coef, &
 218                             prain,dsnow,wsnow,dgr,wgr,            &
 219                             qra,qgr,qsn,qthres,                   &
 220                             pdfrrs,pdfrrg,pdfsrs,pdfgrg           &
 221                             )
 222
 223   call da_radzicevar_prepare_zmm_adj(tlopt,gropt,zmm,zmm_ref)
 224 ! ==================FOR RAIN=============================
 225   call da_radzicevar_parameter_zrx(para10r,para14r,para_pr,rhoair,rhor,prain,   &
 226                      beta_ra,alpha_ra,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
 233     if(tlopt>=1.and.gropt==0) call da_radzicevar_rain_tl(zrh,para_pr,dqra,dqsn,dqgr,dtmk,dqvp,rhoair,  &
 234                                                 prain_coef,zmm_ref,pdfrrs,pdfsrs,pdfgrg,     &
 235                                                 pdfrhot,pdfrhoq,pdfrrg,qra,prain)
 236
 237   endif
 238
 239   if(dualpol_opt==1) then
 240     call da_radzicevar_parameter_zrx(para10r,para14r,para_pr,rhoair,rhor,prain,   &
 241                        beta_rb,alpha_rb,mm3todBZ,lambda,Kw2,pi,     &
 242                        ronv)
 243
 244     if(in0r.eq.1) then
 245      ! for two moment microphysics scheme, not yet completed
 246     else
 247
 248       if(tlopt>=1.and.gropt==0) call da_radzicevar_rain_tl(zrv,para_pr,dqra,dqsn,dqgr,dtmk,dqvp,rhoair,  &
 249                                                   prain_coef,zmm_ref,pdfrrs,pdfsrs,pdfgrg,     &
 250                                                   pdfrhot,pdfrhoq,pdfrrg,qra,prain)
 251
 252     endif
 253   endif
 254   zrs=zrh
 255 ! ==================FOR snow=============================
 256   if(rf_noice==0) then
 257 ! -------------------for wet snow----------------------
 258
 259     call da_radzicevar_waterfraction(qra,qsn,fw)
 260     rhows=rhos*(1.-fw**2)+rhor*fw**2
 261
 262     call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhows,wsnow,cs,ds,alphas,   &
 263                        mm3todBZ,lambda,Kw2,pi,sonv)
 264
 265     call da_radzicevar_sigma_in_abc(qsn,fw,1,sigma)  ! for snow
 266     call da_radzicevar_calc_ice_abc(phimean,sigma,ice_abc)
 267
 268     ice_bac(1)=ice_abc(2)
 269     ice_bac(2)=ice_abc(1)
 270     ice_bac(3)=ice_abc(3)
 271     pxkh=0
 272     pxkv=0
 273     pxkh_tlr=0
 274     pxkv_tlr=0
 275     pxkh_tlx=0
 276     pxkv_tlx=0
 277     call da_radzicevar_cal_tl_fw4wetice(1,qsn,pxabk_all,para_alpha_rxa,para_alpha_rxb,  &
 278                                pxkh,pxkv,ice_abc,ice_bac,fw,pxkh_tlr,          &
 279                                pxkv_tlr,pxkh_tlx,pxkv_tlx,qra,tlopt,           &
 280                                npara_alpharxa)
 281
 282      if(in0s.eq.1) then
 283       ! for two moment microphysics scheme, not yet completed
 284      else
 285
 286        if(tlopt>=1.and.gropt==0) call da_radzicevar_wetice_tl(zsh,para_pdx_dq,para_pdx_df,pxkh,     &
 287                                                 pxkh_tlr,pxkh_tlx, dqsn,dqra,dtmk,    &
 288                                                 dqvp,qsn,wsnow,rhoair,dsnow_coef,     &
 289                                                 zmm_ref,pdfsrs,pdfrrs,pdfrhot,pdfrhoq)
 290
 291
 292      endif
 293
 294      if(dualpol_opt==1) then
 295        if(in0s.eq.1) then
 296          ! for two moment microphysics scheme, not yet completed
 297        else
 298
 299          if(tlopt>=1.and.gropt==0) call da_radzicevar_wetice_tl(zsv,para_pdx_dq,para_pdx_df,pxkv,  &
 300                                                   pxkv_tlr,pxkv_tlx,dqsn,dqra,dtmk,  &
 301                                                   dqvp,qsn,wsnow,rhoair,dsnow_coef,  &
 302                                                   zmm_ref,pdfsrs,pdfrrs,pdfrhot,pdfrhoq)
 303
 304        endif
 305      endif
 306
 307 !  -------------------for dry snow----------------------
 308      call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhos,dsnow,cs,ds,alphas,   &
 309                         mm3todBZ,lambda,Kw2,pi,sonv)
 310      ice_abc_d(1)=Asd
 311      ice_abc_d(2)=Bsd
 312      ice_abc_d(3)=Csd
 313      ice_bac_d(1)=Bsd
 314      ice_bac_d(2)=Asd
 315
 316      pxabk_all(1)=alpha_rdsa**2
 317      pxabk_all(2)=alpha_rdsb**2
 318      pxabk_all(3)=alpha_rdsb*alpha_rdsa
 319      call da_radzicevar_pkx(ice_abc_d,pxabk_all,pxkh)
 320      call da_radzicevar_pkx(ice_bac_d,pxabk_all,pxkv)
 321
 322      if(in0s.eq.1) then
 323        ! for two moment microphysics scheme, not yet completed
 324      else
 325
 326        if(tlopt>=1.and.gropt==0) call da_radzicevar_dryice_tl(zdsh,para_pdx_dq,pxkh,dqsn,dqra,dtmk,dqvp,      &
 327                                                 rhoair,dsnow_coef,zmm_ref,dsnow,qsn,pdfsrs,     &
 328                                                 pdfrrs,pdfrhot,pdfrhoq)
 329
 330      endif
 331
 332      if(dualpol_opt==1) then
 333        if(in0s.eq.1) then
 334         ! for two moment microphysics scheme, not yet completed
 335        else
 336
 337          if(tlopt>=1.and.gropt==0) call da_radzicevar_dryice_tl(zdsv,para_pdx_dq,pxkv,dqsn,dqra,dtmk,dqvp,      &
 338                                                   rhoair,dsnow_coef,zmm_ref,dsnow,qsn,pdfsrs,     &
 339                                                   pdfrrs,pdfrhot,pdfrhoq)
 340
 341        endif
 342      endif
 343      zss=zsh+zdsh
 344
 345 ! ==================FOR graupel==========================
 346 ! -------------------for wet graupel---------------------
 347
 348     call da_radzicevar_waterfraction(qra,qgr,fw)
 349     rhowg=rhog*(1.-fw**2)+rhor*fw**2
 350
 351     call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhowg,wgr,cice,dg,alphag,   &
 352                        mm3todBZ,lambda,Kw2,pi,gonv)
 353
 354     call da_radzicevar_sigma_in_abc(qgr,fw,2,sigma)  ! for snow
 355     call da_radzicevar_calc_ice_abc(phimean,sigma,ice_abc)
 356
 357     ice_bac(1)=ice_abc(2)
 358     ice_bac(2)=ice_abc(1)
 359     ice_bac(3)=ice_abc(3)
 360     pxkh=0
 361     pxkv=0
 362     pxkh_tlr=0
 363     pxkv_tlr=0
 364     pxkh_tlx=0
 365     pxkv_tlx=0
 366     call da_radzicevar_cal_tl_fw4wetice(2,qgr,pxabk_all,para_alpha_rxa,para_alpha_rxb,       &
 367                                pxkh,pxkv,ice_abc,ice_bac,fw,pxkh_tlr,          &
 368                                pxkv_tlr,pxkh_tlx,pxkv_tlx,qra,tlopt,           &
 369                                npara_alpharxa)
 370
 371     if(in0g.eq.1) then
 372        ! for two moment microphysics scheme, not yet completed
 373     else
 374
 375       if(tlopt>=1.and.gropt==0) call da_radzicevar_wetice_tl(zgh,para_pdx_dq,para_pdx_df,pxkh,pxkh_tlr, &
 376                                                pxkh_tlx,dqgr,dqra,dtmk,dqvp,qgr,wgr,      &
 377                                                rhoair,dgr_coef,zmm_ref,pdfgrg,            &
 378                                                pdfrrg,pdfrhot,pdfrhoq)
 379
 380
 381     endif
 382
 383     if(dualpol_opt==1) then
 384        if(in0g.eq.1) then
 385        ! for two moment microphysics scheme, not yet completed
 386        else
 387
 388          if(tlopt>=1.and.gropt==0) call da_radzicevar_wetice_tl(zgv,para_pdx_dq,para_pdx_df,pxkv,           &
 389                                                   pxkv_tlr,pxkv_tlx,dqgr,dqra,dtmk,dqvp,qgr,  &
 390                                                   wgr,rhoair,dgr_coef,zmm_ref,pdfgrg,         &
 391                                                   pdfrrg,pdfrhot,pdfrhoq)
 392
 393        endif
 394      endif
 395
 396 ! --------------------------for dry graupel----------------------
 397
 398      call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhog,dgr,cice,dg,alphag,   &
 399                         mm3todBZ,lambda,Kw2,pi,gonv)
 400
 401      ice_abc_d(1)=Ahd
 402      ice_abc_d(2)=Bhd
 403      ice_abc_d(3)=Chd
 404      ice_bac_d(1)=Bhd
 405      ice_bac_d(2)=Ahd
 406
 407      pxabk_all(1)=alpha_rdha**2
 408      pxabk_all(2)=alpha_rdhb**2
 409      pxabk_all(3)=alpha_rdhb*alpha_rdha
 410      call da_radzicevar_pkx(ice_abc_d,pxabk_all,pxkh)
 411      call da_radzicevar_pkx(ice_bac_d,pxabk_all,pxkv)
 412
 413      if(in0s.eq.1) then
 414       ! for two moment microphysics scheme, not yet completed
 415      else
 416
 417        if(tlopt>=1.and.gropt==0) call da_radzicevar_dryice_tl(zdgh,para_pdx_dq,pxkh,dqgr,dqra,dtmk,dqvp,     &
 418                                                 rhoair,dgr_coef,zmm_ref,dgr,qgr,pdfgrg,        &
 419                                                 pdfrrg,pdfrhot,pdfrhoq)
 420
 421      endif
 422
 423      if(dualpol_opt==1) then
 424        if(in0s.eq.1) then
 425         ! for two moment microphysics scheme, not yet completed
 426        else
 427
 428          if(tlopt>=1.and.gropt==0) call da_radzicevar_dryice_tl(zdgv,para_pdx_dq,pxkv,dqgr,dqra,dtmk,dqvp,   &
 429                                                   rhoair,dgr_coef,zmm_ref,dgr,qgr,pdfgrg,      &
 430                                                   pdfrrg,pdfrhot,pdfrhoq)
 431
 432        endif
 433      endif
 434
 435      zhs=zgh+zdgh
 436 ! ============================done=================================
 437   endif  !(if rf_noice == 0)
 438
 439   z_e =zrh+zsh+zgh+zdsh+zdgh
 440   zv_e=zrv+zsv+zgv+zdsv+zdgv
 441 !-------------------------------------
 442 ! Convert to dBZ
 443 !-------------------------------------
 444
 445   if(tlopt>=1) then
 446      if(zmm>=3.16e-2) then
 447         dbz = 10./(zmm*log(10.0))*z_e !ref
 448      else
 449         dbz =0.0
 450      endif
 451   endif
 452
 453
 454 end subroutine da_radzicevar_tl