var/da/da_radar/da_radzicevar_adj.inc

   1 subroutine da_radzicevar_adj(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
 108   real           :: zh_e,zv_e
 109
 110   ! for dry graupel
 111   REAL,PARAMETER :: sigmahd = 1.0472
 112   REAL,PARAMETER :: Ahd = 0.4308
 113   REAL,PARAMETER :: Bhd = 0.3192
 114   REAL,PARAMETER :: Chd = 0.1250
 115   REAL,PARAMETER :: Dhd = 0.3750
 116   REAL,PARAMETER :: Ckhd = 0.1116
 117   ! for dry snow
 118   REAL,PARAMETER :: sigmas = 0.3491
 119   REAL,PARAMETER :: Asd = 0.8140
 120   REAL,PARAMETER :: Bsd = 0.0303
 121   REAL,PARAMETER :: Csd = 0.0778
 122   REAL,PARAMETER :: Dsd = 0.4221
 123   REAL,PARAMETER :: Cksd = 0.7837
 124
 125   real           :: ice_abc_d(3),ice_bac_d(3)
 126   real           :: pdfrrs,pdfrrg,pdfsrs,pdfgrg    ! partial derivative of F (r&s,r&g) with respect to (r,s,g)
 127   real           :: pdfrhot,pdfrhoq                ! partial derivative of rho for t and qv
 128
 129   real           :: prain_coef,dsnow_coef,dgr_coef ! coefficient to determine the ratio of
 130                                                    !pure rain, dry swno/graupel in total qr and qs/qg, respectively.
 131   real           :: qthres
 132
 133   real           :: temr01,temr02,temr03    ! temporary variable for real
 134   integer        :: temi01,temi02,temi03    ! temporary variable for integer
 135
 136   real           :: para_pr,para_pdx_dq,para_pdx_df        ! for
 137   real           :: tc
 138
 139   !---------------------------------------------------------------------------------
 140
 141   data para_alpha_rxa(:,1)/0.194e-4,7.094e-4,2.135e-4,-5.225e-4,0,0,0/
 142   data para_alpha_rxb(:,1)/0.191e-4,6.916e-4,-2.841e-4,-1.160e-4,0,0,0/
 143   !----hail--------------
 144   !data para_alpha_rxa(:,2)/0.191e-3,2.39e-3,-12.57e-3,38.71e-3,-65.53e-3,56.16e-3,  &
 145   !                        -18.98e-3/
 146   !data para_alpha_rxb(:,2)/0.165e-3,1.72e-3,-9.92e-3,32.15e-3,-56.0e-3,48.83e-3,   &
 147   !                        -16.69e-3/
 148   !----graupel-----------
 149   data para_alpha_rxa(:,2)/1.05E-04,1.82E-03,-3.77E-03,-7.97E-04,1.63E-02,-2.20E-02, &
 150                            8.74E-03/
 151   data para_alpha_rxb(:,2)/9.25E-05,1.93E-03,-9.79E-03,2.92E-02,-4.82E-02,3.93E-02, &
 152                           -1.22E-02/
 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   tc=tmk-273.15
 165   if(tc<=0.0) then
 166     qra=max(qthres*(1+exp(-(tc-5.0)**2/6.25)),qra)
 167     qsn=max(qthres*2,qsn)
 168     qgr=max(qthres*2,qgr)
 169   else
 170     qra=max(qthres*2,qra)
 171     qsn=max(qthres*(1+exp(-(tc-5.0)**2/6.25)),qsn)
 172     qgr=max(qthres*(1+exp(-(tc-5.0)**2/6.25)),qgr)
 173   endif
 174
 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>=1) call da_radzicevar_rain_adj(para_pr,dqra,dqsn,dqgr,dtmk,dqvp,rhoair,  &
 234                                                 prain_coef,zmm_ref,pdfrrs,pdfsrs,pdfgrg, &
 235                                                 pdfrhot,pdfrhoq,pdfrrg,qra,prain)
 236   endif
 237
 238   if(dualpol_opt==1) then
 239     call da_radzicevar_parameter_zrx(para10r,para14r,para_pr,rhoair,rhor,prain,   &
 240                        beta_rb,alpha_rb,mm3todBZ,lambda,Kw2,pi,     &
 241                        ronv)
 242
 243     if(in0r.eq.1) then
 244      ! for two moment microphysics scheme, not yet completed
 245     else
 246
 247       if(tlopt>=1.and.gropt>=1) call da_radzicevar_rain_adj(para_pr,dqra,dqsn,dqgr,dtmk,dqvp,rhoair,  &
 248                                                   prain_coef,zmm_ref,pdfrrs,pdfsrs,pdfgrg, &
 249                                                   pdfrhot,pdfrhoq,pdfrrg,qra,prain)
 250     endif
 251   endif
 252
 253 ! ==================FOR snow=============================
 254   if(rf_noice==0) then
 255 ! -------------------for wet snow----------------------
 256
 257     call da_radzicevar_waterfraction(qra,qsn,fw)
 258     rhows=rhos*(1.-fw**2)+rhor*fw**2
 259
 260     call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhows,wsnow,cs,ds,alphas,   &
 261                        mm3todBZ,lambda,Kw2,pi,sonv)
 262
 263     call da_radzicevar_sigma_in_abc(qsn,fw,1,sigma)  ! for snow
 264     call da_radzicevar_calc_ice_abc(phimean,sigma,ice_abc)
 265
 266     ice_bac(1)=ice_abc(2)
 267     ice_bac(2)=ice_abc(1)
 268     ice_bac(3)=ice_abc(3)
 269     pxkh=0
 270     pxkv=0
 271     pxkh_tlr=0
 272     pxkv_tlr=0
 273     pxkh_tlx=0
 274     pxkv_tlx=0
 275     call da_radzicevar_cal_tl_fw4wetice(1,qsn,pxabk_all,para_alpha_rxa,para_alpha_rxb,  &
 276                                pxkh,pxkv,ice_abc,ice_bac,fw,pxkh_tlr,          &
 277                                pxkv_tlr,pxkh_tlx,pxkv_tlx,qra,tlopt,           &
 278                                npara_alpharxa)
 279
 280      if(in0s.eq.1) then
 281       ! for two moment microphysics scheme, not yet completed
 282      else
 283
 284        if(tlopt>=1.and.gropt>=1) call da_radzicevar_wetice_adj(para_pdx_dq,para_pdx_df,pxkh,pxkh_tlr,  &
 285                                                  pxkh_tlx,dqsn,dqra,dtmk,dqvp,qsn,wsnow, &
 286                                                  rhoair,dsnow_coef,zmm_ref,pdfsrs,       &
 287                                                  pdfrrs,pdfrhot,pdfrhoq)
 288
 289      endif
 290
 291      if(dualpol_opt==1) then
 292        if(in0s.eq.1) then
 293          ! for two moment microphysics scheme, not yet completed
 294        else
 295
 296          if(tlopt>=1.and.gropt>=1) call da_radzicevar_wetice_adj(para_pdx_dq,para_pdx_df,pxkv,pxkv_tlr,   &
 297                                                    pxkv_tlx, dqsn,dqra,dtmk,dqvp,qsn,wsnow, &
 298                                                    rhoair,dsnow_coef,zmm_ref,pdfsrs,        &
 299                                                    pdfrrs,pdfrhot,pdfrhoq)
 300        endif
 301      endif
 302
 303 !  -------------------for dry snow----------------------
 304      call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhos,dsnow,cs,ds,alphas,   &
 305                         mm3todBZ,lambda,Kw2,pi,sonv)
 306      ice_abc_d(1)=Asd
 307      ice_abc_d(2)=Bsd
 308      ice_abc_d(3)=Csd
 309      ice_bac_d(1)=Bsd
 310      ice_bac_d(2)=Asd
 311
 312      pxabk_all(1)=alpha_rdsa**2
 313      pxabk_all(2)=alpha_rdsb**2
 314      pxabk_all(3)=alpha_rdsb*alpha_rdsa
 315      call da_radzicevar_pkx(ice_abc_d,pxabk_all,pxkh)
 316      call da_radzicevar_pkx(ice_bac_d,pxabk_all,pxkv)
 317
 318      if(in0s.eq.1) then
 319        ! for two moment microphysics scheme, not yet completed
 320      else
 321
 322        if(tlopt>=1.and.gropt>=1) call da_radzicevar_dryice_adj(para_pdx_dq,pxkh,dqsn,dqra,dtmk,dqvp,         &
 323                                                  rhoair,dsnow_coef,zmm_ref,dsnow,qsn,pdfsrs,   &
 324                                                  pdfrrs,pdfrhot,pdfrhoq)
 325      endif
 326
 327      if(dualpol_opt==1) then
 328        if(in0s.eq.1) then
 329         ! for two moment microphysics scheme, not yet completed
 330        else
 331
 332          if(tlopt>=1.and.gropt>=1) call da_radzicevar_dryice_adj(para_pdx_dq,pxkv,dqsn,dqra,dtmk,dqvp,         &
 333                                                    rhoair,dsnow_coef,zmm_ref,dsnow,qsn,pdfsrs,   &
 334                                                    pdfrrs,pdfrhot,pdfrhoq)
 335        endif
 336      endif
 337
 338 ! ==================FOR graupel==========================
 339 ! -------------------for wet graupel---------------------
 340
 341     call da_radzicevar_waterfraction(qra,qgr,fw)
 342     rhowg=rhog*(1.-fw**2)+rhor*fw**2
 343
 344     call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhowg,wgr,cice,dg,alphag,   &
 345                        mm3todBZ,lambda,Kw2,pi,gonv)
 346
 347     call da_radzicevar_sigma_in_abc(qgr,fw,2,sigma)  ! for snow
 348     call da_radzicevar_calc_ice_abc(phimean,sigma,ice_abc)
 349
 350     ice_bac(1)=ice_abc(2)
 351     ice_bac(2)=ice_abc(1)
 352     ice_bac(3)=ice_abc(3)
 353     pxkh=0
 354     pxkv=0
 355     pxkh_tlr=0
 356     pxkv_tlr=0
 357     pxkh_tlx=0
 358     pxkv_tlx=0
 359     call da_radzicevar_cal_tl_fw4wetice(2,qgr,pxabk_all,para_alpha_rxa,para_alpha_rxb,       &
 360                                pxkh,pxkv,ice_abc,ice_bac,fw,pxkh_tlr,          &
 361                                pxkv_tlr,pxkh_tlx,pxkv_tlx,qra,tlopt,           &
 362                                npara_alpharxa)
 363
 364     if(in0g.eq.1) then
 365        ! for two moment microphysics scheme, not yet completed
 366     else
 367
 368       if(tlopt>=1.and.gropt>=1) call da_radzicevar_wetice_adj(para_pdx_dq,para_pdx_df,pxkh,pxkh_tlr,      &
 369                                                pxkh_tlx,dqgr,dqra,dtmk,dqvp,qgr,wgr,rhoair, &
 370                                                dgr_coef,zmm_ref,pdfgrg,                     &
 371                                                pdfrrg,pdfrhot,pdfrhoq)
 372
 373     endif
 374
 375     if(dualpol_opt==1) then
 376        if(in0g.eq.1) then
 377        ! for two moment microphysics scheme, not yet completed
 378        else
 379
 380          if(tlopt>=1.and.gropt>=1) call da_radzicevar_wetice_adj(para_pdx_dq,para_pdx_df,pxkv,pxkv_tlr,     &
 381                                                   pxkv_tlx,dqgr,dqra,dtmk,dqvp,qgr,wgr,       &
 382                                                   rhoair,dgr_coef,zmm_ref,pdfgrg,             &
 383                                                   pdfrrg,pdfrhot,pdfrhoq)
 384        endif
 385      endif
 386
 387 ! --------------------------for dry graupel----------------------
 388
 389      call da_radzicevar_parameter_zxx(para1sg,para_pdx_dq,para_pdx_df,rhoair,rhog,dgr,cice,dg,alphag,   &
 390                         mm3todBZ,lambda,Kw2,pi,gonv)
 391
 392      ice_abc_d(1)=Ahd
 393      ice_abc_d(2)=Bhd
 394      ice_abc_d(3)=Chd
 395      ice_bac_d(1)=Bhd
 396      ice_bac_d(2)=Ahd
 397
 398      pxabk_all(1)=alpha_rdha**2
 399      pxabk_all(2)=alpha_rdhb**2
 400      pxabk_all(3)=alpha_rdhb*alpha_rdha
 401      call da_radzicevar_pkx(ice_abc_d,pxabk_all,pxkh)
 402      call da_radzicevar_pkx(ice_bac_d,pxabk_all,pxkv)
 403
 404      if(in0s.eq.1) then
 405       ! for two moment microphysics scheme, not yet completed
 406      else
 407
 408        if(tlopt>=1.and.gropt>=1) call da_radzicevar_dryice_adj(para_pdx_dq,pxkh,dqgr,dqra,dtmk,dqvp,         &
 409                                                 rhoair,dgr_coef,zmm_ref,dgr,qgr,pdfgrg,        &
 410                                                 pdfrrg,pdfrhot,pdfrhoq)
 411      endif
 412
 413      if(dualpol_opt==1) then
 414        if(in0s.eq.1) then
 415         ! for two moment microphysics scheme, not yet completed
 416        else
 417
 418          if(tlopt>=1.and.gropt>=1) call da_radzicevar_dryice_adj(para_pdx_dq,pxkv,dqgr,dqra,dtmk,dqvp,    &
 419                                                   rhoair,dgr_coef,zmm_ref,dgr,qgr,pdfgrg,   &
 420                                                   pdfrrg,pdfrhot,pdfrhoq)
 421        endif
 422      endif
 423
 424 ! ============================done=================================
 425   endif  !(if rf_noice == 0)
 426
 427
 428 end subroutine da_radzicevar_adj