[381] | 1 | MODULE dcmip2016_kessler_physic_mod |
---|
| 2 | |
---|
| 3 | |
---|
| 4 | |
---|
| 5 | CONTAINS |
---|
| 6 | !----------------------------------------------------------------------- |
---|
| 7 | ! |
---|
| 8 | ! Version: 2.0 |
---|
| 9 | ! |
---|
| 10 | ! Date: January 22nd, 2015 |
---|
| 11 | ! |
---|
| 12 | ! Change log: |
---|
| 13 | ! v2 - Added sub-cycling of rain sedimentation so as not to violate |
---|
| 14 | ! CFL condition. |
---|
| 15 | ! |
---|
| 16 | ! The KESSLER subroutine implements the Kessler (1969) microphysics |
---|
| 17 | ! parameterization as described by Soong and Ogura (1973) and Klemp |
---|
| 18 | ! and Wilhelmson (1978, KW). KESSLER is called at the end of each |
---|
| 19 | ! time step and makes the final adjustments to the potential |
---|
| 20 | ! temperature and moisture variables due to microphysical processes |
---|
| 21 | ! occurring during that time step. KESSLER is called once for each |
---|
| 22 | ! vertical column of grid cells. Increments are computed and added |
---|
| 23 | ! into the respective variables. The Kessler scheme contains three |
---|
| 24 | ! moisture categories: water vapor, cloud water (liquid water that |
---|
| 25 | ! moves with the flow), and rain water (liquid water that falls |
---|
| 26 | ! relative to the surrounding air). There are no ice categories. |
---|
| 27 | ! Variables in the column are ordered from the surface to the top. |
---|
| 28 | ! |
---|
| 29 | ! SUBROUTINE KESSLER(theta, qv, qc, qr, rho, pk, dt, z, nz, rainnc) |
---|
| 30 | ! |
---|
| 31 | ! Input variables: |
---|
| 32 | ! theta - potential temperature (K) |
---|
| 33 | ! qv - water vapor mixing ratio (gm/gm) |
---|
| 34 | ! qc - cloud water mixing ratio (gm/gm) |
---|
| 35 | ! qr - rain water mixing ratio (gm/gm) |
---|
| 36 | ! rho - dry air density (not mean state as in KW) (kg/m^3) |
---|
| 37 | ! pk - Exner function (not mean state as in KW) (p/p0)**(R/cp) |
---|
| 38 | ! dt - time step (s) |
---|
| 39 | ! z - heights of thermodynamic levels in the grid column (m) |
---|
| 40 | ! nz - number of thermodynamic levels in the column |
---|
[411] | 41 | ! precl - Precipitation rate (m_water/s) |
---|
[381] | 42 | ! |
---|
| 43 | ! Output variables: |
---|
| 44 | ! Increments are added into t, qv, qc, qr, and rainnc which are |
---|
| 45 | ! returned to the routine from which KESSLER was called. To obtain |
---|
| 46 | ! the total precip qt, after calling the KESSLER routine, compute: |
---|
| 47 | ! |
---|
| 48 | ! qt = sum over surface grid cells of (rainnc * cell area) (kg) |
---|
| 49 | ! [here, the conversion to kg uses (10^3 kg/m^3)*(10^-3 m/mm) = 1] |
---|
| 50 | ! |
---|
| 51 | ! |
---|
| 52 | ! Authors: Paul Ullrich |
---|
| 53 | ! University of California, Davis |
---|
| 54 | ! Email: paullrich@ucdavis.edu |
---|
| 55 | ! |
---|
| 56 | ! Based on a code by Joseph Klemp |
---|
| 57 | ! (National Center for Atmospheric Research) |
---|
| 58 | ! |
---|
| 59 | ! Reference: |
---|
| 60 | ! |
---|
| 61 | ! Klemp, J. B., W. C. Skamarock, W. C., and S.-H. Park, 2015: |
---|
| 62 | ! Idealized Global Nonhydrostatic Atmospheric Test Cases on a Reduced |
---|
| 63 | ! Radius Sphere. Journal of Advances in Modeling Earth Systems. |
---|
| 64 | ! doi:10.1002/2015MS000435 |
---|
| 65 | ! |
---|
| 66 | !======================================================================= |
---|
| 67 | |
---|
[411] | 68 | SUBROUTINE KESSLER(theta, qv, qc, qr, rho, pk, dt, z, nz, precl) |
---|
[381] | 69 | |
---|
| 70 | IMPLICIT NONE |
---|
| 71 | |
---|
| 72 | !------------------------------------------------ |
---|
| 73 | ! Input / output parameters |
---|
| 74 | !------------------------------------------------ |
---|
[901] | 75 | INTEGER, INTENT(IN) :: nz ! Number of thermodynamic levels in the column |
---|
[381] | 76 | REAL(8), DIMENSION(nz), INTENT(INOUT) :: & |
---|
| 77 | theta , & ! Potential temperature (K) |
---|
| 78 | qv , & ! Water vapor mixing ratio (gm/gm) |
---|
| 79 | qc , & ! Cloud water mixing ratio (gm/gm) |
---|
[411] | 80 | qr ! Rain water mixing ratio (gm/gm) |
---|
| 81 | |
---|
| 82 | REAL(8), DIMENSION(nz), INTENT(IN) :: & |
---|
[381] | 83 | rho ! Dry air density (not mean state as in KW) (kg/m^3) |
---|
| 84 | |
---|
[411] | 85 | REAL(8), INTENT(OUT) :: & |
---|
| 86 | precl ! Precipitation rate (m_water / s) |
---|
[381] | 87 | |
---|
| 88 | REAL(8), DIMENSION(nz), INTENT(IN) :: & |
---|
| 89 | z , & ! Heights of thermo. levels in the grid column (m) |
---|
| 90 | pk ! Exner function (p/p0)**(R/cp) |
---|
| 91 | |
---|
| 92 | REAL(8), INTENT(IN) :: & |
---|
| 93 | dt ! Time step (s) |
---|
| 94 | |
---|
| 95 | !------------------------------------------------ |
---|
| 96 | ! Local variables |
---|
| 97 | !------------------------------------------------ |
---|
| 98 | REAL, DIMENSION(nz) :: r, rhalf, velqr, sed, pc |
---|
| 99 | |
---|
| 100 | REAL(8) :: f5, f2x, xk, ern, qrprod, prod, qvs, psl, rhoqr, dt_max, dt0 |
---|
| 101 | |
---|
| 102 | INTEGER :: k, rainsplit, nt |
---|
| 103 | |
---|
| 104 | !------------------------------------------------ |
---|
| 105 | ! Begin calculation |
---|
| 106 | !------------------------------------------------ |
---|
| 107 | f2x = 17.27d0 |
---|
| 108 | f5 = 237.3d0 * f2x * 2500000.d0 / 1003.d0 |
---|
| 109 | xk = .2875d0 ! kappa (r/cp) |
---|
| 110 | psl = 1000.d0 ! pressure at sea level (mb) |
---|
| 111 | rhoqr = 1000.d0 ! density of liquid water (kg/m^3) |
---|
| 112 | |
---|
| 113 | do k=1,nz |
---|
| 114 | r(k) = 0.001d0*rho(k) |
---|
| 115 | rhalf(k) = sqrt(rho(1)/rho(k)) |
---|
| 116 | pc(k) = 3.8d0/(pk(k)**(1./xk)*psl) |
---|
| 117 | |
---|
| 118 | ! Liquid water terminal velocity (m/s) following KW eq. 2.15 |
---|
| 119 | velqr(k) = 36.34d0*(qr(k)*r(k))**0.1364*rhalf(k) |
---|
| 120 | |
---|
| 121 | end do |
---|
| 122 | |
---|
| 123 | ! Maximum time step size in accordance with CFL condition |
---|
[411] | 124 | if (dt .le. 0.d0) then |
---|
| 125 | write(*,*) 'kessler.f90 called with nonpositive dt' |
---|
| 126 | stop |
---|
| 127 | end if |
---|
| 128 | |
---|
[381] | 129 | dt_max = dt |
---|
| 130 | do k=1,nz-1 |
---|
| 131 | if (velqr(k) .ne. 0.d0) then |
---|
| 132 | dt_max = min(dt_max, 0.8d0*(z(k+1)-z(k))/velqr(k)) |
---|
| 133 | end if |
---|
| 134 | end do |
---|
| 135 | |
---|
| 136 | ! Number of subcycles |
---|
| 137 | rainsplit = ceiling(dt / dt_max) |
---|
| 138 | dt0 = dt / real(rainsplit,8) |
---|
| 139 | |
---|
| 140 | ! Subcycle through rain process |
---|
[411] | 141 | precl = 0.d0 |
---|
| 142 | |
---|
[381] | 143 | do nt=1,rainsplit |
---|
| 144 | |
---|
[411] | 145 | ! Precipitation rate (m/s) |
---|
| 146 | precl = precl + rho(1) * qr(1) * velqr(1) / rhoqr |
---|
[381] | 147 | |
---|
| 148 | ! Sedimentation term using upstream differencing |
---|
| 149 | do k=1,nz-1 |
---|
| 150 | sed(k) = dt0*(r(k+1)*qr(k+1)*velqr(k+1)-r(k)*qr(k)*velqr(k))/(r(k)*(z(k+1)-z(k))) |
---|
| 151 | end do |
---|
| 152 | sed(nz) = -dt0*qr(nz)*velqr(nz)/(.5*(z(nz)-z(nz-1))) |
---|
| 153 | |
---|
| 154 | ! Adjustment terms |
---|
| 155 | do k=1,nz |
---|
| 156 | |
---|
| 157 | ! Autoconversion and accretion rates following KW eq. 2.13a,b |
---|
| 158 | qrprod = qc(k) - (qc(k)-dt0*amax1(.001*(qc(k)-.001d0),0.))/(1.d0+dt0*2.2d0*qr(k)**.875) |
---|
| 159 | qc(k) = amax1(qc(k)-qrprod,0.) |
---|
| 160 | qr(k) = amax1(qr(k)+qrprod+sed(k),0.) |
---|
| 161 | |
---|
| 162 | ! Saturation vapor mixing ratio (gm/gm) following KW eq. 2.11 |
---|
| 163 | qvs = pc(k)*exp(f2x*(pk(k)*theta(k)-273.d0) & |
---|
| 164 | /(pk(k)*theta(k)- 36.d0)) |
---|
| 165 | prod = (qv(k)-qvs)/(1.d0+qvs*f5/(pk(k)*theta(k)-36.d0)**2) |
---|
| 166 | |
---|
| 167 | ! Evaporation rate following KW eq. 2.14a,b |
---|
| 168 | ern = amin1(dt0*(((1.6d0+124.9d0*(r(k)*qr(k))**.2046) & |
---|
| 169 | *(r(k)*qr(k))**.525)/(2550000d0*pc(k) & |
---|
| 170 | /(3.8d0 *qvs)+540000d0))*(dim(qvs,qv(k)) & |
---|
| 171 | /(r(k)*qvs)),amax1(-prod-qc(k),0.),qr(k)) |
---|
| 172 | |
---|
| 173 | ! Saturation adjustment following KW eq. 3.10 |
---|
| 174 | theta(k)= theta(k) + 2500000d0/(1003.d0*pk(k))*(amax1( prod,-qc(k))-ern) |
---|
| 175 | qv(k) = amax1(qv(k)-max(prod,-qc(k))+ern,0.) |
---|
| 176 | qc(k) = qc(k)+max(prod,-qc(k)) |
---|
| 177 | qr(k) = qr(k)-ern |
---|
| 178 | end do |
---|
| 179 | |
---|
| 180 | ! Recalculate liquid water terminal velocity |
---|
| 181 | if (nt .ne. rainsplit) then |
---|
| 182 | do k=1,nz |
---|
| 183 | velqr(k) = 36.34d0*(qr(k)*r(k))**0.1364*rhalf(k) |
---|
| 184 | end do |
---|
| 185 | end if |
---|
| 186 | end do |
---|
| 187 | |
---|
[411] | 188 | precl = precl / dble(rainsplit) |
---|
| 189 | |
---|
[381] | 190 | END SUBROUTINE KESSLER |
---|
| 191 | |
---|
| 192 | !======================================================================= |
---|
| 193 | |
---|
[411] | 194 | |
---|
| 195 | |
---|
[381] | 196 | END MODULE dcmip2016_kessler_physic_mod |
---|