PageRenderTime 56ms CodeModel.GetById 10ms RepoModel.GetById 1ms app.codeStats 0ms

/wrfv2_fire/phys/module_mp_wsm3_accel.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1683 lines | 1205 code | 25 blank | 453 comment | 91 complexity | 6e5ab2c544fe92863a73d707643372bd MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

  1. #if ( RWORDSIZE == 4 )
  2. # define VREC vsrec
  3. # define VSQRT vssqrt
  4. #else
  5. # define VREC vrec
  6. # define VSQRT vsqrt
  7. #endif
  8. MODULE module_mp_wsm3
  9. !
  10. !
  11. REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
  12. REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
  13. REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
  14. REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
  15. REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
  16. REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
  17. REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
  18. REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
  19. REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
  20. REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
  21. REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
  22. REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain
  23. REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
  24. REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
  25. REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
  26. REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
  27. REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
  28. REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
  29. REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
  30. REAL, SAVE :: &
  31. qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
  32. g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
  33. precr1,precr2,xmmax,roqimax,bvts1, &
  34. bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
  35. g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
  36. pidn0s,xlv1, &
  37. rslopermax,rslopesmax,rslopegmax, &
  38. rsloperbmax,rslopesbmax,rslopegbmax, &
  39. rsloper2max,rslopes2max,rslopeg2max, &
  40. rsloper3max,rslopes3max,rslopeg3max
  41. !
  42. ! Specifies code-inlining of fpvs function in WSM32D below. JM 20040507
  43. !
  44. CONTAINS
  45. !===================================================================
  46. !
  47. SUBROUTINE wsm3(th, q, qci, qrs &
  48. , w, den, pii, p, delz &
  49. , delt,g, cpd, cpv, rd, rv, t0c &
  50. , ep1, ep2, qmin &
  51. , XLS, XLV0, XLF0, den0, denr &
  52. , cliq,cice,psat &
  53. , rain, rainncv &
  54. , snow, snowncv &
  55. , sr &
  56. , ids,ide, jds,jde, kds,kde &
  57. , ims,ime, jms,jme, kms,kme &
  58. , its,ite, jts,jte, kts,kte &
  59. )
  60. !-------------------------------------------------------------------
  61. #ifdef _OPENMP
  62. use omp_lib
  63. #endif
  64. IMPLICIT NONE
  65. !-------------------------------------------------------------------
  66. !
  67. !
  68. ! This code is a 3-class simple ice microphyiscs scheme (WSM3) of the WRF
  69. ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
  70. ! number concentration is a function of temperature, and seperate assumption
  71. ! is developed, in which ice crystal number concentration is a function
  72. ! of ice amount. A theoretical background of the ice-microphysics and related
  73. ! processes in the WSMMPs are described in Hong et al. (2004).
  74. ! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
  75. ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
  76. !
  77. ! WSM3 cloud scheme
  78. !
  79. ! Coded by Song-You Hong (Yonsei Univ.)
  80. ! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
  81. ! Summer 2002
  82. !
  83. ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
  84. ! Summer 2003
  85. !
  86. ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
  87. ! Dudhia (D89, 1989) J. Atmos. Sci.
  88. ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
  89. !
  90. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
  91. ims,ime, jms,jme, kms,kme , &
  92. its,ite, jts,jte, kts,kte
  93. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
  94. INTENT(INOUT) :: &
  95. th, &
  96. q, &
  97. qci, &
  98. qrs
  99. REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
  100. INTENT(IN ) :: w, &
  101. den, &
  102. pii, &
  103. p, &
  104. delz
  105. REAL, INTENT(IN ) :: delt, &
  106. g, &
  107. rd, &
  108. rv, &
  109. t0c, &
  110. den0, &
  111. cpd, &
  112. cpv, &
  113. ep1, &
  114. ep2, &
  115. qmin, &
  116. XLS, &
  117. XLV0, &
  118. XLF0, &
  119. cliq, &
  120. cice, &
  121. psat, &
  122. denr
  123. REAL, DIMENSION( ims:ime , jms:jme ), &
  124. INTENT(INOUT) :: rain, &
  125. rainncv
  126. REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
  127. INTENT(INOUT) :: snow, &
  128. snowncv, &
  129. sr
  130. ! LOCAL VAR
  131. REAL, DIMENSION( its:ite , kts:kte ) :: t
  132. INTEGER :: i,j,k
  133. #ifdef _ACCEL_PROF
  134. integer :: l
  135. real*8 wsm3_t(8,256), wsm5_t(8,256), t1, t2
  136. common /wsm_times/ wsm3_t(8,256), wsm5_t(8,256)
  137. #endif
  138. !-------------------------------------------------------------------
  139. #ifdef _ACCEL_PROF
  140. call cpu_time(t1)
  141. #endif
  142. #ifdef _ACCEL
  143. CALL wsm32D(th, q, qci, qrs, &
  144. w, den, pii, p, delz, rain, rainncv, &
  145. delt,g, cpd, cpv, rd, rv, t0c, &
  146. ep1, ep2, qmin, &
  147. XLS, XLV0, XLF0, den0, denr, &
  148. cliq,cice,psat, &
  149. ids,ide, jds,jde, kds,kde, &
  150. ims,ime, jms,jme, kms,kme, &
  151. its,ite, jts,jte, kts,kte )
  152. #else
  153. DO j=jts,jte
  154. DO k=kts,kte
  155. DO i=its,ite
  156. t(i,k)=th(i,k,j)*pii(i,k,j)
  157. ENDDO
  158. ENDDO
  159. CALL wsm32D(t, q(ims,kms,j), qci(ims,kms,j) &
  160. ,qrs(ims,kms,j),w(ims,kms,j), den(ims,kms,j) &
  161. ,p(ims,kms,j), delz(ims,kms,j) &
  162. ,delt,g, cpd, cpv, rd, rv, t0c &
  163. ,ep1, ep2, qmin &
  164. ,XLS, XLV0, XLF0, den0, denr &
  165. ,cliq,cice,psat &
  166. ,j &
  167. ,rain(ims,j), rainncv(ims,j) &
  168. ,snow(ims,j),snowncv(ims,j) &
  169. ,sr(ims,j) &
  170. ,ids,ide, jds,jde, kds,kde &
  171. ,ims,ime, jms,jme, kms,kme &
  172. ,its,ite, jts,jte, kts,kte &
  173. )
  174. DO K=kts,kte
  175. DO I=its,ite
  176. th(i,k,j)=t(i,k)/pii(i,k,j)
  177. ENDDO
  178. ENDDO
  179. ENDDO
  180. #endif
  181. #ifdef _ACCEL_PROF
  182. call cpu_time(t2)
  183. #ifdef _OPENMP
  184. l = omp_get_thread_num() + 1
  185. #else
  186. l = 1
  187. #endif
  188. wsm3_t(1,l) = wsm3_t(1,l) + (t2 - t1)
  189. #endif
  190. END SUBROUTINE wsm3
  191. #ifdef _ACCEL
  192. !===================================================================
  193. !{
  194. SUBROUTINE wsm32D(th, q, qci, qrs, &
  195. w, den, pii, p, delz, rain, rainncv, &
  196. delt,g, cpd, cpv, rd, rv, t0c, &
  197. ep1, ep2, qmin, &
  198. XLS, XLV0, XLF0, den0, denr, &
  199. cliq,cice,psat, &
  200. ids,ide, jds,jde, kds,kde, &
  201. ims,ime, jms,jme, kms,kme, &
  202. its,ite, jts,jte, kts,kte )
  203. !-------------------------------------------------------------------
  204. IMPLICIT NONE
  205. !-------------------------------------------------------------------
  206. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
  207. ims,ime, jms,jme, kms,kme , &
  208. its,ite, jts,jte, kts,kte
  209. REAL, DIMENSION( ims:ime , kms:kme, ims:ims ), &
  210. INTENT(INOUT) :: &
  211. th
  212. REAL, DIMENSION( ims:ime , kms:kme, ims:ims ), &
  213. INTENT(IN) :: &
  214. pii
  215. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  216. INTENT(INOUT) :: &
  217. q, &
  218. qci, &
  219. qrs
  220. REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
  221. INTENT(IN ) :: w, &
  222. den, &
  223. p, &
  224. delz
  225. REAL, DIMENSION( ims:ime , jms:jme ), &
  226. INTENT(INOUT) :: rain, &
  227. rainncv
  228. REAL, INTENT(IN ) :: delt, &
  229. g, &
  230. cpd, &
  231. cpv, &
  232. t0c, &
  233. den0, &
  234. rd, &
  235. rv, &
  236. ep1, &
  237. ep2, &
  238. qmin, &
  239. XLS, &
  240. XLV0, &
  241. XLF0, &
  242. cliq, &
  243. cice, &
  244. psat, &
  245. denr
  246. ! LOCAL VAR
  247. REAL, DIMENSION( its:ite , kts:kte ) :: &
  248. rh, qs, denfac, rslope, rslope2, rslope3, rslopeb, &
  249. pgen, paut, pacr, pisd, pres, pcon, fall, falk, &
  250. xl, cpm, work1, work2, xni, qs0, n0sfac
  251. ! LOCAL VAR
  252. REAL, DIMENSION( its:ite , kts:kte, jts:jte ) :: t
  253. REAL, DIMENSION( its:ite , kts:kte ) :: &
  254. falkc, work1c, work2c, fallc
  255. INTEGER :: mstep, numdt
  256. LOGICAL, DIMENSION( its:ite ) :: flgcld
  257. REAL :: pi, &
  258. cpmcal, xlcal, lamdar, lamdas, diffus, &
  259. viscos, xka, venfac, conden, diffac, &
  260. x, y, z, a, b, c, d, e, &
  261. qdt, pvt, qik, delq, facq, qrsci, frzmlt, &
  262. snomlt, hold, holdrs, facqci, supcol, coeres, &
  263. supsat, dtcld, xmi, qciik, delqci, eacrs, satdt, &
  264. qimax, diameter, xni0, roqi0
  265. REAL :: holdc, holdci
  266. INTEGER :: i, k, j, &
  267. iprt, latd, lond, loop, loops, ifsat, kk, n
  268. !
  269. #define INL
  270. #ifdef INL
  271. ! Temporaries used for inlining fpvs function
  272. REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
  273. #endif
  274. !=================================================================
  275. ! compute internal functions
  276. !
  277. cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
  278. xlcal(x) = xlv0-xlv1*(x-t0c)
  279. ! tvcal(x,y) = x+x*ep1*max(y,qmin)
  280. !----------------------------------------------------------------
  281. ! size distributions: (x=mixing ratio, y=air density):
  282. ! valid for mixing ratio > 1.e-9 kg/kg.
  283. !
  284. lamdar(x,y)=(pidn0r/(x*y))**.25
  285. lamdas(x,y,z)=(pidn0s*z/(x*y))**.25
  286. !
  287. !----------------------------------------------------------------
  288. ! diffus: diffusion coefficient of the water vapor
  289. ! viscos: kinematic viscosity(m2s-1)
  290. !
  291. diffus(x,y) = 8.794e-5*x**1.81/y
  292. viscos(x,y) = 1.496e-6*x**1.5/(x+120.)/y
  293. xka(x,y) = 1.414e3*viscos(x,y)*y
  294. diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
  295. venfac(a,b,c) = (viscos(b,c)/diffus(b,a))**(.3333333) &
  296. /viscos(b,c)**(.5)*(den0/c)**0.25
  297. conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
  298. !
  299. pi = 4. * atan(1.)
  300. !
  301. !----------------------------------------------------------------
  302. ! compute the minor time steps.
  303. !
  304. loops = max(int(delt/dtcldcr+0.5),1)
  305. dtcld = delt/loops
  306. if(delt.le.dtcldcr) dtcld = delt
  307. #ifdef INL
  308. cvap = cpv
  309. hvap=xlv0
  310. hsub=xls
  311. ttp=t0c+0.1
  312. dldt=cvap-cliq
  313. xa=-dldt/rv
  314. xb=xa+hvap/(rv*ttp)
  315. dldti=cvap-cice
  316. xai=-dldti/rv
  317. xbi=xai+hsub/(rv*ttp)
  318. #endif
  319. !
  320. !----------------------------------------------------------------
  321. ! paddint 0 for negative values generated by dynamics
  322. !
  323. !$acc region &
  324. !$acc local(t) &
  325. !$acc copyin(delz(:,:,:),p(:,:,:)) &
  326. !$acc copyin(den(:,:,:),w(:,:,:)) &
  327. !$acc copy(q(:,:,:),qci(:,:,:),qrs(:,:,:))
  328. !$acc do &
  329. !$acc private(rh,qs,denfac,rslope,rslope2,rslope3,rslopeb) &
  330. !$acc private(pgen,paut,pacr,pisd,pres,pcon,fall,falk) &
  331. !$acc private(xl,cpm,work1,work2,xni,qs0,n0sfac) &
  332. !$acc private(falkc,work1c,work2c,fallc) &
  333. !$acc parallel
  334. do j = jts, jte
  335. !$acc do &
  336. !$acc private(numdt,mstep) &
  337. !$acc kernel vector(96)
  338. do i = its, ite
  339. do k = kts, kte
  340. t(i,k,j)=th(i,k,j)*pii(i,k,j)
  341. qci(i,k,j) = max(qci(i,k,j),0.0)
  342. qrs(i,k,j) = max(qrs(i,k,j),0.0)
  343. enddo
  344. !
  345. !----------------------------------------------------------------
  346. ! latent heat for phase changes and heat capacity. neglect the
  347. ! changes during microphysical process calculation
  348. ! emanuel(1994)
  349. !
  350. do k = kts, kte
  351. cpm(i,k) = cpmcal(q(i,k,j))
  352. xl(i,k) = xlcal(t(i,k,j))
  353. enddo
  354. !
  355. do loop = 1,loops
  356. !
  357. !----------------------------------------------------------------
  358. ! initialize the large scale variables
  359. !
  360. mstep = 1
  361. flgcld(i) = .true.
  362. !
  363. do k = kts, kte
  364. denfac(i,k) = sqrt(den0/den(i,k,j))
  365. enddo
  366. !
  367. do k = kts, kte
  368. #ifndef INL
  369. qs(i,k) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
  370. qs0(i,k) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
  371. #else
  372. tr=ttp/t(i,k,j)
  373. if(t(i,k,j).lt.ttp) then
  374. qs(i,k) =psat*(tr**xai)*exp(xbi*(1.-tr))
  375. else
  376. qs(i,k) =psat*(tr**xa)*exp(xb*(1.-tr))
  377. endif
  378. qs0(i,k) =psat*(tr**xa)*exp(xb*(1.-tr))
  379. #endif
  380. qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k)
  381. qs(i,k) = ep2 * qs(i,k) / (p(i,k,j) - qs(i,k))
  382. qs(i,k) = max(qs(i,k),qmin)
  383. rh(i,k) = max(q(i,k,j) / qs(i,k),qmin)
  384. enddo
  385. !
  386. !----------------------------------------------------------------
  387. ! initialize the variables for microphysical physics
  388. !
  389. !
  390. do k = kts, kte
  391. pres(i,k) = 0.
  392. paut(i,k) = 0.
  393. pacr(i,k) = 0.
  394. pgen(i,k) = 0.
  395. pisd(i,k) = 0.
  396. pcon(i,k) = 0.
  397. fall(i,k) = 0.
  398. falk(i,k) = 0.
  399. fallc(i,k) = 0.
  400. falkc(i,k) = 0.
  401. xni(i,k) = 1.e3
  402. enddo
  403. !
  404. !----------------------------------------------------------------
  405. ! compute the fallout term:
  406. ! first, vertical terminal velosity for minor loops
  407. !---------------------------------------------------------------
  408. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
  409. !---------------------------------------------------------------
  410. do k = kts, kte
  411. supcol = t0c-t(i,k,j)
  412. n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
  413. if(t(i,k,j).ge.t0c) then
  414. if(qrs(i,k,j).le.qcrmin)then
  415. rslope(i,k) = rslopermax
  416. rslopeb(i,k) = rsloperbmax
  417. rslope2(i,k) = rsloper2max
  418. rslope3(i,k) = rsloper3max
  419. else
  420. rslope(i,k) = 1./lamdar(qrs(i,k,j),den(i,k,j))
  421. rslopeb(i,k) = rslope(i,k)**bvtr
  422. rslope2(i,k) = rslope(i,k)*rslope(i,k)
  423. rslope3(i,k) = rslope2(i,k)*rslope(i,k)
  424. endif
  425. else
  426. if(qrs(i,k,j).le.qcrmin)then
  427. rslope(i,k) = rslopesmax
  428. rslopeb(i,k) = rslopesbmax
  429. rslope2(i,k) = rslopes2max
  430. rslope3(i,k) = rslopes3max
  431. else
  432. rslope(i,k) = 1./lamdas(qrs(i,k,j),den(i,k,j),n0sfac(i,k))
  433. rslopeb(i,k) = rslope(i,k)**bvts
  434. rslope2(i,k) = rslope(i,k)*rslope(i,k)
  435. rslope3(i,k) = rslope2(i,k)*rslope(i,k)
  436. endif
  437. endif
  438. !-------------------------------------------------------------
  439. ! Ni: ice crystal number concentraiton [HDC 5c]
  440. !-------------------------------------------------------------
  441. xni(i,k) = min(max(5.38e7*(den(i,k,j) &
  442. *max(qci(i,k,j),qmin))**0.75,1.e3),1.e6)
  443. enddo
  444. !
  445. numdt = 1
  446. do k = kte, kts, -1
  447. if(t(i,k,j).lt.t0c) then
  448. pvt = pvts
  449. else
  450. pvt = pvtr
  451. endif
  452. work1(i,k) = pvt*rslopeb(i,k)*denfac(i,k)
  453. work2(i,k) = work1(i,k)/delz(i,k,j)
  454. numdt = max(int(work2(i,k)*dtcld+1.),1)
  455. if(numdt.ge.mstep) mstep = numdt
  456. enddo
  457. !
  458. do n = 1, mstep
  459. k = kte
  460. falk(i,k) = den(i,k,j)*qrs(i,k,j)*work2(i,k)/mstep
  461. hold = falk(i,k)
  462. fall(i,k) = fall(i,k)+falk(i,k)
  463. holdrs = qrs(i,k,j)
  464. qrs(i,k,j) = max(qrs(i,k,j)-falk(i,k)*dtcld/den(i,k,j),0.)
  465. do k = kte-1, kts, -1
  466. falk(i,k) = den(i,k,j)*qrs(i,k,j)*work2(i,k)/mstep
  467. hold = falk(i,k)
  468. fall(i,k) = fall(i,k)+falk(i,k)
  469. holdrs = qrs(i,k,j)
  470. qrs(i,k,j) = max(qrs(i,k,j)-(falk(i,k) &
  471. -falk(i,k+1)*delz(i,k+1,j)/delz(i,k,j))*dtcld/den(i,k,j),0.)
  472. enddo
  473. enddo
  474. !---------------------------------------------------------------
  475. ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
  476. !---------------------------------------------------------------
  477. mstep = 1
  478. numdt = 1
  479. do k = kte, kts, -1
  480. if(t(i,k,j).lt.t0c.and.qci(i,k,j).gt.0.) then
  481. xmi = den(i,k,j)*qci(i,k,j)/xni(i,k)
  482. diameter = dicon * sqrt(xmi)
  483. work1c(i,k) = 1.49e4*diameter**1.31
  484. else
  485. work1c(i,k) = 0.
  486. endif
  487. if(qci(i,k,j).le.0.) then
  488. work2c(i,k) = 0.
  489. else
  490. work2c(i,k) = work1c(i,k)/delz(i,k,j)
  491. endif
  492. numdt = max(int(work2c(i,k)*dtcld+1.),1)
  493. if(numdt.ge.mstep) mstep = numdt
  494. enddo
  495. !
  496. do n = 1, mstep
  497. k = kte
  498. falkc(i,k) = den(i,k,j)*qci(i,k,j)*work2c(i,k)/mstep
  499. holdc = falkc(i,k)
  500. fallc(i,k) = fallc(i,k)+falkc(i,k)
  501. holdci = qci(i,k,j)
  502. qci(i,k,j) = max(qci(i,k,j)-falkc(i,k)*dtcld/den(i,k,j),0.)
  503. do k = kte-1, kts, -1
  504. falkc(i,k) = den(i,k,j)*qci(i,k,j)*work2c(i,k)/mstep
  505. holdc = falkc(i,k)
  506. fallc(i,k) = fallc(i,k)+falkc(i,k)
  507. holdci = qci(i,k,j)
  508. qci(i,k,j) = max(qci(i,k,j)-(falkc(i,k) &
  509. -falkc(i,k+1)*delz(i,k+1,j)/delz(i,k,j))*dtcld/den(i,k,j),0.)
  510. enddo
  511. enddo
  512. !
  513. !----------------------------------------------------------------
  514. ! compute the freezing/melting term. [D89 B16-B17]
  515. ! freezing occurs one layer above the melting level
  516. !
  517. mstep = 0
  518. !
  519. do k = kts, kte
  520. if(t(i,k,j).ge.t0c) then
  521. mstep = k
  522. endif
  523. enddo
  524. !
  525. if(mstep.ne.0.and.w(i,mstep,j).gt.0.) then
  526. work1(i,1) = float(mstep + 1)
  527. work1(i,2) = float(mstep)
  528. else
  529. work1(i,1) = float(mstep)
  530. work1(i,2) = float(mstep)
  531. endif
  532. !
  533. k = int(work1(i,1)+0.5)
  534. kk = int(work1(i,2)+0.5)
  535. if(k*kk.ge.1) then
  536. qrsci = qrs(i,k,j) + qci(i,k,j)
  537. if(qrsci.gt.0..or.fall(i,kk).gt.0.) then
  538. frzmlt = min(max(-w(i,k,j)*qrsci/delz(i,k,j),-qrsci/dtcld), &
  539. qrsci/dtcld)
  540. snomlt = min(max(fall(i,kk)/den(i,kk,j),-qrs(i,k,j)/dtcld), &
  541. qrs(i,k,j)/dtcld)
  542. if(k.eq.kk) then
  543. t(i,k,j) = t(i,k,j) - xlf0/cpm(i,k)*(frzmlt+snomlt)*dtcld
  544. else
  545. t(i,k,j) = t(i,k,j) - xlf0/cpm(i,k)*frzmlt*dtcld
  546. t(i,kk,j) = t(i,kk,j) - xlf0/cpm(i,kk)*snomlt*dtcld
  547. endif
  548. endif
  549. endif
  550. !
  551. !----------------------------------------------------------------
  552. ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
  553. !
  554. if(fall(i,1).gt.0.) then
  555. rainncv(i,j) = fall(i,1)*delz(i,1,j)/denr*dtcld*1000.
  556. rain(i,j) = fall(i,1)*delz(i,1,j)/denr*dtcld*1000. &
  557. + rain(i,j)
  558. endif
  559. !
  560. !----------------------------------------------------------------
  561. ! rsloper: reverse of the slope parameter of the rain(m,j)
  562. ! xka: thermal conductivity of air(jm-1s-1k-1)
  563. ! work1: the thermodynamic term in the denominator associated with
  564. ! heat conduction and vapor diffusion
  565. ! (ry88, y93, h85)
  566. ! work2: parameter associated with the ventilation effects(y93)
  567. !
  568. do k = kts, kte
  569. if(t(i,k,j).ge.t0c) then
  570. if(qrs(i,k,j).le.qcrmin)then
  571. rslope(i,k) = rslopermax
  572. rslopeb(i,k) = rsloperbmax
  573. rslope2(i,k) = rsloper2max
  574. rslope3(i,k) = rsloper3max
  575. else
  576. rslope(i,k) = 1./lamdar(qrs(i,k,j),den(i,k,j))
  577. rslopeb(i,k) = rslope(i,k)**bvtr
  578. rslope2(i,k) = rslope(i,k)*rslope(i,k)
  579. rslope3(i,k) = rslope2(i,k)*rslope(i,k)
  580. endif
  581. else
  582. if(qrs(i,k,j).le.qcrmin)then
  583. rslope(i,k) = rslopesmax
  584. rslopeb(i,k) = rslopesbmax
  585. rslope2(i,k) = rslopes2max
  586. rslope3(i,k) = rslopes3max
  587. else
  588. rslope(i,k) = 1./lamdas(qrs(i,k,j),den(i,k,j),n0sfac(i,k))
  589. rslopeb(i,k) = rslope(i,k)**bvts
  590. rslope2(i,k) = rslope(i,k)*rslope(i,k)
  591. rslope3(i,k) = rslope2(i,k)*rslope(i,k)
  592. endif
  593. endif
  594. enddo
  595. !
  596. do k = kts, kte
  597. if(t(i,k,j).ge.t0c) then
  598. work1(i,k) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k))
  599. else
  600. work1(i,k) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k))
  601. endif
  602. work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
  603. enddo
  604. !
  605. do k = kts, kte
  606. supsat = max(q(i,k,j),qmin)-qs(i,k)
  607. satdt = supsat/dtcld
  608. if(t(i,k,j).ge.t0c) then
  609. !
  610. !===============================================================
  611. !
  612. ! warm rain processes
  613. !
  614. ! - follows the processes in RH83 and LFO except for autoconcersion
  615. !
  616. !===============================================================
  617. !---------------------------------------------------------------
  618. ! paut1: auto conversion rate from cloud to rain [HDC 16]
  619. ! (C->R)
  620. !---------------------------------------------------------------
  621. if(qci(i,k,j).gt.qc0) then
  622. paut(i,k) = qck1*qci(i,k,j)**(7./3.)
  623. paut(i,k) = min(paut(i,k),qci(i,k,j)/dtcld)
  624. endif
  625. !---------------------------------------------------------------
  626. ! pracw: accretion of cloud water by rain [D89 B15]
  627. ! (C->R)
  628. !---------------------------------------------------------------
  629. if(qrs(i,k,j).gt.qcrmin.and.qci(i,k,j).gt.qmin) then
  630. pacr(i,k) = min(pacrr*rslope3(i,k)*rslopeb(i,k) &
  631. *qci(i,k,j)*denfac(i,k),qci(i,k,j)/dtcld)
  632. endif
  633. !---------------------------------------------------------------
  634. ! pres1: evaporation/condensation rate of rain [HDC 14]
  635. ! (V->R or R->V)
  636. !---------------------------------------------------------------
  637. if(qrs(i,k,j).gt.0.) then
  638. coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
  639. pres(i,k) = (rh(i,k)-1.)*(precr1*rslope2(i,k) &
  640. +precr2*work2(i,k)*coeres)/work1(i,k)
  641. if(pres(i,k).lt.0.) then
  642. pres(i,k) = max(pres(i,k),-qrs(i,k,j)/dtcld)
  643. pres(i,k) = max(pres(i,k),satdt/2)
  644. else
  645. pres(i,k) = min(pres(i,k),satdt/2)
  646. endif
  647. endif
  648. else
  649. !
  650. !===============================================================
  651. !
  652. ! cold rain processes
  653. !
  654. ! - follows the revised ice microphysics processes in HDC
  655. ! - the processes same as in RH83 and LFO behave
  656. ! following ice crystal hapits defined in HDC, inclduing
  657. ! intercept parameter for snow (n0s), ice crystal number
  658. ! concentration (ni), ice nuclei number concentration
  659. ! (n0i), ice diameter (d)
  660. !
  661. !===============================================================
  662. !
  663. supcol = t0c-t(i,k,j)
  664. ifsat = 0
  665. !-------------------------------------------------------------
  666. ! Ni: ice crystal number concentraiton [HDC 5c]
  667. !-------------------------------------------------------------
  668. xni(i,k) = min(max(5.38e7*(den(i,k,j) &
  669. *max(qci(i,k,j),qmin))**0.75,1.e3),1.e6)
  670. eacrs = exp(0.05*(-supcol))
  671. !
  672. if(qrs(i,k,j).gt.qcrmin.and.qci(i,k,j).gt.qmin) then
  673. pacr(i,k) = min(pacrs*n0sfac(i,k)*eacrs*rslope3(i,k) &
  674. *rslopeb(i,k)*qci(i,k,j)*denfac(i,k),qci(i,k,j)/dtcld)
  675. endif
  676. !-------------------------------------------------------------
  677. ! pisd: Deposition/Sublimation rate of ice [HDC 9]
  678. ! (T<T0: V->I or I->V)
  679. !-------------------------------------------------------------
  680. if(qci(i,k,j).gt.0.) then
  681. xmi = den(i,k,j)*qci(i,k,j)/xni(i,k)
  682. diameter = dicon * sqrt(xmi)
  683. pisd(i,k) = 4.*diameter*xni(i,k)*(rh(i,k)-1.) &
  684. /work1(i,k)
  685. if(pisd(i,k).lt.0.) then
  686. pisd(i,k) = max(pisd(i,k),satdt/2)
  687. pisd(i,k) = max(pisd(i,k),-qci(i,k,j)/dtcld)
  688. else
  689. pisd(i,k) = min(pisd(i,k),satdt/2)
  690. endif
  691. if(abs(pisd(i,k)).ge.abs(satdt)) ifsat = 1
  692. endif
  693. !-------------------------------------------------------------
  694. ! pres2: deposition/sublimation rate of snow [HDC 14]
  695. ! (V->S or S->V)
  696. !-------------------------------------------------------------
  697. if(qrs(i,k,j).gt.0..and.ifsat.ne.1) then
  698. coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
  699. pres(i,k) = (rh(i,k)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k) &
  700. +precs2*work2(i,k)*coeres)/work1(i,k)
  701. if(pres(i,k).lt.0.) then
  702. pres(i,k) = max(pres(i,k),-qrs(i,k,j)/dtcld)
  703. pres(i,k) = max(pres(i,k),satdt/2)
  704. else
  705. pres(i,k) = min(pres(i,k),satdt/2)
  706. endif
  707. if(abs(pisd(i,k)+pres(i,k)).ge.abs(satdt)) ifsat = 1
  708. endif
  709. !-------------------------------------------------------------
  710. ! pgen: generation(nucleation) of ice from vapor [HDC 7-8]
  711. ! (T<T0: V->I)
  712. !-------------------------------------------------------------
  713. if(supsat.gt.0.and.ifsat.ne.1) then
  714. xni0 = 1.e3*exp(0.1*supcol)
  715. roqi0 = 4.92e-11*xni0**1.33
  716. pgen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,j),0.))/dtcld)
  717. pgen(i,k) = min(pgen(i,k),satdt)
  718. endif
  719. !-------------------------------------------------------------
  720. ! paut2: conversion(aggregation) of ice to snow [HDC 12]
  721. ! (T<T0: I->S)
  722. !-------------------------------------------------------------
  723. if(qci(i,k,j).gt.0.) then
  724. qimax = roqimax/den(i,k,j)
  725. paut(i,k) = max(0.,(qci(i,k,j)-qimax)/dtcld)
  726. endif
  727. endif
  728. enddo
  729. !
  730. !----------------------------------------------------------------
  731. ! check mass conservation of generation terms and feedback to the
  732. ! large scale
  733. !
  734. do k = kts, kte
  735. qciik = max(qmin,qci(i,k,j))
  736. delqci = (paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))*dtcld
  737. if(delqci.ge.qciik) then
  738. facqci = qciik/delqci
  739. paut(i,k) = paut(i,k)*facqci
  740. pacr(i,k) = pacr(i,k)*facqci
  741. pgen(i,k) = pgen(i,k)*facqci
  742. pisd(i,k) = pisd(i,k)*facqci
  743. endif
  744. qik = max(qmin,q(i,k,j))
  745. delq = (pres(i,k)+pgen(i,k)+pisd(i,k))*dtcld
  746. if(delq.ge.qik) then
  747. facq = qik/delq
  748. pres(i,k) = pres(i,k)*facq
  749. pgen(i,k) = pgen(i,k)*facq
  750. pisd(i,k) = pisd(i,k)*facq
  751. endif
  752. work2(i,k) = -pres(i,k)-pgen(i,k)-pisd(i,k)
  753. q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
  754. qci(i,k,j) = max(qci(i,k,j)-(paut(i,k)+pacr(i,k)-pgen(i,k) &
  755. -pisd(i,k))*dtcld,0.)
  756. qrs(i,k,j) = max(qrs(i,k,j)+(paut(i,k)+pacr(i,k) &
  757. +pres(i,k))*dtcld,0.)
  758. if(t(i,k,j).lt.t0c) then
  759. t(i,k,j) = t(i,k,j)-xls*work2(i,k)/cpm(i,k)*dtcld
  760. else
  761. t(i,k,j) = t(i,k,j)-xl(i,k)*work2(i,k)/cpm(i,k)*dtcld
  762. endif
  763. enddo
  764. !
  765. do k = kts, kte
  766. #ifndef INL
  767. qs(i,k) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
  768. #else
  769. tr=ttp/t(i,k,j)
  770. qs(i,k)=psat*(tr**xa)*exp(xb*(1.-tr))
  771. #endif
  772. qs(i,k) = ep2 * qs(i,k) / (p(i,k,j) - qs(i,k))
  773. qs(i,k) = max(qs(i,k),qmin)
  774. denfac(i,k) = sqrt(den0/den(i,k,j))
  775. enddo
  776. !
  777. !----------------------------------------------------------------
  778. ! pcon: condensational/evaporational rate of cloud water [RH83 A6]
  779. ! if there exists additional water vapor condensated/if
  780. ! evaporation of cloud water is not enough to remove subsaturation
  781. !
  782. do k = kts, kte
  783. work1(i,k) = conden(t(i,k,j),q(i,k,j),qs(i,k),xl(i,k),cpm(i,k))
  784. work2(i,k) = qci(i,k,j)+work1(i,k)
  785. pcon(i,k) = min(max(work1(i,k),0.),max(q(i,k,j),0.))/dtcld
  786. if(qci(i,k,j).gt.0..and.work1(i,k).lt.0.and.t(i,k,j).gt.t0c) &
  787. pcon(i,k) = max(work1(i,k),-qci(i,k,j))/dtcld
  788. q(i,k,j) = q(i,k,j)-pcon(i,k)*dtcld
  789. qci(i,k,j) = max(qci(i,k,j)+pcon(i,k)*dtcld,0.)
  790. t(i,k,j) = t(i,k,j)+pcon(i,k)*xl(i,k)/cpm(i,k)*dtcld
  791. enddo
  792. !
  793. !----------------------------------------------------------------
  794. ! padding for small values
  795. !
  796. do k = kts, kte
  797. if(qci(i,k,j).le.qmin) qci(i,k,j) = 0.0
  798. enddo
  799. !
  800. enddo ! big loops
  801. DO K=kts,kte
  802. th(i,k,j)=t(i,k,j)/pii(i,k,j)
  803. ENDDO
  804. enddo
  805. enddo
  806. !$acc end region
  807. END SUBROUTINE wsm32D !}
  808. #else
  809. !===================================================================
  810. !
  811. SUBROUTINE wsm32D(t, q, qci, qrs,w, den, p, delz &
  812. ,delt,g, cpd, cpv, rd, rv, t0c &
  813. ,ep1, ep2, qmin &
  814. ,XLS, XLV0, XLF0, den0, denr &
  815. ,cliq,cice,psat &
  816. ,lat &
  817. ,rain, rainncv &
  818. ,snow,snowncv &
  819. ,sr &
  820. ,ids,ide, jds,jde, kds,kde &
  821. ,ims,ime, jms,jme, kms,kme &
  822. ,its,ite, jts,jte, kts,kte &
  823. )
  824. !-------------------------------------------------------------------
  825. IMPLICIT NONE
  826. !-------------------------------------------------------------------
  827. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
  828. ims,ime, jms,jme, kms,kme, &
  829. its,ite, jts,jte, kts,kte, &
  830. lat
  831. REAL, DIMENSION( its:ite , kts:kte ), &
  832. INTENT(INOUT) :: &
  833. t
  834. REAL, DIMENSION( ims:ime , kms:kme ), &
  835. INTENT(INOUT) :: &
  836. q, &
  837. qci, &
  838. qrs
  839. REAL, DIMENSION( ims:ime , kms:kme ), &
  840. INTENT(IN ) :: w, &
  841. den, &
  842. p, &
  843. delz
  844. REAL, INTENT(IN ) :: delt, &
  845. g, &
  846. cpd, &
  847. cpv, &
  848. t0c, &
  849. den0, &
  850. rd, &
  851. rv, &
  852. ep1, &
  853. ep2, &
  854. qmin, &
  855. XLS, &
  856. XLV0, &
  857. XLF0, &
  858. cliq, &
  859. cice, &
  860. psat, &
  861. denr
  862. REAL, DIMENSION( ims:ime ), &
  863. INTENT(INOUT) :: rain, &
  864. rainncv
  865. REAL, DIMENSION( ims:ime ), OPTIONAL, &
  866. INTENT(INOUT) :: snow, &
  867. snowncv, &
  868. sr
  869. ! LOCAL VAR
  870. REAL, DIMENSION( its:ite , kts:kte ) :: &
  871. rh, &
  872. qs, &
  873. denfac, &
  874. rslope, &
  875. rslope2, &
  876. rslope3, &
  877. rslopeb
  878. REAL, DIMENSION( its:ite , kts:kte ) :: &
  879. pgen, &
  880. pisd, &
  881. paut, &
  882. pacr, &
  883. pres, &
  884. pcon
  885. REAL, DIMENSION( its:ite , kts:kte ) :: &
  886. fall, &
  887. falk, &
  888. xl, &
  889. cpm, &
  890. work1, &
  891. work2, &
  892. xni, &
  893. qs0, &
  894. n0sfac
  895. REAL, DIMENSION( its:ite , kts:kte ) :: &
  896. falkc, &
  897. work1c, &
  898. work2c, &
  899. fallc
  900. INTEGER, DIMENSION( its:ite ) :: kwork1,&
  901. kwork2
  902. INTEGER, DIMENSION( its:ite ) :: mstep, &
  903. numdt
  904. LOGICAL, DIMENSION( its:ite ) :: flgcld
  905. REAL :: pi, &
  906. cpmcal, xlcal, lamdar, lamdas, diffus, &
  907. viscos, xka, venfac, conden, diffac, &
  908. x, y, z, a, b, c, d, e, &
  909. fallsum, fallsum_qsi, vt2i,vt2s,acrfac, &
  910. qdt, pvt, qik, delq, facq, qrsci, frzmlt, &
  911. snomlt, hold, holdrs, facqci, supcol, coeres, &
  912. supsat, dtcld, xmi, qciik, delqci, eacrs, satdt, &
  913. qimax, diameter, xni0, roqi0, supice,holdc, holdci
  914. INTEGER :: i, j, k, mstepmax, &
  915. iprt, latd, lond, loop, loops, ifsat, kk, n
  916. ! Temporaries used for inlining fpvs function
  917. REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
  918. ! variables for optimization
  919. REAL, DIMENSION( its:ite ) :: tvec1
  920. !
  921. !=================================================================
  922. ! compute internal functions
  923. !
  924. cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
  925. xlcal(x) = xlv0-xlv1*(x-t0c)
  926. !----------------------------------------------------------------
  927. ! size distributions: (x=mixing ratio, y=air density):
  928. ! valid for mixing ratio > 1.e-9 kg/kg.
  929. !
  930. ! Optimizatin : A**B => exp(log(A)*(B))
  931. lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
  932. lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
  933. !
  934. !----------------------------------------------------------------
  935. ! diffus: diffusion coefficient of the water vapor
  936. ! viscos: kinematic viscosity(m2s-1)
  937. !
  938. diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
  939. viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
  940. xka(x,y) = 1.414e3*viscos(x,y)*y
  941. diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
  942. ! venfac(a,b,c) = (viscos(b,c)/diffus(b,a))**(.3333333) &
  943. ! /viscos(b,c)**(.5)*(den0/c)**0.25
  944. venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
  945. /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
  946. conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
  947. !
  948. pi = 4. * atan(1.)
  949. !
  950. !----------------------------------------------------------------
  951. ! paddint 0 for negative values generated by dynamics
  952. !
  953. do k = kts, kte
  954. do i = its, ite
  955. qci(i,k) = max(qci(i,k),0.0)
  956. qrs(i,k) = max(qrs(i,k),0.0)
  957. enddo
  958. enddo
  959. !
  960. !----------------------------------------------------------------
  961. ! latent heat for phase changes and heat capacity. neglect the
  962. ! changes during microphysical process calculation
  963. ! emanuel(1994)
  964. !
  965. do k = kts, kte
  966. do i = its, ite
  967. cpm(i,k) = cpmcal(q(i,k))
  968. xl(i,k) = xlcal(t(i,k))
  969. enddo
  970. enddo
  971. !
  972. !----------------------------------------------------------------
  973. ! compute the minor time steps.
  974. !
  975. loops = max(nint(delt/dtcldcr),1)
  976. dtcld = delt/loops
  977. if(delt.le.dtcldcr) dtcld = delt
  978. !
  979. do loop = 1,loops
  980. !
  981. !----------------------------------------------------------------
  982. ! initialize the large scale variables
  983. !
  984. do i = its, ite
  985. mstep(i) = 1
  986. flgcld(i) = .true.
  987. enddo
  988. !
  989. ! do k = kts, kte
  990. ! do i = its, ite
  991. ! denfac(i,k) = sqrt(den0/den(i,k))
  992. ! enddo
  993. ! enddo
  994. do k = kts, kte
  995. CALL VREC( tvec1(its), den(its,k), ite-its+1)
  996. do i = its, ite
  997. tvec1(i) = tvec1(i)*den0
  998. enddo
  999. CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
  1000. enddo
  1001. !
  1002. ! Inline expansion for fpvs
  1003. ! qs(i,k) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
  1004. ! qs0(i,k) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
  1005. cvap = cpv
  1006. hvap=xlv0
  1007. hsub=xls
  1008. ttp=t0c+0.01
  1009. dldt=cvap-cliq
  1010. xa=-dldt/rv
  1011. xb=xa+hvap/(rv*ttp)
  1012. dldti=cvap-cice
  1013. xai=-dldti/rv
  1014. xbi=xai+hsub/(rv*ttp)
  1015. do k = kts, kte
  1016. do i = its, ite
  1017. ! tr=ttp/t(i,k)
  1018. ! if(t(i,k).lt.ttp) then
  1019. ! qs(i,k) =psat*(tr**xai)*exp(xbi*(1.-tr))
  1020. ! else
  1021. ! qs(i,k) =psat*(tr**xa)*exp(xb*(1.-tr))
  1022. ! endif
  1023. ! qs0(i,k) =psat*(tr**xa)*exp(xb*(1.-tr))
  1024. tr=ttp/t(i,k)
  1025. if(t(i,k).lt.ttp) then
  1026. qs(i,k) =psat*(exp(log(tr)*(xai)))*exp(xbi*(1.-tr))
  1027. else
  1028. qs(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
  1029. endif
  1030. qs0(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
  1031. qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k)
  1032. qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
  1033. qs(i,k) = max(qs(i,k),qmin)
  1034. rh(i,k) = max(q(i,k) / qs(i,k),qmin)
  1035. enddo
  1036. enddo
  1037. !
  1038. !----------------------------------------------------------------
  1039. ! initialize the variables for microphysical physics
  1040. !
  1041. !
  1042. do k = kts, kte
  1043. do i = its, ite
  1044. pres(i,k) = 0.
  1045. paut(i,k) = 0.
  1046. pacr(i,k) = 0.
  1047. pgen(i,k) = 0.
  1048. pisd(i,k) = 0.
  1049. pcon(i,k) = 0.
  1050. fall(i,k) = 0.
  1051. falk(i,k) = 0.
  1052. fallc(i,k) = 0.
  1053. falkc(i,k) = 0.
  1054. xni(i,k) = 1.e3
  1055. enddo
  1056. enddo
  1057. !
  1058. !----------------------------------------------------------------
  1059. ! compute the fallout term:
  1060. ! first, vertical terminal velosity for minor loops
  1061. !---------------------------------------------------------------
  1062. ! n0s: Intercept parameter for snow [m-4] [HDC 6]
  1063. !---------------------------------------------------------------
  1064. do k = kts, kte
  1065. do i = its, ite
  1066. supcol = t0c-t(i,k)
  1067. n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
  1068. if(t(i,k).ge.t0c) then
  1069. if(qrs(i,k).le.qcrmin)then
  1070. rslope(i,k) = rslopermax
  1071. rslopeb(i,k) = rsloperbmax
  1072. rslope2(i,k) = rsloper2max
  1073. rslope3(i,k) = rsloper3max
  1074. else
  1075. rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
  1076. ! rslopeb(i,k) = rslope(i,k)**bvtr
  1077. rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
  1078. rslope2(i,k) = rslope(i,k)*rslope(i,k)
  1079. rslope3(i,k) = rslope2(i,k)*rslope(i,k)
  1080. endif
  1081. else
  1082. if(qrs(i,k).le.qcrmin)then
  1083. rslope(i,k) = rslopesmax
  1084. rslopeb(i,k) = rslopesbmax
  1085. rslope2(i,k) = rslopes2max
  1086. rslope3(i,k) = rslopes3max
  1087. else
  1088. rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
  1089. ! rslopeb(i,k) = rslope(i,k)**bvts
  1090. rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
  1091. rslope2(i,k) = rslope(i,k)*rslope(i,k)
  1092. rslope3(i,k) = rslope2(i,k)*rslope(i,k)
  1093. endif
  1094. endif
  1095. !-------------------------------------------------------------
  1096. ! Ni: ice crystal number concentraiton [HDC 5c]
  1097. !-------------------------------------------------------------
  1098. ! xni(i,k) = min(max(5.38e7 &
  1099. !

Large files files are truncated, but you can click here to view the full file