var/da/da_physics/da_moist_phys_lin.inc

   1 subroutine da_moist_phys_lin(grid)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Partition of the hydrometeors via the moist explicit scheme.
   5    !           A warm rain process is used in this subroutine.
   6    !           This is the tangent linear code of the scheme.
   7    !
   8    !  Method: The warm rain process is according to Hsie and Anthes (1984)
   9    !          and Dudhia (1989)
  10    !
  11    !  Assumptions: 1) Model level stored top down.
  12    !---------------------------------------------------------------------------
  13
  14    implicit none
  15
  16    type(domain), intent(inout)               :: grid
  17
  18    real :: T_OLD(ims:ime,jms:jme,kms:kme),T_NEW(ims:ime,jms:jme,kms:kme)
  19    real :: Q_OLD(ims:ime,jms:jme,kms:kme),Q_NEW(ims:ime,jms:jme,kms:kme)
  20    real :: QCW_OLD(ims:ime,jms:jme,kms:kme),QCW_NEW(ims:ime,jms:jme,kms:kme)
  21    real :: QRN_OLD(ims:ime,jms:jme,kms:kme),QRN_NEW(ims:ime,jms:jme,kms:kme)
  22
  23    real    :: EES(kms:kme)
  24    real    :: QVSS(kms:kme)
  25    real    :: EES9(kms:kme)
  26    real    :: QVSS9(kms:kme)
  27    real    :: DT(its:ite,jts:jte,kms:kme)
  28    real    :: QVT(kms:kme)
  29    real    :: QCT(kms:kme)
  30    real    :: QRT(kms:kme)
  31    real    :: TTT(kms:kme)
  32    real    :: QVT9(kms:kme)
  33    real    :: QCT9(kms:kme)
  34    real    :: QRT9(kms:kme)
  35    real    :: TTT9(kms:kme)
  36    real    :: SCR2(kms:kme)
  37    real    :: SCR3(kms:kme)
  38    real    :: SCR4(kms:kme)
  39    real    :: SCR5(kms:kme)
  40    real    :: SCR6(kms:kme)
  41    real    :: DUM31(kms:kme)
  42    real    :: PRA(kms:kme)
  43    real    :: PRC(kms:kme)
  44    real    :: PRD(kms:kme)
  45    real    :: PRE(kms:kme)
  46    real    :: SCR31(kms:kme)
  47    real    :: SCR42(kms:kme)
  48    real    :: SCR71(kms:kme)
  49    real    :: DUM112(kms:kme)
  50    real    :: DUM113(kms:kme)
  51    real    :: DUM211(kms:kme)
  52    real    :: DUM411(kms:kme)
  53    real    :: SCR29(kms:kme)
  54    real    :: SCR39(kms:kme)
  55    real    :: SCR49(kms:kme)
  56    real    :: SCR59(kms:kme)
  57    real    :: SCR69(kms:kme)
  58    real    :: DUM319(kms:kme)
  59    real    :: PRA9(kms:kme)
  60    real    :: PRC9(kms:kme)
  61    real    :: PRD9(kms:kme)
  62    real    :: PRE9(kms:kme)
  63    real    :: SCR319(kms:kme)
  64    real    :: SCR429(kms:kme)
  65    real    :: SCR719(kms:kme)
  66    real    :: DUM1129(kms:kme)
  67    real    :: DUM4119(kms:kme)
  68    real    :: TMP(kms:kme)
  69
  70    integer, parameter :: qcth = 0.5e-3
  71    integer, parameter :: qrth = 1.0e-6
  72    integer, parameter :: alpha = 0.001
  73    integer, parameter :: beta  = 0.0486
  74    integer, parameter :: gamma = 0.002
  75
  76    integer :: i, j, k
  77
  78    real    :: qrth1
  79
  80    if (trace_use) call da_trace_entry("da_moist_phys_lin")
  81
  82    qrth1 = (QRTH*1.0e3)**0.875
  83
  84    T_OLD  (its:ite,jts:jte,kts:kte) = grid%xa%t  (its:ite,jts:jte,kts:kte)
  85    Q_OLD  (its:ite,jts:jte,kts:kte) = grid%xa%q  (its:ite,jts:jte,kts:kte)
  86    QCW_OLD(its:ite,jts:jte,kts:kte) = grid%xa%qcw(its:ite,jts:jte,kts:kte)
  87    QRN_OLD(its:ite,jts:jte,kts:kte) = grid%xa%qrn(its:ite,jts:jte,kts:kte)
  88
  89    !  Preparation
  90
  91    grid%xa%q(its:ite,jts:jte,kts:kte) =grid%xa%qt(its:ite,jts:jte,kts:kte) - grid%xa%qcw(its:ite,jts:jte,kts:kte) - grid%xa %qrn(its:ite,jts:jte,kts:kte)
  92    DT(its:ite,jts:jte,kts:kte) = grid%xb%delt(its:ite,jts:jte,kts:kte)
  93
  94    do j=jts,jte
  95       do i=its,ite
  96          do K=kts,kte
  97
  98             if (dt(i,j,k) <= 0.0) cycle
  99
 100             if ( grid%xb%t(I,J,K) > TO )then
 101                EES(K)=SVP1*EXP(SVP2*(grid%xb%t(I,J,K)-SVPT0)/(grid%xb%t(I,J,K)-SVP3))
 102                EES9(K)=EES(K)*SVP2*(SVPT0-SVP3)/((grid%xb%t(I,J,K)-SVP3) * (grid%xb%t(I,J,K)-SVP3))*grid%xa%t(I,J,K)
 103             else
 104                EES(K)=.611*EXP(22.514-6.15E3/grid%xb%t(I,J,K))
 105                EES9(K)=EES(K)*6.15E3/(grid%xb%t(I,J,K)*grid%xb%t(I,J,K))*grid%xa%t(I,J,K)
 106             end if
 107
 108             TMP(K)=622.0/((grid%xb%p(I,J,K)-EES(K))**2)
 109             QVSS9(K)=TMP(K)*grid%xb%p(I,J,K)*EES9(K) - TMP(K)*EES(K)*grid%xa%p(I,J,K)
 110             QVSS(K)=622.0*EES(K)/(grid%xb%p(I,J,K)-EES(K))
 111
 112             SCR49(K)=grid%xa%q(I,J,K)/QVSS(K)-grid%xb%q(I,J,K)/QVSS(K)**2*QVSS9(K)
 113             SCR4(K)=grid%xb%q(I,J,K)/QVSS(K)
 114
 115             if (grid%xb%qcw(I,J,K) > 0.0) then
 116                SCR29(K)=grid%xa%qcw(I,J,K)
 117                SCR2(K)=grid%xb%qcw(I,J,K)
 118             else
 119                SCR29(K)=0.0
 120                SCR2(K)=0.0
 121             end if
 122             if (grid%xb%qrn(I,J,K) > 1.0e-25) then
 123                SCR39(K)=grid%xa%qrn(I,J,K)
 124                SCR3(K)=grid%xb%qrn(I,J,K)
 125             else
 126                SCR39(K)=0.0
 127                SCR3(K)=1.0E-25
 128             end if
 129             SCR59(K)=grid%xa%q(I,J,K)/SCR4(K)-grid%xb%q(I,J,K)/SCR4(K)**2*SCR49(K)
 130             SCR5(K)=grid%xb%q(I,J,K)/SCR4(K)
 131
 132             SCR69(K)=grid%xa%p(I,J,K)/(gas_constant*grid%xb%t(I,J,K))-grid%xb%p(I,J,K)/  &
 133                      (gas_constant*grid%xb%t(I,J,K)**2)*grid%xa%t(I,J,K)
 134             SCR6(K)=grid%xb%p(I,J,K)/(gas_constant*grid%xb%t(I,J,K))
 135
 136             DUM319(K)=-XLV1*grid%xa%t(I,J,K)
 137             DUM31(K)=3.1484E6-XLV1*grid%xb%t(I,J,K)
 138
 139             ! Auto conversion
 140
 141             if (scr2(k) >= qcth) then
 142                prc9(k) = alpha * scr29(k)
 143                prc(k) = alpha * (scr2(k) - qcth)
 144             else
 145                prc9(k) = 0.0
 146                prc(k) = 0.0
 147             end if
 148
 149             ! Accretion
 150
 151             if (SCR2(k) > 0.0 .and. SCR3(k) > QRTH ) then
 152                PRA9(K) = gamma * 0.875 * SCR2(k) * (SCR3(K)*1.0e3)**(-0.125) * 1.0e3 * SCR39(K)  &
 153                             + gamma * SCR29(k) * (SCR3(K)*1.0e3)**0.875
 154                PRA(k) = gamma * SCR2(k) * (SCR3(k)*1.0e3)**0.875
 155             else if (SCR2(k) > 0.0 .and. SCR3(k) <= QRTH ) then
 156                PRA9(K) = gamma * SCR29(k) * qrth1
 157                PRA(k) = gamma * SCR2(k) * qrth1
 158             else
 159                PRA9(K) = 0.0
 160                PRA(k) = 0.0
 161             end if
 162
 163          end do
 164
 165          call da_evapo_lin(DT(i,j,:),SCR3,SCR5,grid%xb%q(I,J,:),PRE,SCR6,  &
 166                            SCR39,SCR59,grid%xa%q(I,J,:),PRE9,SCR69, &
 167                            kts,kte,kms,kme)
 168
 169          do K=kts, kte
 170
 171             if (dt(i,j,k) <= 0.0) cycle
 172
 173             !  Readjust
 174
 175             DUM112(K)=(PRC(k)+PRA(k))*dt(i,j,k)
 176             if (DUM112(K) > SCR2(k)) then
 177                DUM1129(K)=(PRC9(k)+PRA9(k))*dt(i,j,k)
 178                PRC9(K)=SCR29(K)*PRC(K)/DUM112(K)  &
 179                       +PRC9(K)*SCR2(K)/DUM112(K)  &
 180                       -SCR2(K)*PRC(K)/(DUM112(K)*DUM112(K))*DUM1129(K)
 181                PRC(k)=SCR2(K)*PRC(K)/DUM112(K)
 182                PRA9(K)=SCR29(K)*PRA(K)/DUM112(K)  &
 183                       +PRA9(K)*SCR2(K)/DUM112(K)  &
 184                       -SCR2(K)*PRA(K)/(DUM112(K)*DUM112(K))*DUM1129(K)
 185                PRA(k)=SCR2(K)*PRA(K)/DUM112(K)
 186             end if
 187             QVT9(K)=-PRE9(K)
 188             QVT(K)=-PRE(K)
 189             QCT9(K)=-PRC9(K)-PRA9(K)
 190             QCT(K)=-PRC(K)-PRA(K)
 191             QRT9(K)=PRC9(K)+PRA9(K)+PRE9(K)
 192             QRT(K)=PRC(K)+PRA(K)+PRE(K)
 193             if (grid%xb%t(I,J,K).GT.TO)then
 194                DUM4119(K)=DUM319(K)
 195                DUM411(K)=DUM31(K)
 196             else
 197                DUM4119(K)=0.0
 198                DUM411(K)=XLS
 199             end if
 200             PRD9(K)=cp*0.887*grid%xa%q(I,J,K)
 201             PRD(K)=cp*(1.0+0.887*grid%xb%q(I,J,K))
 202             TTT9(K)=-DUM4119(K)*QVT(K)/PRD(K)  &
 203                    -QVT9(K)*DUM411(K)/PRD(K)  &
 204                    +DUM411(K)*QVT(K)/(PRD(K)*PRD(K))*PRD9(K)
 205             TTT(K)=-DUM411(K)*QVT(K)/PRD(K)
 206
 207             DUM113(K)=grid%xb%q(I,J,K)+dt(i,j,k)*QVT(K)
 208             if (DUM113(K) > 1.0e-12 ) then
 209                SCR429(K)=grid%xa%q(I,J,K)+dt(i,j,k)*QVT9(K)
 210                SCR42(K)=DUM113(K)
 211             else
 212                SCR429(K)=0.0
 213                SCR42(K)=1.0e-12
 214             end if
 215             DUM211(K)=grid%xb%qcw(I,J,K)+QCT(K)*dt(i,j,k)
 216             if (DUM211(K) > 0.0) then
 217                SCR319(K)=grid%xa%qcw(I,J,K)+QCT9(K)*dt(i,j,k)
 218                SCR31(K)=DUM211(K)
 219             else
 220                SCR319(K)=0.0
 221                SCR31(K)=0.0
 222             end if
 223             SCR719(K)=grid%xa%t(I,J,K)+TTT9(K)*dt(i,j,k)
 224             SCR71(K)=grid%xb%t(I,J,K)+TTT(K)*dt(i,j,k)
 225          end do
 226
 227          call da_condens_lin(DT(i,j,:),SCR31,SCR42,SCR71,DUM31,PRD,         &
 228                              QVT,QCT,QRT,TTT,                        &
 229                              grid%xb%p(I,J,:),grid%xb%t(I,J,:),grid%xb%q(I,J,:),    &
 230                              grid%xb%qcw(I,J,:),grid%xb%qrn(I,J,:),            &
 231                              SCR319,SCR429,SCR719,DUM319,PRD9,       &
 232                              QVT9,QCT9,QRT9,TTT9,                    &
 233                              grid%xa%p(I,J,:),grid%xa%t(I,J,:),grid%xa%q(I,J,:),    &
 234                              grid%xa%qcw(I,J,:),grid%xa%qrn(I,J,:),kts,kte)
 235       end do
 236    end do
 237
 238    T_NEW  (its:ite,jts:jte,kds:kde) = grid%xa%t   (its:ite,jts:jte,kds:kde) - T_OLD  (its:ite,jts:jte,kds:kde)
 239    Q_NEW  (its:ite,jts:jte,kds:kde) = grid%xa%q   (its:ite,jts:jte,kds:kde) - Q_OLD  (its:ite,jts:jte,kds:kde)
 240    QCW_NEW(its:ite,jts:jte,kds:kde) = grid%xa%qcw (its:ite,jts:jte,kds:kde) - QCW_OLD(its:ite,jts:jte,kds:kde)
 241    QRN_NEW(its:ite,jts:jte,kds:kde) = grid%xa%qrn (its:ite,jts:jte,kds:kde) - QRN_OLD(its:ite,jts:jte,kds:kde)
 242
 243    call da_filter(grid, t_new)
 244    call da_filter(grid, q_new)
 245    call da_filter(grid, qcw_new)
 246    call da_filter(grid, qrn_new)
 247
 248    grid%xa%t   (its:ite,jts:jte,kds:kde) = T_NEW  (its:ite,jts:jte,kds:kde) + T_OLD  (its:ite,jts:jte,kds:kde)
 249    grid%xa%q   (its:ite,jts:jte,kds:kde) = Q_NEW  (its:ite,jts:jte,kds:kde) + Q_OLD  (its:ite,jts:jte,kds:kde)
 250    grid%xa%qcw (its:ite,jts:jte,kds:kde) = QCW_NEW(its:ite,jts:jte,kds:kde) + QCW_OLD(its:ite,jts:jte,kds:kde)
 251    grid%xa%qrn (its:ite,jts:jte,kds:kde) = QRN_NEW(its:ite,jts:jte,kds:kde) + QRN_OLD(its:ite,jts:jte,kds:kde)
 252
 253 #ifdef DM_PARALLEL
 254 #include "HALO_XA_CLOUD.inc"
 255 #endif
 256
 257    if (trace_use) call da_trace_exit("da_moist_phys_lin")
 258
 259 end subroutine da_moist_phys_lin
 260
 261