| [2759] | 1 | !WRF:MODEL_LAYER:PHYSICS |
|---|
| 2 | ! |
|---|
| 3 | |
|---|
| 4 | MODULE module_mp_kessler |
|---|
| 5 | |
|---|
| 6 | CONTAINS |
|---|
| 7 | !---------------------------------------------------------------- |
|---|
| 8 | SUBROUTINE kessler( t, qv, qc, qr, rho, pii & |
|---|
| 9 | ,dt_in, z, xlv, cp & |
|---|
| 10 | ,EP2,SVP1,SVP2,SVP3,SVPT0,rhowater & |
|---|
| 11 | ,dz8w & |
|---|
| 12 | ,RAINNC, RAINNCV & |
|---|
| 13 | ,ids,ide, jds,jde, kds,kde & ! domain dims |
|---|
| 14 | ,ims,ime, jms,jme, kms,kme & ! memory dims |
|---|
| 15 | ,its,ite, jts,jte, kts,kte & ! tile dims |
|---|
| 16 | ) |
|---|
| 17 | !---------------------------------------------------------------- |
|---|
| 18 | IMPLICIT NONE |
|---|
| 19 | !---------------------------------------------------------------- |
|---|
| 20 | ! taken from the COMMAS code - WCS 10 May 1999. |
|---|
| 21 | ! converted from FORTRAN 77 to 90, tiled, WCS 10 May 1999. |
|---|
| 22 | !---------------------------------------------------------------- |
|---|
| 23 | REAL , PARAMETER :: c1 = .001 |
|---|
| 24 | REAL , PARAMETER :: c2 = .001 |
|---|
| 25 | REAL , PARAMETER :: c3 = 2.2 |
|---|
| 26 | REAL , PARAMETER :: c4 = .875 |
|---|
| 27 | REAL , PARAMETER :: fudge = 1.0 |
|---|
| 28 | REAL , PARAMETER :: mxfall = 10.0 |
|---|
| 29 | !---------------------------------------------------------------- |
|---|
| 30 | INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & |
|---|
| 31 | ims,ime, jms,jme, kms,kme, & |
|---|
| 32 | its,ite, jts,jte, kts,kte |
|---|
| 33 | REAL , INTENT(IN ) :: xlv, cp |
|---|
| 34 | REAL , INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0 |
|---|
| 35 | REAL , INTENT(IN ) :: rhowater |
|---|
| 36 | |
|---|
| 37 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 38 | INTENT(INOUT) :: & |
|---|
| 39 | t , & |
|---|
| 40 | qv, & |
|---|
| 41 | qc, & |
|---|
| 42 | qr |
|---|
| 43 | |
|---|
| 44 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 45 | INTENT(IN ) :: & |
|---|
| 46 | rho, & |
|---|
| 47 | pii, & |
|---|
| 48 | dz8w |
|---|
| 49 | |
|---|
| 50 | REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), & |
|---|
| 51 | INTENT(IN ) :: z |
|---|
| 52 | |
|---|
| 53 | REAL, INTENT(IN ) :: dt_in |
|---|
| 54 | |
|---|
| 55 | REAL, DIMENSION( ims:ime , jms:jme ), & |
|---|
| 56 | INTENT(INOUT) :: RAINNC, & |
|---|
| 57 | RAINNCV |
|---|
| 58 | |
|---|
| 59 | |
|---|
| 60 | |
|---|
| 61 | ! local variables |
|---|
| 62 | |
|---|
| 63 | REAL :: qrprod, ern, gam, rcgs, rcgsi |
|---|
| 64 | REAL, DIMENSION( its:ite , kts:kte, jts:jte ) :: prod |
|---|
| 65 | REAL, DIMENSION(kts:kte) :: vt, prodk, vtden,rdzk,rhok,factor,rdzw |
|---|
| 66 | INTEGER :: i,j,k |
|---|
| 67 | INTEGER :: nfall, n, nfall_new |
|---|
| 68 | REAL :: qrr, pressure, temp, es, qvs, dz, dt |
|---|
| 69 | REAL :: f5, dtfall, rdz, product |
|---|
| 70 | REAL :: max_heating, max_condense, max_rain, maxqrp |
|---|
| 71 | REAL :: vtmax, ernmax, crmax, factorn, time_sediment |
|---|
| 72 | REAL :: qcr, factorr, ppt |
|---|
| 73 | REAL, PARAMETER :: max_cr_sedimentation = 0.75 |
|---|
| 74 | !---------------------------------------------------------------- |
|---|
| 75 | |
|---|
| 76 | INTEGER :: imax, kmax |
|---|
| 77 | |
|---|
| 78 | dt = dt_in |
|---|
| 79 | |
|---|
| 80 | ! f5 = 237.3 * 17.27 * 2.5e6 / cp |
|---|
| 81 | f5 = svp2*(svpt0-svp3)*xlv/cp |
|---|
| 82 | ernmax = 0. |
|---|
| 83 | maxqrp = -100. |
|---|
| 84 | |
|---|
| 85 | !------------------------------------------------------------------------------ |
|---|
| 86 | ! parameters for the time split terminal advection |
|---|
| 87 | !------------------------------------------------------------------------------ |
|---|
| 88 | |
|---|
| 89 | max_heating = 0. |
|---|
| 90 | max_condense = 0. |
|---|
| 91 | max_rain = 0. |
|---|
| 92 | |
|---|
| 93 | !----------------------------------------------------------------------------- |
|---|
| 94 | ! outer J loop for entire microphysics, outer i loop for sedimentation |
|---|
| 95 | !----------------------------------------------------------------------------- |
|---|
| 96 | |
|---|
| 97 | microphysics_outer_j_loop: DO j = jts, jte |
|---|
| 98 | |
|---|
| 99 | sedimentation_outer_i_loop: DO i = its,ite |
|---|
| 100 | |
|---|
| 101 | ! vtmax = 0. |
|---|
| 102 | crmax = 0. |
|---|
| 103 | |
|---|
| 104 | |
|---|
| 105 | !------------------------------------------------------------------------------ |
|---|
| 106 | ! Terminal velocity calculation and advection, set up coefficients and |
|---|
| 107 | ! compute stable timestep |
|---|
| 108 | !------------------------------------------------------------------------------ |
|---|
| 109 | |
|---|
| 110 | DO k = 1, kte |
|---|
| 111 | prodk(k) = qr(i,k,j) |
|---|
| 112 | rhok(k) = rho(i,k,j) |
|---|
| 113 | qrr = amax1(0.,qr(i,k,j)*0.001*rhok(k)) |
|---|
| 114 | vtden(k) = sqrt(rhok(1)/rhok(k)) |
|---|
| 115 | vt(k) = 36.34*(qrr**0.1364) * vtden(k) |
|---|
| 116 | ! vtmax = amax1(vt(k), vtmax) |
|---|
| 117 | rdzw(k) = 1./dz8w(i,k,j) |
|---|
| 118 | crmax = amax1(vt(k)*dt*rdzw(k),crmax) |
|---|
| 119 | ENDDO |
|---|
| 120 | DO k = 1, kte-1 |
|---|
| 121 | rdzk(k) = 1./(z(i,k+1,j) - z(i,k,j)) |
|---|
| 122 | ENDDO |
|---|
| 123 | rdzk(kte) = 1./(z(i,kte,j) - z(i,kte-1,j)) |
|---|
| 124 | |
|---|
| 125 | nfall = max(1,nint(0.5+crmax/max_cr_sedimentation)) ! courant number for big timestep. |
|---|
| 126 | dtfall = dt / float(nfall) ! splitting so courant number for sedimentation |
|---|
| 127 | time_sediment = dt ! is stable |
|---|
| 128 | |
|---|
| 129 | !------------------------------------------------------------------------------ |
|---|
| 130 | ! Terminal velocity calculation and advection |
|---|
| 131 | ! Do a time split loop on this for stability. |
|---|
| 132 | !------------------------------------------------------------------------------ |
|---|
| 133 | |
|---|
| 134 | column_sedimentation: DO WHILE ( nfall > 0 ) |
|---|
| 135 | |
|---|
| 136 | time_sediment = time_sediment - dtfall |
|---|
| 137 | DO k = 1, kte-1 |
|---|
| 138 | factor(k) = dtfall*rdzk(k)/rhok(k) |
|---|
| 139 | ENDDO |
|---|
| 140 | factor(kte) = dtfall*rdzk(kte) |
|---|
| 141 | |
|---|
| 142 | ppt=0. |
|---|
| 143 | |
|---|
| 144 | k = 1 |
|---|
| 145 | ppt=rhok(k)*prodk(k)*vt(k)*dtfall/rhowater |
|---|
| 146 | RAINNCV(i,j)=ppt*1000. |
|---|
| 147 | RAINNC(i,j)=RAINNC(i,j)+ppt*1000. ! unit = mm |
|---|
| 148 | |
|---|
| 149 | !------------------------------------------------------------------------------ |
|---|
| 150 | ! Time split loop, Fallout done with flux upstream |
|---|
| 151 | !------------------------------------------------------------------------------ |
|---|
| 152 | |
|---|
| 153 | DO k = kts, kte-1 |
|---|
| 154 | prodk(k) = prodk(k) - factor(k) & |
|---|
| 155 | * (rhok(k)*prodk(k)*vt(k) & |
|---|
| 156 | -rhok(k+1)*prodk(k+1)*vt(k+1)) |
|---|
| 157 | ENDDO |
|---|
| 158 | |
|---|
| 159 | k = kte |
|---|
| 160 | prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k) |
|---|
| 161 | |
|---|
| 162 | !------------------------------------------------------------------------------ |
|---|
| 163 | ! compute new sedimentation velocity, and check/recompute new |
|---|
| 164 | ! sedimentation timestep if this isn't the last split step. |
|---|
| 165 | !------------------------------------------------------------------------------ |
|---|
| 166 | |
|---|
| 167 | IF( nfall > 1 ) THEN ! this wasn't the last split sedimentation timestep |
|---|
| 168 | |
|---|
| 169 | nfall = nfall - 1 |
|---|
| 170 | crmax = 0. |
|---|
| 171 | DO k = kts, kte |
|---|
| 172 | qrr = amax1(0.,prodk(k)*0.001*rhok(k)) |
|---|
| 173 | vt(k) = 36.34*(qrr**0.1364) * vtden(k) |
|---|
| 174 | ! vtmax = amax1(vt(k), vtmax) |
|---|
| 175 | crmax = amax1(vt(k)*time_sediment*rdzw(k),crmax) |
|---|
| 176 | ENDDO |
|---|
| 177 | |
|---|
| 178 | nfall_new = max(1,nint(0.5+crmax/max_cr_sedimentation)) |
|---|
| 179 | if (nfall_new /= nfall ) then |
|---|
| 180 | nfall = nfall_new |
|---|
| 181 | dtfall = time_sediment/nfall |
|---|
| 182 | end if |
|---|
| 183 | |
|---|
| 184 | ELSE ! this was the last timestep |
|---|
| 185 | |
|---|
| 186 | DO k=kts,kte |
|---|
| 187 | prod(i,k,j) = prodk(k) |
|---|
| 188 | ENDDO |
|---|
| 189 | nfall = 0 ! exit condition for sedimentation loop |
|---|
| 190 | |
|---|
| 191 | END IF |
|---|
| 192 | |
|---|
| 193 | ENDDO column_sedimentation |
|---|
| 194 | |
|---|
| 195 | ENDDO sedimentation_outer_i_loop |
|---|
| 196 | |
|---|
| 197 | !------------------------------------------------------------------------------ |
|---|
| 198 | ! Production of rain and deletion of qc |
|---|
| 199 | ! Production of qc from supersaturation |
|---|
| 200 | ! Evaporation of QR |
|---|
| 201 | !------------------------------------------------------------------------------ |
|---|
| 202 | |
|---|
| 203 | DO k = kts, kte |
|---|
| 204 | DO i = its, ite |
|---|
| 205 | factorn = 1.0 / (1.+c3*dt*amax1(0.,qr(i,k,j))**c4) |
|---|
| 206 | qrprod = qc(i,k,j) * (1.0 - factorn) & |
|---|
| 207 | + factorn*c1*dt*amax1(qc(i,k,j)-c2,0.) |
|---|
| 208 | rcgs = 0.001*rho(i,k,j) |
|---|
| 209 | |
|---|
| 210 | qc(i,k,j) = amax1(qc(i,k,j) - qrprod,0.) |
|---|
| 211 | qr(i,k,j) = (qr(i,k,j) + prod(i,k,j)-qr(i,k,j)) |
|---|
| 212 | qr(i,k,j) = amax1(qr(i,k,j) + qrprod,0.) |
|---|
| 213 | |
|---|
| 214 | temp = pii(i,k,j)*t(i,k,j) |
|---|
| 215 | pressure = 1.000e+05 * (pii(i,k,j)**(1004./287.)) |
|---|
| 216 | gam = 2.5e+06/(1004.*pii(i,k,j)) |
|---|
| 217 | ! qvs = 380.*exp(17.27*(temp-273.)/(temp- 36.))/pressure |
|---|
| 218 | es = 1000.*svp1*exp(svp2*(temp-svpt0)/(temp-svp3)) |
|---|
| 219 | qvs = ep2*es/(pressure-es) |
|---|
| 220 | ! prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+qvs*f5/(temp-36.)**2) |
|---|
| 221 | prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+pressure/(pressure-es)*qvs*f5/(temp-svp3)**2) |
|---|
| 222 | ern = amin1(dt*(((1.6+124.9*(rcgs*qr(i,k,j))**.2046) & |
|---|
| 223 | *(rcgs*qr(i,k,j))**.525)/(2.55e8/(pressure*qvs) & |
|---|
| 224 | +5.4e5))*(dim(qvs,qv(i,k,j))/(rcgs*qvs)), & |
|---|
| 225 | amax1(-prod(i,k,j)-qc(i,k,j),0.),qr(i,k,j)) |
|---|
| 226 | |
|---|
| 227 | ! Update all variables |
|---|
| 228 | |
|---|
| 229 | product = amax1(prod(i,k,j),-qc(i,k,j)) |
|---|
| 230 | t (i,k,j) = t(i,k,j) + gam*(product - ern) |
|---|
| 231 | qv(i,k,j) = amax1(qv(i,k,j) - product + ern,0.) |
|---|
| 232 | qc(i,k,j) = qc(i,k,j) + product |
|---|
| 233 | qr(i,k,j) = qr(i,k,j) - ern |
|---|
| 234 | |
|---|
| 235 | ENDDO |
|---|
| 236 | ENDDO |
|---|
| 237 | |
|---|
| 238 | ENDDO microphysics_outer_j_loop |
|---|
| 239 | |
|---|
| 240 | RETURN |
|---|
| 241 | |
|---|
| 242 | END SUBROUTINE kessler |
|---|
| 243 | |
|---|
| 244 | END MODULE module_mp_kessler |
|---|
