| blob | cb5147f0bb57abda1eebef083a81eaf78106eb7c |
1 module module_mp_jensen_ishmael
3 !------------------------------------------------------------------------------------------------------!
4 ! This is the initial release of the Ice-Spheroids Habit Model with Aspect-ratio Evolution (ISHMAEL). !
5 ! This code calculates the moisture tendencies of water vapor, cloud water, rainwater, and three ice !
6 ! species which are planar-nucleated particles, columnar-nucleated particles, and aggregates. Planar- !
7 ! and columnar-nucleated ice particles are parameterized as spheroids (two growth dimensions are !
8 ! tracked), and ice particle shape and density evolution are predicted for vapor growth, sublimation, !
9 ! riming, and melting. The model can produce various vapor-grown habits (which currently depend on the !
10 ! environmental temperature and supersaturation), and the model predicts the change in ice particle !
11 ! shape and density during riming. Ice particle fall speeds are diagnosed based on ice particle !
12 ! properties. Aggregates are currently parameterized as spheroids with a fixed shape and density. Ice !
13 ! is initiated as planar or columnar based on temperature. After nucleation, the shapes of particles !
14 ! are not constrained, i.e., planar-nucleated particles can attain columnar shapes if the growth !
15 ! conditions allow it. !
16 ! !
17 ! Predicting two ice particle dimensions requires two prognostic variables. Thus mass, number, and two !
18 ! size-related mixing ratio are predicted for each ice species. That makes 12 prognostic variables !
19 ! for ice. Just throw more processors at it. Rain mass and number and cloud mass are also predicted !
20 ! See Jensen et al. 2017 (JAS), Jensen et al. 2018 (MWR), and Jensen et al. 2018 (JAS) for details. !
21 ! !
22 ! Author: Anders A. Jensen, NCAR-RAL, ajensen@ucar.edu, 303-497-8484 !
23 ! Last modified: 07 June 2019 (ISHMAEL version 0.1) !
24 !------------------------------------------------------------------------------------------------------!
26 use module_wrf_error
28 implicit none !.. You are welcome
30 !.. Constants
31 real, private, parameter :: &
32 PI = 3.14159265, & !.. Pi = 3.14159265... until forever
33 G_HOME = 9.8, & !.. Gravity (m s^-1)
34 RD = 287.15, & !.. Dry air gas constant (J kg^-1 K^-1)
35 RV = 461.5, & !.. Moist air gas constant (J kg^-1 K^-1)
36 CP = 1004.0, & !.. Heat capacity Air (J kg^-1 K^-1)
37 CPW = 4218.0, & !.. Heat capacity Water (J kg^-1 K^-1)
38 RCP = 0.285856574, & !.. RD / CP
39 RHOW = 1000.0, & !.. Water density (kg m^-3)
40 T0 = 273.15, & !.. STP temperature (K)
41 P0 = 101325.0, & !.. STP pressure (Pa)
42 RHOI = 920.0, & !.. Bulk ice density (kg m^-3)
43 R0 = 1.27494 !.. Air density at 1000 mb, T=T0 (kg m^-3)
45 !.. The lookup table arrays are itab for ice-cloud water collection
46 !.. and itabr for ice-rain collection. The lookup tables are read in
47 !.. the subroutine jensen_ishmael_init which occurs once.
48 real, private, allocatable, dimension(:,:,:,:,:) :: itab, itabr
50 !.. Inherent growth ratio data is stored in the igrdata array
51 !.. The gamma_tab is a tabulation of the gamma function needed
52 !.. because ice particle shape evolves.
53 real, private, allocatable, dimension(:) :: igrdata, gamma_tab
55 !.. Aggregation lookup tables arrays
56 real, private, allocatable, dimension(:,:,:) :: coltab, coltabn
58 !.. Added for double precision
59 integer, parameter :: LUT_KIND_R4 = selected_real_kind(6) !.. 4 byte real
60 real(kind=LUT_KIND_R4), private, allocatable, dimension(:,:,:,:,:) :: itab_o, itabr_o
61 real(kind=LUT_KIND_R4), private, allocatable, dimension(:) :: gamma_tab_o
63 contains
64 !--------------------------------------------------------------------------------------------------------------!
65 !--------------------------------------------------------------------------------------------------------------!
66 !.. This subroutine currently reads in the ice-cloud collection,
67 !.. ice-rain collecion, gamma calculations, and inherent growth ratio
68 !.. data lookup tables and calculates the ice-ice aggregation lookup table.
70 subroutine jensen_ishmael_init
71 implicit none
73 integer :: i1, i2, i3, i4, icat
74 integer, parameter :: ncat=7,npair=35,ndn=60
75 integer :: ipair(ncat,ncat)
76 real, dimension(ncat,9) :: dstprms
77 real, dimension(ncat) :: gnu,cfmas,pwmas,cfvt,pwvt,tablo,tabhi
79 !.. Collection pairs
80 ipair=&
81 reshape((/0,0,3,7,11,15,19 ,1,0,5,9,13,17,21 &
82 ,2,4,22,24,0,0,0 ,6,8,23,28,0,0,0 &
83 ,10,12,25,29,35,0,0 ,14,16,26,30,32,0,0 &
84 ,18,20,27,31,33,34,0/), &
85 (/ncat,ncat/))
87 !.. Ice-ice collection lookup table properties (JYH)
88 dstprms=&
89 !-----------------------------------------------------------------------
90 !-----------------------------------------------------------------------
91 ! cloud rain planar columnar aggreg graup hail
92 ! 1 2 3 4 5 6 7 ---> ncat
93 !-----------------------------------------------------------------------
94 reshape((/.5, .5, .318, .318, .5, .5, .5, & !.. shape (capacitance shape)
95 524., 524., .333, .333, .496, 157., 471., & !.. cfmas (m=cfmas*D**pwmas)
96 3., 3., 2.4, 2.4, 2.4, 3., 3., & !.. pwmas
97 3173., 149., 4.836, 4.836, 3.084, 93.3, 161., & !.. cfvt (v = cfvt*D**pwvt)
98 2., .5, 0.25, .25, .2, .5, .5, & !.. pwvt
99 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, & !.. tablo (size limits for table)
100 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, & !.. tabhi
101 .1e-5, .1e-4, .1e-5, .1e-5, .1e-4, .1e-4, .1e-4, & !.. dnmin (characteristic diam
102 .1e-2, .1e-1, .1e-1, .1e-1, .1e-1, .1e-1, .1e-1/), & !.. dnmax limits)
103 (/ncat,9/))
105 do icat=1,ncat
106 cfmas(icat)=dstprms(icat,2)
107 pwmas(icat)=dstprms(icat,3)
108 cfvt (icat)=dstprms(icat,4)
109 pwvt (icat)=dstprms(icat,5)
110 tablo(icat)=dstprms(icat,6)
111 tabhi(icat)=dstprms(icat,7)
112 gnu(icat) = 4.0
113 enddo
115 !.. Read in ice-cloud water collection lookup table
116 if(.not.allocated(itab_o)) then
117 allocate(itab_o(51,51,51,11,2))
119 open(20, FILE='ishmael-qi-qc.bin', FORM='unformatted')
120 do i1 = 1,51
121 do i2 = 1,51
122 do i3 = 1,51
123 do i4 = 1,11
124 read(20) itab_o(i1,i2,i3,i4,1), &
125 itab_o(i1,i2,i3,i4,2)
126 enddo
127 enddo
128 enddo
129 enddo
130 close(20)
131 print*, 'Jensen_ISHMAEL lookup table 1 of 3'
132 endif
134 !.. Read in ice-rain collection lookup table
135 if(.not.allocated(itabr_o)) then
136 allocate(itabr_o(51,51,51,11,6))
138 open(30, FILE='ishmael-qi-qr.bin', FORM='unformatted')
139 do i1 = 1,51
140 do i2 = 1,51
141 do i3 = 1,51
142 do i4 = 1,11
143 read(30) itabr_o(i1,i2,i3,i4,1), &
144 itabr_o(i1,i2,i3,i4,2), &
145 itabr_o(i1,i2,i3,i4,3), &
146 itabr_o(i1,i2,i3,i4,4), &
147 itabr_o(i1,i2,i3,i4,5), &
148 itabr_o(i1,i2,i3,i4,6)
149 enddo
150 enddo
151 enddo
152 enddo
153 close(30)
154 print*, 'Jensen_ISHMAEL lookup table 2 of 3'
155 endif
157 !.. Read in gamma function lookup table for speedy code
158 if(.not.allocated(gamma_tab_o)) then
159 allocate(gamma_tab_o(505001))
161 open(40, FILE='ishmael-gamma-tab.bin', FORM='unformatted')
162 do i1 = 1,505001
163 read(40) gamma_tab_o(i1)
164 enddo
165 close(40)
166 print*, 'Jensen_ISHMAEL lookup table 3 of 3'
167 endif
169 !.. Read in inhernent growth ratio data (-1C to -60C)
170 !.. This currently assumes columns below -20C (plates are used later in the code)
171 !.. From Chen and Lamb 1994
172 if(.not.allocated(igrdata)) then
173 allocate(igrdata(60))
175 igrdata = &
176 (/ 0.910547, 0.81807, 0.6874, 0.60127, 1.59767, 2.32423, 2.08818, &
177 1.61921, 1.15865, 0.863071, 0.617586, 0.453917, 0.351975, 0.28794, &
178 0.269298, 0.28794, 0.333623, 0.418883, 0.56992, 0.796458, 1.14325, &
179 1.64103, 1.90138, 1.82653, 1.61921, 1.47436, 1.32463, 1.25556, &
180 1.22239, 1.206, 1.11522, 1.10751, 1.10738, 1.11484, 1.12234, &
181 1.12221, 1.14529, 1.16884, 1.20104, 1.22573, 1.25094, 1.27666, &
182 1.31183, 1.3388, 1.35704, 1.37553, 1.38479, 1.39411, 1.40349, &
183 1.41294, 1.42245, 1.43202, 1.44166, 1.45137, 1.46114, 1.47097, &
184 1.48087, 1.50105, 1.50087, 1.51098 /)
185 endif
187 !.. make collection table (JYH)
188 !.. called only on initialization of the model
189 !.. coltab = mass mixing ratio collection
190 !.. coltabn = number mixing ratio collection
191 !.. ndn = table elements
192 !.. ncat = # of hydrometeor categories
193 if(.not.allocated(coltab)) then
194 allocate(coltab(ndn,ndn,npair))
195 endif
197 if(.not.allocated(coltabn)) then
198 allocate(coltabn(ndn,ndn,npair))
199 endif
201 !.. Call to build the aggregation lookup table (JYH)
202 call mkcoltb(ndn,ncat,coltab,coltabn,ipair,gnu,tablo,tabhi,cfmas,pwmas,cfvt,pwvt)
203 print*, 'Jensen_ISHMAEL aggregation lookup table built'
205 !.. Check that tables are allocated (Thanks to JYH and Shaofeng)
206 if(.not.allocated(itab)) allocate(itab(51,51,51,11,2))
207 if(.not.allocated(itabr)) allocate(itabr(51,51,51,11,6))
208 if(.not.allocated(gamma_tab)) allocate(gamma_tab(505001))
210 itab = itab_o; itabr = itabr_o; gamma_tab = gamma_tab_o
212 end subroutine jensen_ishmael_init
214 !--------------------------------------------------------------------------------------------------------------!
215 !--------------------------------------------------------------------------------------------------------------!
216 !.. The main subroutine.
218 subroutine MP_JENSEN_ISHMAEL(ITIMESTEP, DT_IN, P, DZ, &
219 TH, QV, QC, QR, NR, QI1, NI1, AI1, CI1, &
220 QI2, NI2, AI2, CI2, &
221 QI3, NI3, AI3, CI3, &
222 IDS,IDE, JDS,JDE, KDS,KDE, &
223 IMS,IME, JMS,JME, KMS,KME, &
224 ITS,ITE, JTS,JTE, KTS,KTE, &
225 RAINNC, RAINNCV, SNOWNC, SNOWNCV, &
226 diag_effc3d, diag_effi3d, diag_dbz3d, &
227 diag_vmi3d_1, diag_di3d_1, diag_rhopo3d_1, diag_phii3d_1, &
228 diag_vmi3d_2, diag_di3d_2, diag_rhopo3d_2, diag_phii3d_2, &
229 diag_vmi3d_3, diag_di3d_3, diag_rhopo3d_3, diag_phii3d_3, &
230 diag_itype_1,diag_itype_2,diag_itype_3 &
231 )
233 !.. ITIMESTEP - time step counter (integer)
234 !.. DT_IN - model time step (s)
235 !.. P - air pressure (Pa)
236 !.. DZ - difference in height over interface (m)
237 !.. TH - potential temperature (K)
238 !.. QV - water vapor mixing ratio (kg kg^-1)
239 !.. QC - cloud water mixing ratio (kg kg^-1)
240 !.. QR - rain water mixing ratio (kg kg^-1)
241 !.. NR - rain number concentration (kg^-1)
242 !.. QI1 - planar-nucleated mixing ratio (kg kg^-1)
243 !.. NI1 - planar-nucleated number mixing ratio (kg^-1)
244 !.. AI1 - planar-nucleated a^2c mixing ratio (m^3 kg^-1)
245 !.. CI1 - planar-nucleated c^2a mixing ratio (m^3 kg^-1)
246 !.. QI2 - columnar-nucleated mixing ratio (kg kg^-1)
247 !.. NI2 - columnar-nucleated number mixing ratio (kg^-1)
248 !.. AI2 - columnar-nucleated a^2c mixing ratio (m^3 kg^-1)
249 !.. CI2 - columnar-nucleated c^2a mixing ratio (m^3 kg^-1)
250 !.. QI3 - aggregate mixing ratio (kg kg^-1)
251 !.. NI3 - aggregate number mixing ratio (kg^-1)
252 !.. AI3 - aggregate a^2c mixing ratio (m^3 kg^-1)
253 !.. CI3 - aggregate c^2a mixing ratio (m^3 kg^-1)
254 !.. IDS,IDE - I domain start/end locations
255 !.. JDS,JDE - J domain start/end locations
256 !.. KDS,KDE - K domain start/end locations
257 !.. IMS,IME - I memory start/end locations
258 !.. JMS,JME - J memory start/end locations
259 !.. KMS,KME - K memory start/end locations
260 !.. ITS,ITE - I tile start/end locations
261 !.. JTS,JTE - J tile start/end locations
262 !.. KTS,KTE - K tile start/end locations
263 !.. RAINNC - total accumulated precipitation (mm)
264 !.. RAINNCV - total accumulated precipitation in a time step (mm)
265 !.. SNOWNC - total ice accumulated precipitation (mm)
266 !.. SNOWNCV - total ice accumulated precipitation in a time step (mm)
267 !.. DIAG_EFFC3D - Effective cloud water radius for coupling to radiation (m)
268 !.. DIAG_EFFI3D - Effective ice radius for coupling to radiation (m)
269 !.. DIAG_DBZ3D - 10-cm reflectivity (dBZ)
270 !.. DIAG_VMI3D_1 - planar-nucleated mass-weighted fall speeds (m s^-1)
271 !.. DIAG_VMI3D_2 - columnar-nucleated mass-weighted fall speeds (m s^-1)
272 !.. DIAG_VMI3D_3 - aggregate mass-weighted fall speeds (m s^-1)
273 !.. DIAG_DI3D_1 - planar-nucleated mass-weighted maximum diameter (m)
274 !.. DIAG_DI3D_2 - columnar-nucleated mass-weighted maximum diameter (m)
275 !.. DIAG_DI3D_3 - aggregate mass-weighted maximum diameter (m)
276 !.. DIAG_RHOPO3D_1 - planar-nucleated mass-weighted effective density (kg m^-3)
277 !.. DIAG_RHOPO3D_2 - columnar-nucleated mass-weighted effective density (kg m^-3)
278 !.. DIAG_RHOPO3D_3 - aggregate mass-weighted effective density (kg m^-3)
279 !.. DIAG_PHII3D_1 - planar-nucleated number-weighted aspect ratio
280 !.. DIAG_PHII3D_2 - columnar-nucleated number-weighted aspect ratio
281 !.. DIAG_PHII3D_3 - aggregate number-weighted aspect ratio
282 !.. DIAG_ITYPE_1 - currently not used
283 !.. DIAG_ITYPE_2 - currently not used
284 !.. DIAG_ITYPE_3 - currently not used
286 implicit none
288 !.. Grid variables
289 integer, intent(in) :: &
290 ids, ide, jds, jde, kds, kde , &
291 ims, ime, jms, jme, kms, kme , &
292 its, ite, jts, jte, kts, kte
294 !.. Time and timestep
295 integer, intent(in) :: ITIMESTEP
296 real, intent(in) :: dt_in
298 !.. 3D variables
299 real, dimension(ims:ime, kms:kme, jms:jme), intent(in) :: p, dz
301 real, dimension(ims:ime, kms:kme, jms:jme), intent(inout) :: &
302 th, qv, qc, qr, nr, qi1, ni1, ai1, ci1, qi2, ni2, ai2, ci2, &
303 qi3, ni3, ai3, ci3, &
304 diag_effc3d, diag_effi3d, diag_dbz3d, &
305 diag_vmi3d_1, diag_di3d_1, diag_rhopo3d_1, diag_phii3d_1, &
306 diag_vmi3d_2, diag_di3d_2, diag_rhopo3d_2, diag_phii3d_2, &
307 diag_vmi3d_3, diag_di3d_3, diag_rhopo3d_3, diag_phii3d_3, &
308 diag_itype_1, diag_itype_2, diag_itype_3
310 !.. 2D variables
311 real, dimension(ims:ime, jms:jme), intent(inout) :: &
312 rainnc, rainncv, snownc, snowncv
314 integer :: I,K,J
316 !.. Define the ice species
317 !.. Do not change these without changing/knowing how nucleation works
318 integer, parameter :: CAT = 3 !.. Number of ice categories (species)
319 integer, parameter :: ICE1 = 1 !.. Ice initiated as planar
320 integer, parameter :: ICE2 = 2 !.. Ice initiated as columnar
321 integer, parameter :: ICE3 = 3 !.. Aggregates
323 !.. Use tile locations for our variables
324 !.. Make everthing 1D for our loop
325 !.. These variables go in and out of newmicro
326 real, dimension(kts:kte) :: &
327 th1d, qv1d, qc1d, qr1d, nr1d, p1d, dz1d, &
328 effc1d, effi1d, dbz1d
330 real, dimension(cat, kts:kte) :: &
331 qi1d, ni1d, ai1d, ci1d, &
332 vmi1d, di1d, rhopo1d, phii1d, icetype1d
334 !.. Precipitation
335 real :: qrpre1d
336 real, dimension(cat) :: qipre1d
338 !--------------------------------------------------------------------------------------------------------------!
339 !--------------------------------------------------------------------------------------------------------------!
340 !.. Accumulated precipitation in one timestep (zero out before microphysics)
341 qrpre1d =0.
342 qipre1d(ICE1)=0.
343 qipre1d(ICE2)=0.
344 qipre1d(ICE3)=0.
346 !.. Loop over the grid and send each vertical column into the microphysics routine
347 do j=jts,jte !.. j loop (north-south)
348 do i=its,ite !.. i loop (horizontal)
350 !.. Get the 1D variables from WRF and zero out process rates
351 do k=kts,kte !.. k loop (vertical)
353 !.. Thermo,rain,cloud,vapor,etc
354 TH1D(k) = TH(i,k,j)
355 QV1D(k) = QV(i,k,j)
356 QC1D(k) = QC(i,k,j)
357 QR1D(k) = QR(i,k,j)
358 NR1D(k) = NR(i,k,j)
359 P1D(k) = P(i,k,j)
360 DZ1D(k) = DZ(i,k,j)
361 EFFC1D(k) = 0.
362 EFFI1D(k) = 0.
363 DBZ1D(k) = -35.
364 !.. Planar-nucleated
365 QI1D(ICE1,k) = QI1(i,k,j)
366 NI1D(ICE1,k) = NI1(i,k,j)
367 AI1D(ICE1,k) = AI1(i,k,j)
368 CI1D(ICE1,k) = CI1(i,k,j)
369 VMI1D(ICE1,k) = 0.
370 DI1D(ICE1,k) = 0.
371 RHOPO1D(ICE1,k) = 0.
372 PHII1D(ICE1,k) = 0.
373 ICETYPE1D(ICE1,k) = 0.
374 !.. Columnar-nucleated
375 QI1D(ICE2,k) = QI2(i,k,j)
376 NI1D(ICE2,k) = NI2(i,k,j)
377 AI1D(ICE2,k) = AI2(i,k,j)
378 CI1D(ICE2,k) = CI2(i,k,j)
379 VMI1D(ICE2,k) = 0.
380 DI1D(ICE2,k) = 0.
381 RHOPO1D(ICE2,k) = 0.
382 PHII1D(ICE2,k) = 0.
383 ICETYPE1D(ICE2,k) = 0.
384 !.. Aggregates
385 QI1D(ICE3,k) = QI3(i,k,j)
386 NI1D(ICE3,k) = NI3(i,k,j)
387 AI1D(ICE3,k) = AI3(i,k,j)
388 CI1D(ICE3,k) = CI3(i,k,j)
389 VMI1D(ICE3,k) = 0.
390 DI1D(ICE3,k) = 0.
391 RHOPO1D(ICE3,k) = 0.
392 PHII1D(ICE3,k) = 0.
393 ICETYPE1D(ICE3,k) = 0.
394 enddo !.. (vertical k-loop)
396 !.. Call me Ishmael.
397 call me_ishmael(itimestep, kts, kte, i, j, dt_in, p1d, dz1d, &
398 th1d, qv1d, qc1d, qr1d, nr1d, qi1d, ni1d, ai1d, ci1d, &
399 CAT, ICE1, ICE2, ICE3, effc1d, effi1d, dbz1d, icetype1d, &
400 vmi1d, di1d, rhopo1d, phii1d, qrpre1d, qipre1d)
402 !.. Send the 1D columns from newmicro back to WRF
403 do k=kts,kte !.. k loop (vertical)
404 TH(i,k,j) = TH1D(k)
405 QV(i,k,j) = QV1D(k)
406 QC(i,k,j) = QC1D(k)
407 QR(i,k,j) = QR1D(k)
408 NR(i,k,j) = NR1D(k)
409 DIAG_EFFC3D(i,k,j) = max((min(EFFC1D(k),50.e-6)),2.51e-6)
410 DIAG_EFFI3D(i,k,j) = max((min(EFFI1D(k),999.e-6)),5.e-6)
411 DIAG_DBZ3D(i,k,j) = max(dbz1d(k),-35.)
412 !.. Planar-nucleated
413 QI1(i,k,j) = QI1D(ICE1,k)
414 NI1(i,k,j) = NI1D(ICE1,k)
415 AI1(i,k,j) = AI1D(ICE1,k)
416 CI1(i,k,j) = CI1D(ICE1,k)
417 DIAG_VMI3D_1(i,k,j) = vmi1d(ICE1,k)
418 DIAG_DI3D_1(i,k,j) = di1d(ICE1,k)
419 DIAG_RHOPO3D_1(i,k,j)= rhopo1d(ICE1,k)
420 DIAG_PHII3D_1(i,k,j) = phii1d(ICE1,k)
421 DIAG_ITYPE_1(i,k,j) = icetype1d(ICE1,k)
422 !.. Columnar-nucleated
423 QI2(i,k,j) = QI1D(ICE2,k)
424 NI2(i,k,j) = NI1D(ICE2,k)
425 AI2(i,k,j) = AI1D(ICE2,k)
426 CI2(i,k,j) = CI1D(ICE2,k)
427 DIAG_VMI3D_2(i,k,j) = vmi1d(ICE2,k)
428 DIAG_DI3D_2(i,k,j) = di1d(ICE2,k)
429 DIAG_RHOPO3D_2(i,k,j)= rhopo1d(ICE2,k)
430 DIAG_PHII3D_2(i,k,j) = phii1d(ICE2,k)
431 DIAG_ITYPE_2(i,k,j) = icetype1d(ICE2,k)
432 !.. Aggregates
433 QI3(i,k,j) = QI1D(ICE3,k)
434 NI3(i,k,j) = NI1D(ICE3,k)
435 AI3(i,k,j) = AI1D(ICE3,k)
436 CI3(i,k,j) = CI1D(ICE3,k)
437 DIAG_VMI3D_3(i,k,j) = vmi1d(ICE3,k)
438 DIAG_DI3D_3(i,k,j) = di1d(ICE3,k)
439 DIAG_RHOPO3D_3(i,k,j)= rhopo1d(ICE3,k)
440 DIAG_PHII3D_3(i,k,j) = phii1d(ICE3,k)
441 DIAG_ITYPE_3(i,k,j) = icetype1d(ICE3,k)
442 enddo
444 !.. Accumulated precipitation and precipitation rates
445 !.. Frozen precipitation
446 SNOWNC(i,j) = SNOWNC(i,j) + qipre1d(ICE1) + qipre1d(ICE2) + qipre1d(ICE3)
447 SNOWNCV(i,j) = qipre1d(ICE1) + qipre1d(ICE2) + qipre1d(ICE3)
449 !.. Total precipitation
450 RAINNC(i,j) = RAINNC(i,j) + qrpre1d + qipre1d(ICE1) + qipre1d(ICE2) + qipre1d(ICE3)
451 RAINNCV(i,j) = qrpre1d + qipre1d(ICE1) + qipre1d(ICE2) + qipre1d(ICE3)
453 enddo ! i loop (horizontal)
454 enddo ! j loop (north-south)
456 end subroutine MP_JENSEN_ISHMAEL
458 !--------------------------------------------------------------------------------------------------------------!
459 !--------------------------------------------------------------------------------------------------------------!
460 !.. The microphysics
461 subroutine me_ishmael(it, kts, kte, i, j, dt, pres_e, dzmic, &
462 theta, qv, qc, qr, nr, qi, ni, ai, ci, CAT, ICE1, ICE2, ICE3, &
463 effc1d, effi1d, dbz1d, icetype1d, vmi1d, di1d, rhopo1d, phii1d, &
464 qrpre1d, qipre1d)
466 implicit none
468 !.. User logicals (all should be TRUE except for testing)
469 logical, parameter :: RIMING = .TRUE. ! Riming
470 logical, parameter :: RAIN_ICE = .TRUE. ! Rain-ice collection (for testing)
471 logical, parameter :: FREEZE_QC = .TRUE. ! Freeze cloud drops (Homogeneous)
472 logical, parameter :: SPLINTERS = .TRUE. ! Rime Splinters
473 logical, parameter :: WET_GROW = .TRUE. ! Wet growth
475 !.. Input/Output variables from WRF
476 integer, intent(in) :: it, i, j !.. Integer timestep, model i and j
477 integer, intent(in) :: kts, kte !.. Tile k-start and k-end
478 integer, intent(in) :: CAT, ICE1, ICE2, ICE3 !.. Number of categories, planar and columnar index
479 real, intent(in) :: dt !.. Timestep
481 real, dimension(kts:kte), intent(in) :: pres_e, dzmic !.. Pressure and dz
482 real, dimension(kts:kte), intent(inout) :: theta !.. Theta (K)
483 real, dimension(kts:kte), intent(inout) :: qv !.. Vapor mixing ratio (kg kg^-1)
484 real, dimension(kts:kte), intent(inout) :: qc !.. Cloud water mixing ratio (kg kg^-1)
485 real, dimension(kts:kte), intent(inout) :: qr !.. Rain mixing ratio (kg kg^-1)
486 real, dimension(kts:kte), intent(inout) :: nr !.. Rain number mixing ratio (# kg^-1)
488 real, dimension(cat, kts:kte), intent(inout) :: qi !.. Ice mass mixing ratio (kg kg^-1)
489 real, dimension(cat, kts:kte), intent(inout) :: ni !.. Ice number mixing ratio (# kg^-1)
490 real, dimension(cat, kts:kte), intent(inout) :: ai !.. Ice a^2c mixing ratio (m^3 kg^-1)
491 real, dimension(cat, kts:kte), intent(inout) :: ci !.. Ice c^2a mixing ratio (m^3 kg^-1)
493 !.. Ice diagnostics
494 real, dimension(kts:kte), intent(inout) :: effc1d, effi1d, dbz1d
495 real, dimension(cat, kts:kte), intent(inout) :: rhopo1d
496 real, dimension(cat, kts:kte), intent(inout) :: phii1d
497 real, dimension(cat, kts:kte), intent(inout) :: vmi1d
498 real, dimension(cat, kts:kte), intent(inout) :: di1d
499 real, dimension(cat, kts:kte), intent(inout) :: icetype1d
500 real :: dbzr,dbzsum
502 real :: qrpre1d !.. Rain precipitation accumulation (mm)
503 real, dimension(cat) :: qipre1d !.. Ice liquid-equivalent precipitation accumulation (mm)
505 !.. The Parameters
506 real, parameter :: NU = 4.0 !.. Ice distribution parameter (keep at 4)
507 !.. If changed, build new lookup tables
509 real, parameter :: QSMALL = 1.e-12 !.. Smallest ice mass (kg kg^-1)
510 real, parameter :: QNSMALL= 1.25e-7 !.. Smallest ice number (# kg^-1)
511 real, parameter :: QASMALL= 1.e-24 !.. Smallest ice volume (m^-3 kg^-1)
512 real, parameter :: ao = 0.1e-6 !.. alphstr=ao^(1-deltastr)
513 real, parameter :: AR = 149.1 !.. 'a' parameter in rain fallspeed-size relation (m^(1-br)/s)
514 real, parameter :: BR = 0.5 !.. 'b' parameter in rain fallspeed-size relation
515 real, parameter :: F1R = 0.78 !.. Ventilation coefficient for rain
516 real, parameter :: F2R = 0.308 !.. Ventilation coefficient for rain
517 real, parameter :: LAMMAXR = 1./20.E-6 !.. Max allowed value for rain slope parameter (m^-1)
518 real, parameter :: LAMMINR = 1./2800.E-6 !.. Min allowed value for rain slope parameter (m^-1)
520 !.. Local/Dummy variables
521 logical :: vgflag !.. Vapor growth flag to determine if vapor growth
522 !.. is producing subsaturated conditions
524 logical, dimension(cat) :: has_ice !.. true if ice exists for an ice species
525 logical, dimension(cat) :: dry_growth !.. Wet growth/dry growth riming flag
526 integer :: k, cc, ccvar !.. k is the k-loop index, cc is the ice species index
527 integer :: gi !.. Index for the Gamma(NU) lookup
528 integer :: gi2, gi3 !.. Index for the Gamma(NU) lookup
529 integer :: numice !.. Number of ice species currently with ice
530 integer :: current_index !.. Used to determine number of ice species wi th mass > QSMALL
531 real, allocatable, dimension(:) :: icearray !.. Array with size of number
532 !.. of ice species with mass > QSMALL
533 real :: wetg !.. wet growth check
534 real :: gamma_arg !.. Gamma argument: Gamma(gamma_arg)
535 real, dimension(cat) :: ani, cni, rni !.. Characteristic a-,c-,r-axis sizes
536 real, dimension(cat) :: aniold, cniold, dsold
537 real :: alphstr !.. alpha_* = ao^(1-delta_*)
538 real :: alphv !.. alphv = 4/3*PI*alpha_*
539 real :: betam !.. betam = 2 + delta_*
540 real :: anf, cnf, rnf, iwcf !.. a-,c-,r-axis and iwc after vapor growth
541 real :: anfr, cnfr, rnfr, iwcfr !.. a-,c-,r-axis and iwc after riming
542 real :: phibr, phifr !.. aspect ratio before and after riming
543 real :: capgam !.. distribution-weighted capacitance (m)
544 real :: i_dt !.. inverse timestep
545 real :: i_cp !.. inverse heat capacity of air
546 real :: sink, source, ratio !.. Terms used to prevent over-depletion
547 real :: rrr !.. rain characteristic radius
549 !.. Fall speed parameters
550 real :: phiivt, bl, al, aa, ba, qe, xn, bx, xm, f_c1, f_c2, bm, am, nre
552 !.. Local variables for riming
553 real :: vir, vif, vfr, rhobarrime, rimedr, rimedrr, rimec1, rhobarwg
555 !.. Local variables for ice nucleation, multiplication
556 real :: anuc, nibnuc(CAT), amass, cmass, fmult, nucfrac
557 real :: a_demott, b_demott, c_demott, d_demott
558 real :: curnum, ratel, ratekg, inrate
560 !.. Local variables for vapor growth
561 real :: vi, iwci, nim3dum, prdsum, qcrimesum, qcrimesumr
563 !.. Local variables for melting
564 real :: mfact, nimelt
566 !.. Local variables for evaporation
567 real :: epsr, dumqv, dumqvs, arn
569 !.. Some dummy variables
570 real :: dum, dum1, dumt, dumgam, dumgam2
571 real :: tmpsum, tmpsumr, tmpmelt
573 !.. Local variables for Radiation
574 real :: lrsig, sig, lwc, ncm3dum, r_n, phirad, radw1, radw2, radw3, totaliwc
576 !.. Local variables for Aggregation
577 real :: dn1, dn2, dn3, phiagg1, phiagg2, phiagg3
579 !.. Fall speeds (some not used) and effectice ice radius variables
580 real, dimension(kts:kte) :: vtrm, vtrn, vtrmc, vtrnc, nc
581 real, dimension(cat, kts:kte) :: vtrni1, effi, vtrmi1, vtrzi1
583 !.. Lookup table values
584 real :: rrri, rrni, rqci, rdsi, rrho, proc(cat,2), procr(cat,6)
585 integer :: irri, irni, iqci, idsi, irho, iti
586 real :: qi_qc_nrm(cat), qi_qc_nrd(cat)
587 real :: qi_qr_nrm(cat), qi_qr_nrd(cat), qi_qr_nrn(cat)
589 !.. Sedimentation variables
590 integer :: nn, nstep
591 real, dimension(kts:kte) :: fluxqr, fluxnr
592 real, dimension(cat, kts:kte) :: fluxqi, fluxni, fluxai, fluxci
593 real, dimension(cat) :: falltndqi, falltndni, falltndai, falltndci
594 real :: falltndqr, falltndnr, maxfall
595 logical :: sedi
596 logical :: domicro
598 !.. The rest of the variables
599 real, dimension(kts:kte) :: rhoair, i_rhoair !.. Air density (kg m^-3) and inverse air density
600 real :: qcrimefrac(cat) !.. Qc_Rime/(Qc_Rime + Qr_Rime)
601 real :: gammnu !.. Gamma(NU)
602 real :: i_gammnu !.. 1/Gamma(NU)
603 real :: fourthirdspi !.. 4/3 * PI
604 real :: mnuccd !.. Deposition nucleation rate (kg kg^-1 s^-1)
605 real :: nnuccd !.. Deposition nucleation rate (# kg^-1)
606 real :: pcc !.. Cloud mass growth rate (kg kg^-1 s^-1)
607 real :: qagg(cat) !.. Aggregation mass tendency (kg kg^-1)
608 real :: nagg(cat) !.. Aggregation number tendency (# kg^-1)
609 real :: prd(cat) !.. Vapor growth rate (kg kg^-1 s^-1)
610 real :: nrd(cat) !.. Vapor growth number rate sublimation (# kg^-1 s^-1)
611 real :: ard(cat) !.. A-axis vapor growth rate mixing ratio (m kg^-1 s^-1)
612 real :: crd(cat) !.. C-axis vapor growth rate mixing ratio (m kg^-1 s^-1)
613 real :: prdr(cat) !.. Riming growth rate (kg kg^-1 s^-1)
614 real :: ardr(cat) !.. A-axis riming growth rate mixing ratio (m kg^-1 s^-1)
615 real :: crdr(cat) !.. C-axis riming growth rate mixing ratio (m kg^-1 s^-1)
616 real :: qmlt(cat) !.. Mass melting rate (kg kg^-1 s^-1)
617 real :: nmlt(cat) !.. Number melting rate (# kg^-1 s^-1)
618 real :: amlt(cat) !.. A-axis melting rate (m kg^-1 s^-1)
619 real :: cmlt(cat) !.. C-axis melting rate (m kg^-1 s^-1)
620 real :: rimesum(cat) !.. Ice-cloud rime mass rate (kg m^-3 s^-1)
621 real :: rimesumr(cat) !.. Ice-rain rime mass rate (kg m^-3 s^-1)
622 real :: rimetotal(cat) !.. Total rime mass rate (kg m^-3 s^-1)
623 real :: dQRfzri(cat) !.. Change in rain mass from ice-rain (kg kg^-1 s^-1)
624 real :: dQIfzri(cat) !.. Change in ice mass from ice-rain (kg kg^-1 s^-1)
625 real :: dNfzri(cat) !.. Change in ice number from ice-rain (# kg^-1 s^-1)
626 real :: dQImltri(cat) !.. Change in ice mass from ice-rain melt (kg kg^-1 s^-1)
627 real :: dNmltri(cat) !.. Change in ice number from ice-rain melt (# kg^-1 s^-1)
628 real :: numrateri(cat) !.. Number rate from ice-rain (value from lookkup table)
629 real :: rainrateri(cat) !.. Rain mass rate from ice-rain (value from lookkup table)
630 real :: icerateri(cat) !.. Ice mass rate from ice-rain (value from lookkup table)
631 real :: totalnuc(cat) !.. Total mass nucleated to ice (kg kg^-1 s^-1)
632 real :: totalnucn(cat) !.. Total number nucleated to ice (# kg^-1 s^-1)
633 real :: nim !.. Number of homogeneously frozen cloud drops (# kg^-1 s^-1)
634 real :: mim !.. Mass of homogeneously frozen cloud drops (kg kg^-1 s^-1)
635 real :: nimr !.. Number of homogeneously frozen rain drops (# kg^-1 s^-1)
636 real :: mimr !.. Mass of homogeneously frozen rain drops (kg kg^-1 s^-1)
637 real :: deltastr(cat) !.. Delta_* (Average inherent growth ratio)
638 real :: rhobar(cat) !.. Ice effective density (kg m^-3)
639 real :: rhodepout(cat) !.. Ice effective density (kg m^-3)
640 real :: dsdepout(cat) !.. Ice effective density (kg m^-3)
641 real :: rhorimeout(cat) !.. Ice effective density (kg m^-3)
642 real :: gdentotal !.. Total rime density from ice-cloud and ice-rain (kg m^-3)
643 real :: gdenavg(cat) !.. Rime density from ice-cloud (kg m^-3)
644 real :: gdenavgr(cat) !.. Rime density from ice-rain (kg m^-3)
645 real :: fv(cat) !.. Ventilation coefficient for diffusion
646 real :: fh(cat) !.. Ventilation coefficient for heat
647 real :: temp, i_temp !.. Air temperature (K) and inverse
648 real :: qvs !.. Liquid saturation vapor mixing ratio (kg kg^-1)
649 real :: qvi !.. Ice saturation vapor mixing ratio (kg kg^-1)
650 real :: sup !.. Liquid supersaturation ratio
651 real :: sui !.. Ice supersaturation ratio
652 real :: svpi !.. Ice saturation vapor pressure
653 real :: svpl !.. Liquid saturation vapor pressure
654 real :: qs0 !.. Saturation vapor mixing ratio at 0C (kg kg^-1)
655 real :: xxls !.. Latent heat of sublimation (J kg^-1)
656 real :: xxlv !.. Latent heat of vaporization (J kg^-1)
657 real :: xxlf !.. Latent heat of fusion (J kg^-1)
658 real :: ab !.. Latent heat correction parameter for growth
659 real :: abi !.. Latent heat correction parameter for growth (ice)
660 real :: mu !.. Dynamic visocity of air (kg m^-1 s^-1)
661 real :: dv !.. Diffusivity of water vapor in air (m^2 s^-1)
662 real :: kt !.. Thermal conductivity (W m^-1 K^-1)
663 real :: nsch !.. Schmidt number
664 real :: npr !.. Prandlt number
665 real :: lamr !.. Rain slope parameter n(D)=n0rr*exp(-lamr*D)dD
666 real :: n0rr !.. Rain distribution parameter
667 real :: prc !.. Autoconversion cloud water to rain (kg kg^-1 s^-1)
668 real :: nprc !.. Autoconversion cloud number loss (# kg^-1 s^-1)
669 real :: nprc1 !.. Autoconversion rain number gain (# kg^-1 s^-1)
670 real :: pre !.. Rain evpaoration rate (kg kg^-1 s^-1)
671 real :: npre !.. Rain evaporation number (# kg^-1 s^-1)
672 real :: pra !.. Rain accretion of cloud (kg kg^-1 s^-1)
673 real :: npra !.. Rain accretion of cloud number (# kg^-1 s^-1)
674 real :: nragg !.. Rain self-collection (# kg^-1 s^-1)
675 real :: igr !.. Inherent growth ratio
676 real :: qmult(cat) !.. Ice multiplication rate (kg kg^-1 s^-1)
677 real :: nmult(cat) !.. Ice number multiplication rate (# kg^-1 s^-1)
678 real :: mbiggr !.. Biggs freezing (kg kg^-1 s^-1)
679 real :: nbiggr !.. Biggs number freezing (# kg^-1 s^-1)
681 !--------------------------------------------------------------------------------------------------------------!
682 !--------------------------------------------------------------------------------------------------------------!
683 !.. NU=4 is constant so gammnu=6
684 gammnu = gamma(NU)
685 i_gammnu = 1./gammnu !.. Inverse gamma
686 i_dt = 1./dt !.. Inverse timestep
687 i_cp = 1./CP !.. Inverse heat capacity
688 fourthirdspi = 4./3.*PI !.. Four divided by three times pi
690 !.. No need to loop through microphysics and sedimentation unless
691 !.. condensate exists
692 sedi = .false.
694 !--------------------------------------------------------------------------------------------------------------!
695 !--------------------------------------------------------------------------------------------------------------!
696 do k=kts,kte !.. Grid K loop
698 !.. Check the column for condensate or supersaturation wrt liquid.
699 !.. If it is there, then do the microphysics
700 domicro= .false.
701 temp = theta(k)/(100000./pres_e(k))**(RCP)
702 svpl = polysvp(temp,0)
703 qvs = 0.622*svpl/(pres_e(k)-svpl)
704 sup = qv(k)/qvs-1.
706 svpi = polysvp(temp,1)
707 qvi = 0.622*svpi/(pres_e(k)-svpi)
708 if(temp.gt.T0) qvi=qvs
709 sui = qv(k)/qvi-1.
711 if(qc(k).gt.QSMALL.or.qr(k).gt.QSMALL.or.qi(ICE1,k).gt.QSMALL.or.&
712 qi(ICE2,k).gt.QSMALL.or.qi(ICE3,k).gt.QSMALL.or.sup.ge.0.) then
713 domicro= .true.
714 endif
716 !.. If condensate or supersaturation wrt liquid, do microphysics
717 if(domicro) then
719 if(allocated(icearray)) deallocate(icearray)
721 numice = 0 !.. Number of species with ice
722 do cc = 1, cat !.. Loop over all ice species
723 has_ice(cc) = .false.
724 if(qi(cc,k).gt.QSMALL) then !.. Check for ice
726 !.. Get the number of species with ice
727 has_ice(cc) = .true.
728 numice = numice + 1
729 endif
730 enddo
732 !.. Determine the number of ice species to loop over
733 !.. For vapor growth/riming, etc
734 if(numice.gt.0.and..not.allocated(icearray)) then
735 allocate(icearray(numice))
736 current_index = 1
737 do cc = 1, cat
738 if(has_ice(cc)) then
739 icearray(current_index) = cc
740 current_index = current_index + 1
741 endif
742 enddo
743 endif
745 !.. Initialize variables for all ice species and process rates
746 do cc = 1, cat !.. Ice variables
747 dry_growth(cc)=.true.
748 qcrimefrac(cc)=0.
749 rimetotal(cc)=0.
750 nmult(cc)=0.
751 qmult(cc)=0.
752 ani(cc)=0.
753 cni(cc)=0.
754 rni(cc)=0.
755 qagg(cc)=0.
756 nagg(cc)=0.
757 prd(cc) =0.
758 nrd(cc) =0.
759 ard(cc) =0.
760 crd(cc) =0.
761 prdr(cc)=0.
762 ardr(cc)=0.
763 crdr(cc)=0.
764 qmlt(cc)=0.
765 nmlt(cc)=0.
766 amlt(cc)=0.
767 cmlt(cc)=0.
768 deltastr(cc)=1.
769 dsold(cc)=1.
770 dsdepout(cc)=1.
771 rhobar(cc)=RHOI
772 rhodepout(cc)=RHOI
773 rhorimeout(cc)=RHOI
774 if(cc.eq.3) rhobar(cc)=50.
775 fv(cc)=1.
776 fh(cc)=1.
777 rimesum(cc)=0.
778 rimesumr(cc)=0.
779 numrateri(cc)=0.
780 rainrateri(cc)=0.
781 icerateri(cc)=0.
782 dQRfzri(cc)=0.
783 dQIfzri(cc)=0.
784 dNfzri(cc)=0.
785 dQImltri(cc)=0.
786 dNmltri(cc)=0.
787 gdenavg(cc)=400.
788 gdenavgr(cc)=400.
789 totalnuc(cc)=0.
790 totalnucn(cc)=0.
791 nibnuc(cc)=0.
792 aniold(cc)=0.
793 cniold(cc)=0.
794 proc(cc,1)=0.
795 proc(cc,2)=0.
796 qi_qc_nrm(cc)=0.
797 qi_qc_nrd(cc)=0.
798 qi_qr_nrm(cc)=0.
799 qi_qr_nrn(cc)=0.
800 qi_qr_nrd(cc)=0.
801 procr(cc,1)=0.
802 procr(cc,2)=0.
803 procr(cc,3)=0.
804 procr(cc,4)=0.
805 procr(cc,5)=0.
806 procr(cc,6)=0.
807 fluxqi(cc,k)=0.
808 fluxni(cc,k)=0.
809 fluxai(cc,k)=0.
810 fluxci(cc,k)=0.
811 falltndqi(cc)=0.
812 falltndni(cc)=0.
813 falltndai(cc)=0.
814 falltndci(cc)=0.
815 di1d(cc,k)=0.
816 vmi1d(cc,k)=0.
817 vtrmi1(cc,k)=0.
818 vtrni1(cc,k)=0.
819 vtrzi1(cc,k)=0.
820 rhopo1d(cc,k)=0.
821 phii1d(cc,k)=0.
822 qipre1d(cc)=0.
823 effi(cc,k)=0.
824 icetype1d(cc,k)=0.
825 enddo
826 dbz1d(k)=-35.
827 qrpre1d=0.
828 fluxqr(k)=0.
829 fluxnr(k)=0.
830 falltndqr=0.; falltndnr=0.
831 vtrm(k)=0.; vtrn(k)=0.
832 effc1d(k) = 0.
833 effi1d(k) = 0.
834 nprc =0.; npre=0.; npra=0.; pra=0.; prc=0.
835 nprc1=0.; nragg=0.
836 nim =0.; mim =0.
837 nimr=0.; mimr=0.
838 mbiggr=0.; nbiggr=0.
839 mnuccd=0.; nnuccd=0.; pcc=0.; fmult=0.
840 gdentotal=400.
841 rimedr=0.; rimedrr=0.
842 dbzr=0.; dbzsum=0.
844 !--------------------------------------------------------------------------------------------------------------!
845 !--------------------------------------------------------------------------------------------------------------!
846 !.. Temperature and temperature-dependent variables
847 i_temp= 1./temp
848 rhoair(k)= pres_e(k)/(RD*temp)
849 i_rhoair(k) = 1./rhoair(k)
850 svpi = polysvp(temp,1)
851 qvi = 0.622*svpi/(pres_e(k)-svpi)
852 if(temp.gt.T0) qvi=qvs
853 qs0 = 0.622*polysvp(T0,0)/(pres_e(k)-polysvp(T0,0))
854 sui = qv(k)/qvi-1.
855 xxls = 3.15e6-2370.*temp+0.3337e6
856 xxlv = 3.1484E6-2370.*temp
857 xxlf = xxls - xxlv
858 ab = 1.+(xxlv*qvs/(RV*temp**2))*xxlv*i_cp
859 abi = 1.+(xxls*qvi/(RV*temp**2))*xxls*i_cp
860 mu = 1.496E-6*temp**1.5/(temp+120.)
861 arn = ar*(R0*i_rhoair(k))**0.5
862 dv = 8.794E-5*temp**1.81/pres_e(k)
863 kt = 2.3823e-2 +7.1177e-5*(temp-T0)
864 nsch = mu*i_rhoair(k)/dv
865 npr = mu*i_rhoair(k)/kt
867 !.. Inherent growth ratio data from Lamb and Scott 1972 and Chen and Lamb 1994
868 !.. Assumes planar particles at temperatures below -20C (see Bailey and Hallett)
869 !.. Function of temperature only (thanks to JHY and GHB for moving this fuction to where it belongs)
870 igr = get_igr(igrdata, temp)
871 if((temp-T0).lt.-20.) then
872 igr=0.7
873 endif
875 !.. Assume constant cloud drop number concentration (200 cm^-3)
876 nc(k) = 200.e6*i_rhoair(k)
878 !--------------------------------------------------------------------------------------------------------------!
879 !--------------------------------------------------------------------------------------------------------------!
880 !.. Ice Microphysics
882 if(numice.gt.0) then !.. If one or more ice species have ice
883 tmpsum = 0.
884 tmpsumr= 0.
885 do ccvar = 1, numice
886 cc = icearray(ccvar)
888 !.. Variable check
889 ni(cc,k)=max(ni(cc,k),QNSMALL)
890 ai(cc,k)=max(ai(cc,k),QASMALL)
891 ci(cc,k)=max(ci(cc,k),QASMALL)
892 ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
893 cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
895 !.. Set incoming volume limit after advection
896 !.. This is the smallest bulk volume for which we care about shape
897 !.. If smaller than this limit, assume spherical
898 if(ai(cc,k).lt.1.e-12.or.ci(cc,k).lt.1.e-12) then
899 ai(cc,k)=min(ai(cc,k),ci(cc,k))
900 ci(cc,k)=ai(cc,k)
901 ani(cc)=(ai(cc,k)/ni(cc,k))**0.3333333333
902 cni(cc)=(ci(cc,k)/ni(cc,k))**0.3333333333
903 endif
905 if(ani(cc).lt.2.e-6.or.cni(cc).lt.2.e-6) then
906 ani(cc)=2.e-6
907 cni(cc)=2.e-6
908 endif
910 ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k)
911 ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k)
913 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
915 !.. Check to keep ice properties in bounds
916 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
917 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
918 alphstr, alphv, betam)
920 dsold(cc) = deltastr(cc)
922 !--------------------------------------------------------------------------------------------------------------!
923 !--------------------------------------------------------------------------------------------------------------!
924 !.. Ice-cloud riming rates from lookup table
925 if(RIMING.and.qc(k).gt.1.e-7) then
927 !.. The lookup table is in log space. Do not change the following 4 lines without
928 !.. changing the offline-built lookup table
929 rrni = 13.498*log10(0.5e6*rni(cc))
930 rqci = 8.776*log10(1.0e7*(exp(qc(k))-1.))
931 rdsi = 50.*(deltastr(cc) - 0.5)
932 rrho = 7.888*log10(0.02*rhobar(cc))
933 irni = int(rrni)
934 iqci = int(rqci)
935 idsi = int(rdsi)
936 irho = int(rrho)
938 !.. Small limit (stay in the bounds of the lookup table)
939 rrni = max(rrni,1.)
940 rqci = max(rqci,1.)
941 rdsi = max(rdsi,1.)
942 rrho = max(rrho,1.)
944 irni = max(irni,1)
945 iqci = max(iqci,1)
946 idsi = max(idsi,1)
947 irho = max(irho,1)
949 !.. Large limit (stay in the bounds of the lookup table)
950 rrni = min(rrni,real(size(itab,1))-1.)
951 rqci = min(rqci,real(size(itab,2))-1.)
952 rdsi = min(rdsi,real(size(itab,3))-1.)
953 rrho = min(rrho,real(size(itab,4))-1.)
955 irni = min(irni,size(itab,1)-1)
956 iqci = min(iqci,size(itab,2)-1)
957 idsi = min(idsi,size(itab,3)-1)
958 irho = min(irho,size(itab,4)-1)
960 do iti = 1, 2
961 call access_lookup_table(itab,irni,iqci,idsi,irho,iti, &
962 rdsi,rrho,rqci,rrni,proc(cc,iti))
963 enddo
965 !.. proc(cc,1) is the normalized riming rate
966 !.. proc(cc,2) is related to the normalized rime density
967 !.. rimesum is the riming rate in kg m^-3 s^-1
968 rimesum(cc) = max((proc(cc,1)*ni(cc,k)*nc(k)*rhoair(k)**2),0.)
969 qi_qc_nrm(cc)= proc(cc,1)
970 qi_qc_nrd(cc)= proc(cc,2)
972 !.. Limit riming rate if too small
973 if((rimesum(cc)*i_rhoair(k)).lt.QSMALL) then
974 rimesum(cc) =0.
975 qi_qc_nrm(cc)=0.
976 qi_qc_nrd(cc)=0.
977 endif
978 endif !.. End ice-cloud riming
980 !--------------------------------------------------------------------------------------------------------------!
981 !--------------------------------------------------------------------------------------------------------------!
982 !.. Ice-rain riming rates from lookup table
983 if(RAIN_ICE.and.qr(k).gt.QSMALL) then
985 nr(k) = max(nr(k),QNSMALL)
986 !.. Get rain size distribution properties
987 lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333
988 n0rr = nr(k)*lamr
989 if(lamr.LT.LAMMINR) then
990 lamr = LAMMINR
991 n0rr = lamr**4*qr(k)/(PI*RHOW)
992 nr(k) = n0rr/lamr
993 elseif(lamr.gt.LAMMAXR) then
994 lamr = LAMMAXR
995 n0rr = lamr**4*qr(k)/(PI*RHOW)
996 nr(k) = n0rr/lamr
997 endif
998 rrr = 0.5 * (1./lamr)
1000 !.. The lookup table is in log space. Do not change the following 4 lines without
1001 !.. changing the offline-built lookup table
1002 rrni = 13.498*log10(0.5e6*rni(cc))
1003 rrri= 23.273*log10(1.e5*rrr)
1004 rdsi = 50.*(deltastr(cc) - 0.5)
1005 rrho = 7.888*log10(0.02*rhobar(cc))
1007 irni = int(rrni)
1008 irri = int(rrri)
1009 idsi = int(rdsi)
1010 irho = int(rrho)
1012 !.. Small lookup table limit
1013 rrni = max(rrni,1.)
1014 rrri = max(rrri,1.)
1015 rdsi = max(rdsi,1.)
1016 rrho = max(rrho,1.)
1018 irni = max(irni,1)
1019 irri = max(irri,1)
1020 idsi = max(idsi,1)
1021 irho = max(irho,1)
1023 !.. Large lookup table limit
1024 rrni = min(rrni,real(size(itabr,1))-1.)
1025 rrri = min(rrri,real(size(itabr,2))-1.)
1026 rdsi = min(rdsi,real(size(itabr,3))-1.)
1027 rrho = min(rrho,real(size(itabr,4))-1.)
1029 irni = min(irni,size(itabr,1)-1)
1030 irri = min(irri,size(itabr,2)-1)
1031 idsi = min(idsi,size(itabr,3)-1)
1032 irho = min(irho,size(itabr,4)-1)
1034 do iti = 1, 6
1035 call access_lookup_table(itabr,irni,irri,idsi,irho,iti, &
1036 rdsi,rrho,rrri,rrni,procr(cc,iti))
1037 enddo
1039 !.. procr(cc,4) is the normalized ice number tendency from ice-rain
1040 !.. procr(cc,5) is the normalized rain mass tendency from ice-rain
1041 !.. procr(cc,6) is the normalized ice mass tendency from ice-rain
1042 numrateri(cc) = max((procr(cc,4)*ni(cc,k)*nr(k)*rhoair(k)),0.)
1043 rainrateri(cc)= max((procr(cc,5)*ni(cc,k)*nr(k)*rhoair(k)),0.)
1044 icerateri(cc) = max((procr(cc,6)*ni(cc,k)*nr(k)*rhoair(k)),0.)
1046 !.. Limit ice-rain rates if to small
1047 if(rainrateri(cc).lt.QSMALL.or.icerateri(cc).lt.QSMALL) then
1048 rainrateri(cc)=0.
1049 icerateri(cc) =0.
1050 numrateri(cc) =0.
1051 endif
1053 if(temp.le.T0) then
1054 !.. For temperatures less than melting point ice-rain freezes rain
1055 if(qr(k).gt.0.1e-3.and.qi(cc,k).gt.0.1e-3) then
1056 dQRfzri(cc) = rainrateri(cc)
1057 dQIfzri(cc) = icerateri(cc)
1058 dNfzri(cc) = numrateri(cc)
1059 endif
1060 else
1061 !.. For temperatures greater than melting point ice-rain melts ice
1062 dQImltri(cc) = icerateri(cc)
1063 dNmltri(cc) = numrateri(cc)
1064 endif
1066 !.. procr(cc,1) is the normalized riming rate
1067 !.. procr(cc,2) is trelated to the he normalized rime density
1068 !.. procr(cc,3) is the normalized loss in rain number from riming
1069 !.. rimesumr is in kg m^-3 s^-1
1071 !.. Ice-rain leads to riming when the fall speed of the collected rain is slower than
1072 !.. the fall speed of the ice collector. Otherwise, ice-rain freezes rain
1073 rimesumr(cc) = max((procr(cc,1)*ni(cc,k)*nr(k)*rhoair(k)**2),0.)
1074 qi_qr_nrm(cc)= procr(cc,1)
1075 qi_qr_nrd(cc)= procr(cc,2)
1076 qi_qr_nrn(cc)= procr(cc,3)
1078 !.. Limit ice-rain riming if small
1079 if((rimesumr(cc)*i_rhoair(k)).lt.QSMALL) then
1080 rimesumr(cc) =0.
1081 qi_qr_nrm(cc)=0.
1082 qi_qr_nrd(cc)=0.
1083 qi_qr_nrn(cc)=0.
1084 endif
1086 endif !.. End ice-rain collection
1088 !.. Do not over-deplete cloud water or rain from ice-cloud and ice-rain
1089 tmpsum = tmpsum + rimesum(cc)
1090 tmpsumr= tmpsumr+ rimesumr(cc)
1092 if(ccvar.eq.numice) then !.. Check at end
1094 !.. Do not let ice-cloud collection overdeplete cloud water
1095 sink = tmpsum*i_rhoair(k)
1096 source = qc(k)*i_dt
1097 if(sink.gt.source.and.qc(k).gt.QSMALL) then
1098 ratio = source / sink
1099 rimesum = rimesum*ratio
1100 qi_qc_nrm= qi_qc_nrm*ratio
1101 qi_qc_nrd= qi_qc_nrd*ratio
1102 endif
1104 !.. Do not let ice-rain collection overdeplete rainwater
1105 sink = tmpsum*i_rhoair(k)
1106 source = qr(k)*i_dt
1107 if(sink.gt.source.and.qr(k).gt.QSMALL) then
1108 ratio = source / sink
1109 rimesumr = rimesumr*ratio
1110 qi_qr_nrm = qi_qr_nrm*ratio
1111 qi_qr_nrd = qi_qr_nrd*ratio
1112 qi_qr_nrn = qi_qr_nrn*ratio
1113 endif
1114 endif
1116 !--------------------------------------------------------------------------------------------------------------!
1117 !--------------------------------------------------------------------------------------------------------------!
1118 !.. Inital volume, iwc, #/m^3 for vapor growth
1119 gamma_arg = NU+2.+deltastr(cc)
1120 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1122 vi = fourthirdspi*rni(cc)**3*gamma_tab(gi)*i_gammnu
1123 iwci = ni(cc,k)*rhobar(cc)*vi*rhoair(k)
1124 nim3dum=ni(cc,k)*rhoair(k)
1126 !.. Total riming rate (the heating goes into vapor growth)
1127 rimetotal(cc) = rimesum(cc)+rimesumr(cc)
1128 alphstr=ao**(1.-deltastr(cc))
1129 !.. Distribtuion-averaged capacitance (see Harrington et al. 2013)
1130 capgam = capacitance_gamma(ani(cc), deltastr(cc), NU, alphstr, i_gammnu)
1132 !.. This subroutine first calculates the bulk fall speeds, then the ventilation,
1133 !.. and finally the vapor growth (or sublimation) rate. (See Harrington et al. 2013)
1134 call vaporgrow(dt, ani(cc), cni(cc), rni(cc), igr, nim3dum, temp, rimetotal(cc), &
1135 pres_e(k), NU, alphstr, sui, sup, qvs, qvi, mu, iwci, rhoair(k), qi(cc,k), &
1136 dv, kt, ao, nsch, npr, gammnu, i_gammnu, fourthirdspi, svpi, xxls, xxlv, xxlf, &
1137 capgam, vtrni1(cc,k), vtrmi1(cc,k), vtrzi1(cc,k), anf, cnf, rnf, iwcf, fv(cc), &
1138 fh(cc), rhobar(cc), deltastr(cc),rhodepout(cc),dsdepout(cc))
1140 !.. Vapor growth mass growth and axes rates
1141 prd(cc)=(iwcf-iwci)*i_rhoair(k)*i_dt
1142 prd(cc)=max(prd(cc),-qi(cc,k)*i_dt)
1143 ard(cc)=(2.*(anf-ani(cc))*cni(cc)+(cnf-cni(cc))*ani(cc))*ani(cc)*ni(cc,k)*i_dt
1144 crd(cc)=(2.*(cnf-cni(cc))*ani(cc)+(anf-ani(cc))*cni(cc))*cni(cc)*ni(cc,k)*i_dt
1146 !.. Sublimation number loss
1147 if(prd(cc).lt.0..and.qi(cc,k).gt.QSMALL) then
1148 nrd(cc)=prd(cc)*ni(cc,k)/qi(cc,k)
1149 endif
1151 !.. Limit vapor growth if values small
1152 if(abs(prd(cc)*dt).lt.(QSMALL*0.01)) then
1153 prd(cc)=0.
1154 nrd(cc)=0.
1155 ard(cc)=0.
1156 crd(cc)=0.
1157 endif
1159 !--------------------------------------------------------------------------------------------------------------!
1160 !--------------------------------------------------------------------------------------------------------------!
1161 !.. Now calculated the mass growth and axes growth rates and density for riming
1163 !.. Check for wet-growth riming conditions
1164 !.. See Lamb and Verlinde (2011)
1165 if(WET_GROW) then
1166 rimetotal(cc)=rimesum(cc)+rimesumr(cc)
1167 call wet_growth_check(NU, temp, rhoair(k), xxlv, xxlf, qv(k), dv, kt, qs0, &
1168 fv(cc), fh(cc), rimetotal(cc), rni(cc), ni(cc,k), dry_growth(cc))
1169 endif
1171 !.. If above the melting point, assume high-density riming
1172 if(temp.gt.T0) then
1173 dry_growth(cc) = .false.
1174 endif
1176 !.. Initial rime volume and r-axis length
1177 vir = vi
1178 rnfr= rni(cc)
1180 !.. Calculate the ice-cloud rime density
1181 !.. This is parameterized from Macklin (1962) and is described in
1182 !.. Jensen et al. 2017 (JAS)
1183 if(RIMING.and.qc(k).gt.QSMALL.and.qi_qc_nrm(cc).gt.0.) then
1184 rimec1 = 0.0066 !.. Value at -10C
1185 if((temp-T0).lt.-30.) then
1186 rimec1 = 0.0036
1187 elseif((temp-T0).lt.-20.and.(temp-T0).ge.-30.) then
1188 dum = (abs((temp-T0))-20.) / 10.
1189 rimec1 = dum*(0.0036) + (1.-dum)*(0.004)
1190 elseif((temp-T0).lt.-15.and.(temp-T0).ge.-20.) then
1191 dum = (abs((temp-T0))-15.) / 5.
1192 rimec1 = dum*(0.004) + (1.-dum)*(0.005)
1193 elseif((temp-T0).lt.-10.and.(temp-T0).ge.-15.) then
1194 dum = (abs((temp-T0))-10.) / 5.
1195 rimec1 = dum*(0.005) + (1.-dum)*(0.0066)
1196 elseif((temp-T0).lt.-5.and.(temp-T0).ge.-10.) then
1197 dum = (abs((temp-T0))-5.) / 5.
1198 rimec1 = dum*(0.0066) + (1.-dum)*(0.012)
1199 elseif((temp-T0).ge.-5.) then
1200 rimec1 = 0.012
1201 endif
1202 gdenavg(cc) = (1000.*(0.8*tanh(rimec1*qi_qc_nrd(cc)/qi_qc_nrm(cc))+0.1))
1204 if((temp-T0).gt.-5..and.(temp-T0).le.0.) then
1205 dum = (abs((temp-T0))-0.) / 5.
1206 gdenavg(cc) = dum*gdenavg(cc) + (1.-dum)*900.
1207 endif
1209 if((temp-T0).gt.0..or..NOT.dry_growth(cc)) then
1210 gdenavg(cc)=900.
1211 endif
1212 !.. Keep rime density between 50 and 900 kg m^-3
1213 gdenavg(cc) = max(gdenavg(cc),50.)
1214 gdenavg(cc) = min(gdenavg(cc),900.)
1216 gamma_arg = NU+2.+deltastr(cc)
1217 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1219 !.. Change in r-axis from riming (based on the rime density)
1220 rimedr = max(((qi_qc_nrm(cc)/gdenavg(cc))/((gamma_tab(gi)*i_gammnu)* &
1221 4.*PI*rni(cc)**2)),0.)
1222 rnfr = rnfr + max((rimedr*nc(k)*rhoair(k)),0.)*dt
1224 endif !.. End ice-cloud rime density calculation
1226 !.. Calculate the ice-rain rime density
1227 !.. This is parameterized from Macklin (1962)
1228 if(RAIN_ICE.and.qr(k).gt.QSMALL.and.qi_qr_nrm(cc).gt.0.) then
1229 rimec1 = 0.0066 !.. Value at -10C
1230 if((temp-T0).lt.-30.) then
1231 rimec1 = 0.0036
1232 elseif((temp-T0).lt.-20.and.(temp-T0).ge.-30.) then
1233 dum = (abs((temp-T0))-20.) / 10.
1234 rimec1 = dum*(0.0036) + (1.-dum)*(0.004)
1235 elseif((temp-T0).lt.-15.and.(temp-T0).ge.-20.) then
1236 dum = (abs((temp-T0))-15.) / 5.
1237 rimec1 = dum*(0.004) + (1.-dum)*(0.005)
1238 elseif((temp-T0).lt.-10.and.(temp-T0).ge.-15.) then
1239 dum = (abs((temp-T0))-10.) / 5.
1240 rimec1 = dum*(0.005) + (1.-dum)*(0.0066)
1241 elseif((temp-T0).lt.-5.and.(temp-T0).ge.-10.) then
1242 dum = (abs((temp-T0))-5.) / 5.
1243 rimec1 = dum*(0.0066) + (1.-dum)*(0.012)
1244 elseif((temp-T0).ge.-5.) then
1245 rimec1 = 0.012
1246 endif
1247 gdenavgr(cc) = (1000.*(0.8*tanh(rimec1*qi_qr_nrd(cc)/qi_qr_nrm(cc))+0.1))
1249 if((temp-T0).gt.-5..and.(temp-T0).le.0.) then
1250 dum = (abs((temp-T0))-0.) / 5.
1251 gdenavgr(cc) = dum*gdenavgr(cc) + (1.-dum)*900.
1252 endif
1254 if((temp-T0).gt.0..or..NOT.dry_growth(cc)) then
1255 gdenavgr(cc)=900.
1256 endif
1257 !.. Keep rime density between 50 and 900 kg m^03
1258 gdenavgr(cc) = max(gdenavgr(cc),50.)
1259 gdenavgr(cc) = min(gdenavgr(cc),900.)
1261 gamma_arg = NU+2.+deltastr(cc)
1262 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1264 !.. Change in r-axs from riming (based on rime density)
1265 rimedrr = max(((qi_qr_nrm(cc)/gdenavgr(cc))/((gamma_tab(gi)*i_gammnu)* &
1266 4.*PI*rni(cc)**2)),0.)
1267 rnfr = rnfr + max((rimedrr*nr(k)*rhoair(k)),0.)*dt
1269 endif !.. End ice-rain rime density calculation
1271 !.. Caculate the volume after riming
1272 gamma_arg = NU+2.+deltastr(cc)
1273 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1274 vfr = fourthirdspi*rnfr**3*gamma_tab(gi)*i_gammnu
1275 vfr = max(vfr,vi)
1276 rnfr= max(rnfr,rni(cc))
1278 !.. Determine the fraction of rime mass from cloud versus small rain drops
1279 !.. Update the ice density and iwc
1280 if(rimetotal(cc).gt.0.and.vfr.gt.vi) then
1281 qcrimefrac(cc) = rimesum(cc)/rimetotal(cc)
1282 qcrimefrac(cc) = max(qcrimefrac(cc),0.)
1283 qcrimefrac(cc) = min(qcrimefrac(cc),1.)
1284 gdentotal = qcrimefrac(cc)*gdenavg(cc) + (1.-qcrimefrac(cc))*gdenavgr(cc)
1285 rhobarrime = rhobar(cc)
1286 rhorimeout(cc) = (rhobar(cc)*(vi/vfr))+(gdentotal*(1.-(vi/vfr)))
1287 rhorimeout(cc) = min(rhorimeout(cc),RHOI)
1288 iwcfr = rhorimeout(cc) * vfr * nim3dum
1289 rhorimeout(cc) = gdentotal
1290 rhorimeout(cc) = min(rhorimeout(cc),RHOI)
1291 else
1292 iwcfr= iwci
1293 rnfr = rni(cc)
1294 endif
1296 !.. For dry growth the drops freeze on contace and one axis (the axis perpendicular to the
1297 !.. fall direction) grows
1298 if(dry_growth(cc)) then
1299 cnfr =cni(cc)
1300 anfr =ani(cc)
1301 phibr= 1.
1302 if(rimetotal(cc).gt.0.and.vfr.gt.vi) then
1303 gamma_arg = NU-1.+deltastr(cc)
1304 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1305 phibr = cni(cc)/ani(cc)*gamma_tab(gi)*i_gammnu
1306 if(phibr.eq.1.) then !.. Spherical ice
1307 phifr= 1.
1308 anfr = rnfr
1309 cnfr = rnfr
1310 elseif(phibr.gt.1.25) then !.. Prolate ice
1311 phifr = phibr*(((rnfr/rni(cc))**3)**(-0.5))
1312 anfr = (ani(cc)/rni(cc)**(1.5))*rnfr**(1.5)
1313 elseif(phibr.lt.0.8) then !.. Oblate ice
1314 phifr = phibr*(rnfr/rni(cc))**3
1315 gamma_arg = NU-1.+deltastr(cc)
1316 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1317 cnfr = phifr*anfr*gammnu/gamma_tab(gi)
1318 else
1319 phifr = phibr !.. Constant aspect ratio if quasi-spherical
1320 gamma_arg = NU-1.+deltastr(cc)
1321 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1322 gamma_arg = NU+2.+deltastr(cc)
1323 gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1324 anfr = ((vfr*gamma_tab(gi))/(fourthirdspi*phifr*gamma_tab(gi2)))**0.3333333333
1325 cnfr = phifr*anfr*gammnu/gamma_tab(gi)
1326 endif
1328 !.. For long time steps, riming can cause oblate ice to become prolate or vise versa
1329 !.. This prevents that from occuring
1330 if(phibr.le.1..and.phifr.gt.1.) then
1331 phifr = 0.99
1332 gamma_arg = NU-1.+deltastr(cc)
1333 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1334 gamma_arg = NU+2.+deltastr(cc)
1335 gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1336 anfr = ((vfr*gamma_tab(gi))/(fourthirdspi*phifr*gamma_tab(gi2)))**0.3333333333
1337 cnfr = phifr*anfr*gammnu/gamma_tab(gi)
1338 endif
1340 if(phibr.ge.1..and.phifr.lt.1.) then
1341 phifr = 1.01
1342 gamma_arg = NU-1.+deltastr(cc)
1343 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1344 gamma_arg = NU+2.+deltastr(cc)
1345 gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1346 anfr = ((vfr*gamma_tab(gi))/(fourthirdspi*phifr*gamma_tab(gi2)))**0.3333333333
1347 cnfr = phifr*anfr*gammnu/gamma_tab(gi)
1348 endif
1349 endif !.. End axis evolution from riming
1351 !.. Don't let the axes shrink from riming
1352 cnfr=max(cnfr,cni(cc))
1353 anfr=max(anfr,ani(cc))
1355 !.. Riming mass growth and axes rates
1356 prdr(cc)=(iwcfr-iwci)*i_rhoair(k)*i_dt
1357 ardr(cc)=(2.*(anfr-ani(cc))*cni(cc)+(cnfr-cni(cc))*ani(cc))*ani(cc)*ni(cc,k)*i_dt
1358 crdr(cc)=(2.*(cnfr-cni(cc))*ani(cc)+(anfr-ani(cc))*cni(cc))*cni(cc)*ni(cc,k)*i_dt
1360 !.. These can never be negative
1361 prdr(cc) = max(prdr(cc),0.)
1362 ardr(cc) = max(ardr(cc),0.)
1363 crdr(cc) = max(crdr(cc),0.)
1365 !.. If no rime mass, no axis growth
1366 if(prdr(cc).eq.0.0) then
1367 ardr(cc) = 0.
1368 crdr(cc) = 0.
1369 endif
1371 else !.. Wet Growth
1373 !.. For wet growth, mass is added but soaks into the particle and thus
1374 !.. the density increases but the size does not
1375 prdr(cc)=(iwcfr-iwci)*i_rhoair(k)*i_dt
1376 ardr(cc)=0.
1377 crdr(cc)=0.
1379 endif !.. End dry versus wet growth
1381 if(temp.gt.T0) then
1383 !.. Melting: ice shapes become more spherical
1384 qmlt(cc)=2.*PI*(kt*fh(cc)*(T0-temp)+rhoair(k)*xxlv*dv*fv(cc)*(qs0-qv(k)))/ &
1385 xxlf*(ni(cc,k)*NU*max(ani(cc),cni(cc))) - &
1386 (CPW/xxlf*(temp-T0)*(rimetotal(cc)/rhoair(k) + dQImltri(cc)))
1388 qmlt(cc)=min(qmlt(cc),0.)
1389 qmlt(cc)=max(qmlt(cc),(-qi(cc,k)*i_dt))
1391 !.. Don't let small very number mixing ratios (large sizes) cause ice precipitation
1392 if(qmlt(cc).lt.0.) then
1393 if(ai(cc,k).lt.1.e-12.or.ci(cc,k).lt.1.e-12) then
1394 qmlt(cc)=-qi(cc,k)*i_dt
1395 endif
1396 endif
1398 nmlt(cc) = max((-ni(cc,k)*i_dt),(ni(cc,k)*qmlt(cc)/qi(cc,k)-dNmltri(cc)))
1400 gamma_arg = NU+2.+deltastr(cc)
1401 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1402 tmpmelt=fourthirdspi*alphstr*gamma_tab(gi)*i_gammnu
1404 amlt(cc) = (1./tmpmelt)*(1./rhobar(cc))*qmlt(cc) + ai(cc,k)*nmlt(cc)/ni(cc,k)
1405 cmlt(cc) = amlt(cc)*cni(cc)/ani(cc)*(1.+(2.*deltastr(cc)))/(2.+deltastr(cc))
1407 endif
1408 enddo !.. End the loop over ice speices
1409 endif !.. End if ice
1411 !--------------------------------------------------------------------------------------------------------------!
1412 !--------------------------------------------------------------------------------------------------------------!
1413 !.. Nucleation: this occurs outside of the ice species loop
1415 !.. Freeze evertying at temperatures below -35C
1416 !.. That is how it happens in nature, right?
1417 if(FREEZE_QC.and.qc(k).gt.QSMALL.and.temp.lt.(T0-35.)) then
1418 mim = qc(k)*i_dt
1419 nim = nc(k)*i_dt
1420 endif
1422 if(FREEZE_QC.and.qr(k).gt.QSMALL.and.temp.lt.(T0-35.)) then
1423 mimr = qr(k)*i_dt
1424 nimr = nr(k)*i_dt
1425 endif
1427 !.. Bigg (1953) freezing below -4C
1428 if(qr(k).gt.QSMALL.and.temp.lt.(T0-4.)) then
1429 nr(k) = max(nr(k),QNSMALL)
1430 lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333
1431 n0rr = nr(k)*lamr
1432 if(lamr.LT.LAMMINR) then
1433 lamr = LAMMINR
1434 n0rr = lamr**4*qr(k)/(PI*RHOW)
1435 nr(k) = n0rr/lamr
1436 elseif(lamr.gt.LAMMAXR) then
1437 lamr = LAMMAXR
1438 n0rr = lamr**4*qr(k)/(PI*RHOW)
1439 nr(k) = n0rr/lamr
1440 endif
1442 mbiggr = 20.*PI*PI*RHOW*100.*nr(k)*(exp(0.66*(T0-temp))-1.)/lamr**3/lamr**3
1443 nbiggr = PI*100.*nr(k)*(exp(0.66*(T0-temp))-1.)/lamr**3
1444 mbiggr = min(mbiggr,qr(k)*i_dt)
1445 nbiggr = min(nbiggr,nr(k)*i_dt)
1446 endif
1448 !.. DeMott (2010) PNAS, Heterogeneous nuclation
1449 if(temp.lt.T0.and.sup.ge.0.) then
1451 a_demott = 0.0000594
1452 b_demott = 3.33
1453 c_demott = 0.0264
1454 d_demott = 0.0033
1456 !.. The 0.03 is the number of large aerosol (>0.5 microns)
1457 !.. See Chagnon and Junge 1962
1458 inrate = a_demott * (273.16 - temp)**b_demott * &
1459 0.03**((c_demott*(273.16 - temp)) + d_demott)
1460 inrate = inrate*1000.0 !.. #/L to #/m^3
1461 inrate = inrate*i_rhoair(k) !.. # kg^-1
1463 if((((ni(ICE1,k)+ni(ICE2,k)+inrate)*rhoair(k))/1000.).le.10000.) then
1464 mnuccd=(inrate)*(fourthirdspi*RHOI*(2.e-6)**3)*i_dt
1465 nnuccd=inrate*i_dt
1466 else
1467 curnum = (ni(ICE1,k)+ni(ICE2,k))*rhoair(k)/1000.
1468 ratel = max(0.,(10000.-curnum))
1469 ratekg = ratel*1000.*i_rhoair(k)
1470 mnuccd=(ratekg)*(fourthirdspi*RHOI*(2.e-6)**3)*i_dt
1471 nnuccd=ratekg*i_dt
1472 endif
1473 endif
1475 !.. Rime splinter (Hallett and Mossop 1974, Nature)
1476 if(SPLINTERS) then
1477 if(temp.lt.270.16.and.temp.gt.265.16) then
1478 if (temp.GT.270.16) then
1479 fmult = 0.
1480 elseif(temp.le.270.16.and.temp.gt.268.16) then
1481 fmult = (270.16-temp)/2.
1482 elseif(temp.ge.265.16.and.temp.le.268.16) then
1483 fmult = (temp-265.16)/3.
1484 elseif(temp.lt.265.16) then
1485 fmult = 0.
1486 endif
1487 endif
1489 do cc = 1, cat
1490 if(prdr(cc).gt.0.) then
1491 nmult(cc) = 35.E4*prdr(cc)*fmult*1000.
1492 qmult(cc) = nmult(cc)*fourthirdspi*RHOI*(5.e-6)**3
1493 qmult(cc) = min(qmult(cc),prdr(cc))
1494 prdr(cc) = prdr(cc)-qmult(cc)
1495 endif
1496 enddo
1497 endif
1499 !.. Aggregation (currently only for T < T0)
1500 if(temp.le.T0) then
1502 do cc = 1, cat
1503 ni(cc,k) = max(ni(cc,k),QNSMALL)
1504 ai(cc,k) = max(ai(cc,k),QASMALL)
1505 ci(cc,k) = max(ci(cc,k),QASMALL)
1506 enddo
1508 !.. Get characteristic diameters for aggregation
1509 dn1 = 2.*((ai(ICE1,k)**2)/(ci(ICE1,k)*ni(ICE1,k)))**0.333333333333
1510 dn2 = 2.*((ci(ICE2,k)**2)/(ai(ICE2,k)*ni(ICE2,k)))**0.333333333333
1511 dn3 = 2.*((ai(ICE3,k)**2)/(ci(ICE3,k)*ni(ICE3,k)))**0.333333333333
1513 dn1 = MIN(dn1,1.e-2)
1514 dn2 = MIN(dn2,1.e-2)
1515 dn3 = MIN(dn3,1.e-2)
1517 dn1 = MAX(dn1,1.e-6)
1518 dn2 = MAX(dn2,1.e-6)
1519 dn3 = MAX(dn3,1.e-6)
1521 !.. Get ice-one and ice-two aspect ratios
1522 !.. Spherical paritcles aggregate at a slower rate than less spherical ones
1523 gamma_arg = NU-1.+deltastr(1)
1524 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1525 phiagg1=ci(1,k)/ai(1,k)*gamma_tab(gi)*i_gammnu
1527 gamma_arg = NU-1.+deltastr(2)
1528 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1529 phiagg2=ci(2,k)/ai(2,k)*gamma_tab(gi)*i_gammnu
1531 phiagg1=MIN(phiagg1,100.)
1532 phiagg1=MAX(phiagg1,0.01)
1534 phiagg2=MIN(phiagg2,100.)
1535 phiagg2=MAX(phiagg2,0.01)
1537 !.. Ice-1 and ice-2 mass and number loss to aggregates occurs inside of aggregation
1538 !.. After aggregation, ice-1 and ice-2 bulk distribution properties are adjusted
1539 !.. assuming that aggregation did not cause the original shape and density of the
1540 !.. unaggregated particles to change
1541 if(qi(ICE1,k).gt.1.e-8.or.qi(ICE2,k).gt.1.e-8.or.qi(ICE3,k).gt.1.e-8) then
1542 call aggregation(dt,rhoair(k),temp,qi(ICE1,k),ni(ICE1,k),dn1,qi(ICE2,k),ni(ICE2,k),dn2, &
1543 qi(ICE3,k),ni(ICE3,k),dn3,rhobar(ICE1),rhobar(ICE2),phiagg1,phiagg2,coltab, &
1544 coltabn,qagg(ICE1),qagg(ICE2),qagg(ICE3),nagg(ICE1),nagg(ICE2),nagg(ICE3))
1545 endif !.. If ice mass mixing ratios are high enough to aggregate
1546 endif
1548 !--------------------------------------------------------------------------------------------------------------!
1549 !--------------------------------------------------------------------------------------------------------------!
1550 !.. Warm Microphysics
1552 !.. Get raindrop size distribution parameters
1553 if (qr(k).gt.QSMALL) then
1554 nr(k) = max(nr(k),QNSMALL)
1555 lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333
1556 n0rr = nr(k)*lamr
1557 if(lamr.LT.LAMMINR) then
1558 lamr = LAMMINR
1559 n0rr = lamr**4*qr(k)/(PI*RHOW)
1560 nr(k) = n0rr/lamr
1561 elseif(lamr.gt.LAMMAXR) then
1562 lamr = LAMMAXR
1563 n0rr = lamr**4*qr(k)/(PI*RHOW)
1564 nr(k) = n0rr/lamr
1565 endif
1566 endif
1568 !.. note: nprc1 is change in Nr,
1569 !.. nprc is change in Nc
1570 !.. NOTE: NPRC IS CALCULATED BUT NOT CURRENTLY USED SINCE NC IS NOT PROGNOSED
1571 !.. Auto conversion of cloud water to rain (Khairoutdinov and Kogan 2000)
1572 if(qc(k).ge.1.E-6) then
1573 PRC=1350.*qc(k)**2.47* &
1574 (nc(k)/1.e6*rhoair(k))**(-1.79)
1575 NPRC1 = PRC/(4./3.*PI*RHOW*(25.E-6)**3)
1576 NPRC = PRC/(qc(k)/nc(k))
1577 NPRC = min(NPRC,nc(k)*i_dt)
1578 NPRC1 = min(NPRC1,NPRC)
1579 endif
1581 !.. Accretion of cloud water by rain Khairoutdinov and Kogan 2000)
1582 !.. NOTE: NPRA CURRENTLY NOT USED SINCE NC IS NOT PROGNOSED
1583 if(qr(k).ge.1.e-8.and.qc(k).ge.1.e-8) then
1584 DUM=(qc(k)*qr(k))
1585 PRA = 67.*(DUM)**1.15
1586 NPRA = PRA/(qc(k)/nc(k))
1587 endif
1589 !.. Self-collection (Beheng 1994), Breakup (Verlinde and Cotton 1993)
1590 if(qr(k).ge.1.e-8) then
1591 dum1=300.e-6
1592 if(1./lamr.lt.dum1) then
1593 dum=1.
1594 elseif(1./lamr.ge.dum1) then
1595 dum=2.-exp(2300.*(1./lamr-dum1))
1596 endif
1597 nragg = -5.78*dum*nr(k)*qr(k)*rhoair(k)
1598 endif
1600 !.. Evaporation of rain (Rutledge and Hobbs 1983)
1601 if(qr(k).gt.QSMALL) then
1602 epsr = 2.*PI*n0rr*rhoair(k)*dv* &
1603 (F1R/(LAMR*LAMR)+ &
1604 F2R*(arn*rhoair(k)/mu)**0.5* &
1605 nsch**0.333333333*(gamma(5./2.+BR/2.))/ &
1606 (lamr**(5./2.+BR/2.)))
1607 else
1608 epsr = 0.
1609 endif
1611 if(qv(k).lt.qvs) then
1612 pre = epsr*(qv(k)-qvs)/ab
1613 pre = MIN(PRE,0.)
1614 dumt = temp+xxlv*i_cp*pre*dt
1615 dumqv = qv(k)-pre*dt
1616 dumqvs = 0.622*polysvp(dumt,0)/ &
1617 (pres_e(k)-polysvp(dumt,0))
1618 if(pre.lt.0..and.dumqv/dumqvs.gt.1.) then
1619 pre=(qv(k)-qvs)*i_dt/ab
1620 endif
1621 else
1622 pre = 0.
1623 endif
1625 !.. Rain number mixing ratio loss from evaporation
1626 if(pre.lt.0.) then
1627 dum = pre*dt/qr(k)
1628 dum = max(-1.,dum)
1629 npre = dum*nr(k)/dt
1630 endif
1632 !--------------------------------------------------------------------------------------------------------------!
1633 !--------------------------------------------------------------------------------------------------------------!
1634 !.. Updates after process-rate calculations
1635 !.. For ice, updates to mass, number, volume, and volume
1636 !.. times aspect ratio occur first for vapor growth,
1637 !.. riming, and melting (shape is predicted to evolve).
1638 !.. Next updates the bulk distribution shapes are done for
1639 !.. nucleation and aggregation. For nucleation, ice particle shape
1640 !.. itself does not change but the distribution shape becomes
1641 !.. more spherical as nucleated particles are added.
1643 !.. Do not over-deplete
1645 !.. Cloud water
1646 !.. First sum the total loss from riming
1647 qcrimesum=0.
1648 do cc = 1, cat
1649 qcrimesum = qcrimesum + prdr(cc)*qcrimefrac(cc)
1650 enddo
1652 sink =(prc+pra+qcrimesum+mim)*dt
1653 source=qc(k)
1654 if(sink.gt.source.and.qc(k).gt.QSMALL) then
1655 ratio = source/sink
1656 prc = prc*ratio
1657 pra = pra*ratio
1658 mim = mim*ratio
1659 do cc = 1, cat
1660 prdr(cc) = prdr(cc)*ratio
1661 ardr(cc) = ardr(cc)*ratio
1662 crdr(cc) = crdr(cc)*ratio
1663 enddo
1664 endif
1666 !.. Rain
1667 !.. First sum the total loss from riming
1668 qcrimesumr=0.
1669 do cc = 1, cat
1670 qcrimesumr = qcrimesumr + prdr(cc)*(1.-qcrimefrac(cc))
1671 enddo
1673 sink = (-pre+qcrimesumr+mimr+mbiggr)*dt
1674 source = qr(k) + (pra+prc)*dt
1675 do cc = 1, cat
1676 source = source + (-qmlt(cc)*dt)
1677 enddo
1678 do cc = 1, cat
1679 sink = sink + (dQRfzri(cc))*dt
1680 enddo
1681 if(sink.gt.source.and.qr(k).gt.QSMALL) then
1682 ratio= source/sink
1683 do cc = 1, cat
1684 prdr(cc) = prdr(cc)*ratio
1685 ardr(cc) = ardr(cc)*ratio
1686 crdr(cc) = crdr(cc)*ratio
1687 qi_qr_nrn(cc)=qi_qr_nrn(cc)*ratio
1688 dQRfzri(cc) = dQRfzri(cc)*ratio
1689 dQIfzri(cc) = dQIfzri(cc)*ratio
1690 dNfzri(cc) = dNfzri(cc)*ratio
1691 enddo
1692 mimr = mimr*ratio
1693 nimr = nimr*ratio
1694 mbiggr =mbiggr*ratio
1695 nbiggr =nbiggr*ratio
1696 pre = pre*ratio
1697 endif
1699 !.. Ice
1700 if(igr.le.1.) then !.. Nucleation initiation to ice 1
1701 !.. Ice 1
1702 sink =(-qmlt(ICE1)-qagg(ICE1))*dt
1703 source=qi(ICE1,k) + (prdr(ICE1))*dt
1704 if(prd(ICE1).gt.0.) then
1705 source = source + prd(ICE1)*dt
1706 else
1707 sink = sink + (-prd(ICE1))*dt
1708 endif
1710 !.. Added because of ice initiation to ice 1
1711 source = source + (mnuccd+mim+mimr+mbiggr+dQIfzri(ICE2)+dQIfzri(ICE3))*dt
1712 do cc = 1, cat
1713 source = source + (qmult(cc)+dQRfzri(cc))*dt
1714 enddo
1716 if(sink.gt.source.and.qi(ICE1,k).gt.QSMALL) then
1717 ratio= source/sink
1718 qmlt(ICE1) = qmlt(ICE1)*ratio
1719 qagg(ICE1) = qagg(ICE1)*ratio
1720 nagg(ICE1) = nagg(ICE1)*ratio
1721 if(prd(ICE1).lt.0.) then
1722 prd(ICE1) = prd(ICE1)*ratio
1723 endif
1724 endif
1726 !.. Ice 2
1727 sink =(dQIfzri(ICE2)-qmlt(ICE2)-qagg(ICE2))*dt
1728 source=qi(ICE2,k) + (prdr(ICE2))*dt
1729 if(prd(ICE2).gt.0.) then
1730 source = source + prd(ICE2)*dt
1731 else
1732 sink = sink + (-prd(ICE2))*dt
1733 endif
1735 if(sink.gt.source.and.qi(ICE2,k).gt.QSMALL) then
1736 ratio= source/sink
1737 dQIfzri(ICE2) = dQIfzri(ICE2)*ratio
1738 dQRfzri(ICE2) = dQRfzri(ICE2)*ratio
1739 dNfzri(ICE2) = dNfzri(ICE2)*ratio
1740 qmlt(ICE2) = qmlt(ICE2)*ratio
1741 qagg(ICE2) = qagg(ICE2)*ratio
1742 nagg(ICE2) = nagg(ICE2)*ratio
1743 if(prd(ICE2).lt.0.) then
1744 prd(ICE2) = prd(ICE2)*ratio
1745 endif
1746 endif
1747 !.. Ice 3
1748 sink =(dQIfzri(ICE3)-qmlt(ICE3))*dt
1749 source=qi(ICE3,k) + (prdr(ICE3)+qagg(ICE3))*dt
1750 if(prd(ICE3).gt.0.) then
1751 source = source + prd(ICE3)*dt
1752 else
1753 sink = sink + (-prd(ICE3))*dt
1754 endif
1756 if(sink.gt.source.and.qi(ICE3,k).gt.QSMALL) then
1757 ratio= source/sink
1758 dQIfzri(ICE3) = dQIfzri(ICE3)*ratio
1759 dQRfzri(ICE3) = dQRfzri(ICE3)*ratio
1760 dNfzri(ICE3) = dNfzri(ICE3)*ratio
1761 qmlt(ICE3) = qmlt(ICE3)*ratio
1762 if(prd(ICE3).lt.0.) then
1763 prd(ICE3) = prd(ICE3)*ratio
1764 endif
1765 endif
1766 else !.. nucleation to ice 2
1767 !.. Ice 2
1768 sink =(-qmlt(ICE2)-qagg(ICE2))*dt
1769 source=qi(ICE2,k) + (prdr(ICE2))*dt
1770 if(prd(ICE2).gt.0.) then
1771 source = source + prd(ICE2)*dt
1772 else
1773 sink = sink + (-prd(ICE2))*dt
1774 endif
1776 !.. Added because of ice initiation to ice 2
1777 source = source + (mnuccd+mim+mimr+mbiggr+dQIfzri(ICE1)+dQIfzri(ICE3))*dt
1778 do cc = 1, cat
1779 source = source + (qmult(cc)+dQRfzri(cc))*dt
1780 enddo
1782 if(sink.gt.source.and.qi(ICE2,k).gt.QSMALL) then
1783 ratio= source/sink
1784 qmlt(ICE2) = qmlt(ICE2)*ratio
1785 qagg(ICE2) = qagg(ICE2)*ratio
1786 nagg(ICE2) = nagg(ICE2)*ratio
1787 if(prd(ICE2).lt.0.) then
1788 prd(ICE2) = prd(ICE2)*ratio
1789 endif
1790 endif
1792 !.. Ice 1
1793 sink =(dQIfzri(ICE1)-qmlt(ICE1)-qagg(ICE1))*dt
1794 source=qi(ICE1,k) + (prdr(ICE1))*dt
1795 if(prd(ICE1).gt.0.) then
1796 source = source + prd(ICE1)*dt
1797 else
1798 sink = sink + (-prd(ICE1))*dt
1799 endif
1801 if(sink.gt.source.and.qi(ICE1,k).gt.QSMALL) then
1802 ratio= source/sink
1803 dQIfzri(ICE1) = dQIfzri(ICE1)*ratio
1804 dQRfzri(ICE1) = dQRfzri(ICE1)*ratio
1805 dNfzri(ICE1) = dNfzri(ICE1)*ratio
1806 qmlt(ICE1) = qmlt(ICE1)*ratio
1807 qagg(ICE1) = qagg(ICE1)*ratio
1808 nagg(ICE1) = nagg(ICE1)*ratio
1809 if(prd(ICE1).lt.0.) then
1810 prd(ICE1) = prd(ICE1)*ratio
1811 endif
1812 endif
1813 !.. Ice 3
1814 sink =(dQIfzri(ICE3)-qmlt(ICE3))*dt
1815 source=qi(ICE3,k) + (prdr(ICE3)+qagg(ICE3))*dt
1816 if(prd(ICE3).gt.0.) then
1817 source = source + prd(ICE3)*dt
1818 else
1819 sink = sink + (-prd(ICE3))*dt
1820 endif
1822 if(sink.gt.source.and.qi(ICE3,k).gt.QSMALL) then
1823 ratio= source/sink
1824 dQIfzri(ICE3) = dQIfzri(ICE3)*ratio
1825 dQRfzri(ICE3) = dQRfzri(ICE3)*ratio
1826 dNfzri(ICE3) = dNfzri(ICE3)*ratio
1827 qmlt(ICE3) = qmlt(ICE3)*ratio
1828 if(prd(ICE3).lt.0.) then
1829 prd(ICE3) = prd(ICE3)*ratio
1830 endif
1831 endif
1832 endif
1834 !.. Do not let vapor growth of ice create subsaturated conditions
1835 prdsum=0.
1836 do cc = 1, cat
1837 prdsum = prdsum + prd(cc)
1838 enddo
1839 vgflag = .false.
1840 sink = prdsum
1841 source=0.99*sui*qvi/abi*i_dt
1842 if(sink.gt.source.and.sui.gt.0.) then
1843 vgflag = .true.
1844 ratio = source/sink
1845 do cc = 1, cat
1846 prd(cc) = prd(cc)*ratio
1847 enddo
1848 endif
1850 !.. Update to ice
1851 do cc = 1, cat
1853 !.. Vapor growth update
1854 if((qi(cc,k)+(prd(cc))*dt).gt.QSMALL.and.qi(cc,k).gt.QSMALL) then
1855 rhobar(cc) = rhobar(cc)*(qi(cc,k)/(qi(cc,k)+(prd(cc))*dt)) + &
1856 rhodepout(cc)*(1.-(qi(cc,k)/(qi(cc,k)+(prd(cc))*dt)))
1857 else
1858 rhobar(cc) = RHOI
1859 endif
1860 qi(cc,k)=qi(cc,k)+(prd(cc))*dt
1861 ni(cc,k)=ni(cc,k)+(nrd(cc))*dt
1862 ai(cc,k)=ai(cc,k)+(ard(cc))*dt
1863 ci(cc,k)=ci(cc,k)+(crd(cc))*dt
1865 ni(cc,k) = max(ni(cc,k),QNSMALL)
1866 ai(cc,k) = max(ai(cc,k),QASMALL)
1867 ci(cc,k) = max(ci(cc,k),QASMALL)
1869 !.. If vapor growth of ice is going to cause
1870 !.. subsaturated condtions, reduce vapor growth
1871 !.. rates (done above) and rediagnose the axes
1872 if(vgflag.and.qi(cc,k).gt.QSMALL) then
1873 alphstr=ao**(1.-deltastr(cc))
1874 alphv=fourthirdspi*alphstr
1875 betam=2.+deltastr(cc)
1876 gamma_arg = NU+2.+deltastr(cc)
1877 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
1878 ani(cc)=((qi(cc,k)*gammnu)/(rhobar(cc)*ni(cc,k)*alphv* &
1879 gamma_tab(gi)))**(1./betam)
1880 cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc)
1881 ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k)
1882 ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k)
1883 endif
1884 enddo
1886 !.. Update to ice
1887 do cc = 1, cat
1888 !.. Riming update
1889 if((qi(cc,k)+(prdr(cc))*dt).gt.QSMALL) then
1890 if(dry_growth(cc)) then
1891 rhobar(cc) = rhobar(cc)*(qi(cc,k)/(qi(cc,k)+(prdr(cc))*dt)) + &
1892 rhorimeout(cc)*(1.-(qi(cc,k)/(qi(cc,k)+(prdr(cc))*dt)))
1893 else
1894 rhobar(cc) = rhobar(cc)*(qi(cc,k)/(qi(cc,k)+(prdr(cc))*dt)) + &
1895 RHOW*(1.-(qi(cc,k)/(qi(cc,k)+(prdr(cc))*dt)))
1896 rhobar(cc) = max(rhobar(cc),RHOI)
1897 endif
1898 else
1899 rhobar(cc) = RHOI
1900 endif
1902 qi(cc,k)=qi(cc,k)+(prdr(cc))*dt
1903 ai(cc,k)=ai(cc,k)+(ardr(cc))*dt
1904 ci(cc,k)=ci(cc,k)+(crdr(cc))*dt
1906 enddo
1908 !.. Keep everything in bounds after vapor growth and riming
1909 do cc = 1, cat
1910 if(qi(cc,k).gt.QSMALL) then
1911 ni(cc,k)=max(ni(cc,k),QNSMALL)
1912 ai(cc,k)=max(ai(cc,k),QASMALL)
1913 ci(cc,k)=max(ci(cc,k),QASMALL)
1914 ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
1915 cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
1917 ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k)
1918 ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k)
1920 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
1922 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
1923 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
1924 alphstr, alphv, betam)
1926 dsold(cc) = deltastr(cc)
1927 endif
1928 enddo
1930 !--------------------------------------------------------------------------------------------------------------!
1931 !--------------------------------------------------------------------------------------------------------------!
1932 !.. Update to ice after melting
1933 if(temp.gt.T0) then
1934 do cc = 1, cat
1936 if(qi(cc,k).gt.QSMALL) then
1937 rhobar(cc) = rhobar(cc)*((qi(cc,k)+(qmlt(cc))*dt)/qi(cc,k)) + &
1938 RHOW*((-qmlt(cc)*dt)/qi(cc,k))
1939 rhobar(cc) = min(rhobar(cc),RHOI)
1940 else
1941 rhobar(cc) = RHOI
1942 endif
1944 qi(cc,k)=qi(cc,k)+qmlt(cc)*dt
1945 ni(cc,k)=ni(cc,k)+nmlt(cc)*dt
1946 ai(cc,k)=ai(cc,k)+amlt(cc)*dt
1947 ci(cc,k)=ci(cc,k)+cmlt(cc)*dt
1949 !.. Small masses of leftover ice can become rain
1950 if(qi(cc,k).le.QSMALL) then
1951 qr(k) = qr(k) + qi(cc,k)
1952 qi(cc,k)=0.
1953 endif
1955 !.. Keep everyting in bounds after melting
1956 if(qi(cc,k).gt.QSMALL) then
1957 ni(cc,k)=max(ni(cc,k),QNSMALL)
1958 ai(cc,k)=max(ai(cc,k),QASMALL)
1959 ci(cc,k)=max(ci(cc,k),QASMALL)
1960 ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
1961 cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
1963 ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k)
1964 ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k)
1966 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
1968 !.. Melting should only make the average distribution more spherical
1969 if(dsold(cc).le.1.) then
1970 if(deltastr(cc).lt.dsold(cc)) then
1971 deltastr(cc)=dsold(cc)
1972 endif
1973 if(deltastr(cc).gt.1.) then
1974 deltastr(cc)=1.
1975 endif
1976 else
1977 if(deltastr(cc).gt.dsold(cc)) then
1978 deltastr(cc)=dsold(cc)
1979 endif
1980 if(deltastr(cc).lt.1.) then
1981 deltastr(cc)=1.
1982 endif
1983 endif
1984 cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc)
1985 ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k)
1986 ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k)
1988 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
1989 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
1990 alphstr, alphv, betam)
1992 dsold(cc) = deltastr(cc)
1993 endif !.. End if ice exists after melting
1994 enddo !.. Loop over all ices species
1995 endif !.. End temperature above melting
1997 !--------------------------------------------------------------------------------------------------------------!
1998 !--------------------------------------------------------------------------------------------------------------!
1999 !.. Update to vapor
2000 qv(k)=qv(k)+(-pre-mnuccd)*dt
2001 do cc = 1, cat
2002 qv(k)=qv(k)+(-prd(cc))*dt
2003 enddo
2005 !.. Update to cloud
2006 qc(k)=qc(k)+(-prc-pra-mim)*dt
2007 do cc = 1, cat
2008 qc(k) = qc(k)-(prdr(cc)*qcrimefrac(cc)*dt)
2009 enddo
2011 !.. Update to rain
2012 qr(k)=qr(k)+(pre+prc+pra-mimr-mbiggr)*dt
2013 nr(k)=nr(k)+(nprc1+nragg+npre-nimr-nbiggr)*dt
2014 do cc = 1, cat
2015 qr(k)=qr(k)+((-prdr(cc)*(1.-qcrimefrac(cc)))-qmlt(cc)-dQRfzri(cc))*dt
2016 nr(k)=max((nr(k)+((-qi_qr_nrn(cc)*ni(cc,k)*nr(k)*rhoair(k))-nmlt(cc)-dNfzri(cc))*dt),0.)
2017 enddo
2019 !--------------------------------------------------------------------------------------------------------------!
2020 !--------------------------------------------------------------------------------------------------------------!
2021 !.. Store the number before nucleation for later use
2022 do cc = 1, cat
2023 nibnuc(cc)=ni(cc,k)
2024 enddo
2026 if(igr.le.1.) then
2027 !.. Nucleation to ice1 (planar-nuclated)
2028 totalnuc(ICE1) =(mnuccd+mim+mimr+mbiggr+dQIfzri(ICE2)+dQIfzri(ICE3))*dt
2029 totalnucn(ICE1)=(nnuccd+nim+nimr+nbiggr+dNfzri(ICE2)+dNfzri(ICE3))*dt
2031 do cc = 1, cat
2032 totalnuc(ICE1) =totalnuc(ICE1) + (qmult(cc)+dQRfzri(cc))*dt
2033 totalnucn(ICE1)=totalnucn(ICE1) + (nmult(cc))*dt
2034 enddo
2036 qi(ICE1,k)=qi(ICE1,k)+totalnuc(ICE1)
2037 ni(ICE1,k)=ni(ICE1,k)+totalnucn(ICE1)
2039 !.. Loss of ice2 from ice2-rain producing ice1
2040 qi(ICE2,k)=qi(ICE2,k)-(dQIfzri(ICE2))*dt
2041 ni(ICE2,k)=ni(ICE2,k)-(dNfzri(ICE2))*dt
2043 !.. Loss of ice3 from ice3-rain producing ice1
2044 qi(ICE3,k)=qi(ICE3,k)-(dQIfzri(ICE3))*dt
2045 ni(ICE3,k)=ni(ICE3,k)-(dNfzri(ICE3))*dt
2046 else
2048 !.. Nucleation to ice2 (columnar-nucleated)
2049 totalnuc(ICE2) =(mnuccd+mim+mimr+mbiggr+dQIfzri(ICE1)+dQIfzri(ICE3))*dt
2050 totalnucn(ICE2)=(nnuccd+nim+nimr+nbiggr+dNfzri(ICE1)+dNfzri(ICE3))*dt
2052 do cc = 1, cat
2053 totalnuc(ICE2) =totalnuc(ICE2) + (qmult(cc)+dQRfzri(cc))*dt
2054 totalnucn(ICE2)=totalnucn(ICE2) + (nmult(cc))*dt
2055 enddo
2057 qi(ICE2,k)=qi(ICE2,k)+totalnuc(ICE2)
2058 ni(ICE2,k)=ni(ICE2,k)+totalnucn(ICE2)
2060 !.. Loss of ice1 from ice1-rain producing ice2
2061 qi(ICE1,k)=qi(ICE1,k)-(dQIfzri(ICE1))*dt
2062 ni(ICE1,k)=ni(ICE1,k)-(dNfzri(ICE1))*dt
2064 !.. Loss of ice3 from ice3-rain producing ice2
2065 qi(ICE3,k)=qi(ICE3,k)-(dQIfzri(ICE3))*dt
2066 ni(ICE3,k)=ni(ICE3,k)-(dNfzri(ICE3))*dt
2068 endif
2070 !.. Limit ice number concentrations to 1000/L
2071 do cc = 1, cat
2072 ni(cc,k)=min(ni(cc,k),(1000.*1000.*i_rhoair(k)))
2073 enddo
2075 !.. Final check on rain and fall speed calculation
2076 if (qr(k).gt.QSMALL) then
2077 nr(k) = max(nr(k),QNSMALL)
2078 lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333
2079 n0rr = nr(k)*lamr
2080 if(lamr.LT.LAMMINR) then
2081 lamr = LAMMINR
2082 n0rr = lamr**4*qr(k)/(PI*RHOW)
2083 nr(k) = n0rr/lamr
2084 elseif(lamr.gt.LAMMAXR) then
2085 lamr = LAMMAXR
2086 n0rr = lamr**4*qr(k)/(PI*RHOW)
2087 nr(k) = n0rr/lamr
2088 endif
2089 vtrn(k) = ARN*(gamma(1.+BR))/LAMR**BR
2090 vtrm(k) = ARN*(gamma(4.+BR))/6./LAMR**BR
2091 vtrn(k) = min(vtrn(k),9.1)
2092 vtrm(k) = min(vtrm(k),9.1)
2094 dbzr = rhoair(k)*nr(k)/lamr**3/lamr**3*720.
2095 dbzr = 1.e18*dbzr
2097 else
2098 !.. If only small mass of rain, push to vapor
2099 qv(k)=qv(k)+qr(k)
2100 theta(k)=theta(k)+theta(k)*i_temp*(qr(k)*xxlv)*i_cp
2101 qr(k) =0.
2102 nr(k) =0.
2103 vtrm(k) =0.
2104 vtrn(k) =0.
2105 endif
2107 !--------------------------------------------------------------------------------------------------------------!
2108 !--------------------------------------------------------------------------------------------------------------!
2109 do cc = 1, cat
2110 !.. Check for ice after nucleation and then update density and bulk distribution sizes
2111 if(qi(cc,k).ge.QSMALL) then
2112 ni(cc,k) =max(ni(cc,k),QNSMALL)
2113 nucfrac = totalnuc(cc)/qi(cc,k)
2114 nucfrac = min(nucfrac, 1.0)
2115 nucfrac = max(nucfrac, 0.0)
2117 !.. Do not nucleate to aggregates
2118 if(cc.eq.ICE3) then
2119 nucfrac=0.
2120 endif
2122 !.. Update density from nucleation
2123 rhobar(cc) = RHOI*nucfrac + (1.-nucfrac)*rhobar(cc)
2125 !.. Keep everything in bounds
2126 ai(cc,k) =max(ai(cc,k),QASMALL)
2127 ci(cc,k) =max(ci(cc,k),QASMALL)
2128 nibnuc(cc)=max(nibnuc(cc),QNSMALL)
2130 ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*nibnuc(cc)))**0.333333333333),2.e-6)
2131 cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*nibnuc(cc)))**0.333333333333),2.e-6)
2133 ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k)
2134 ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k)
2136 aniold(cc) = ani(cc)
2137 cniold(cc) =ao**(1.-dsold(cc))*aniold(cc)**dsold(cc)
2139 !.. Get the size distribution properties of nucleated particles
2140 if(nucfrac.gt.1.e-5.and.totalnucn(cc).gt.0.) then
2141 gamma_arg = NU+3.
2142 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2143 gamma_arg = NU+4.
2144 gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2145 anuc=(((totalnuc(cc)*gammnu)/(RHOI*totalnucn(cc)*fourthirdspi* &
2146 gamma_tab(gi)))**0.3333333)*gamma_tab(gi2)/gamma_tab(gi)
2147 anuc = MAX(anuc,2.e-6)
2148 anuc = MIN(anuc,3.e-3)
2150 gamma_arg = NU+2.+deltastr(cc)
2151 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2152 gamma_arg = NU+3.+deltastr(cc)
2153 gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2154 gamma_arg = NU+2.+(2.*deltastr(cc))
2155 gi3=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2157 amass=ani(cc)*gamma_tab(gi2)/gamma_tab(gi)
2158 cmass=cni(cc)*gamma_tab(gi3)/gamma_tab(gi)
2160 !.. Mass-weighted size of old distribution and nucleated distribution
2161 amass = amass*(1.-nucfrac) + anuc*nucfrac
2162 cmass = cmass*(1.-nucfrac) + anuc*nucfrac
2164 !.. New ani and cni after nucleation
2165 ani(cc) = max((amass*gamma_tab(gi)/gamma_tab(gi2)),2.e-6)
2166 cni(cc) = max((cmass*gamma_tab(gi)/gamma_tab(gi3)),2.e-6)
2168 !.. New delta* (shape) after nucleation
2169 if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then
2170 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
2171 else
2172 deltastr(cc)=1.
2173 endif
2175 !.. Nucleation should only make the average distribution more spherical
2176 if(dsold(cc).le.1.) then
2177 if(deltastr(cc).lt.dsold(cc)) then
2178 deltastr(cc)=dsold(cc)
2179 endif
2180 if(deltastr(cc).gt.1.) then
2181 deltastr(cc)=1.
2182 endif
2183 else
2184 if(deltastr(cc).gt.dsold(cc)) then
2185 deltastr(cc)=dsold(cc)
2186 endif
2187 if(deltastr(cc).lt.1.) then
2188 deltastr(cc)=1.
2189 endif
2190 endif
2191 endif !.. End nucleation update distribution properties
2193 !--------------------------------------------------------------------------------------------------------------!
2194 !--------------------------------------------------------------------------------------------------------------!
2195 !.. This needs to be done because of transfer from e.g.
2196 !.. ice1 to ice2 during nucleation.
2197 !.. Assume constant density/delta*, loss of mass and number
2199 !..With the new delta* and density after nucleation re-diagnose a and c
2200 !..Use the old density and deltastr for loss values e.g. q1 = q1 - (transfered to q2)
2202 alphstr=ao**(1.-deltastr(cc))
2203 alphv=fourthirdspi*alphstr
2204 betam=2.+deltastr(cc)
2205 gamma_arg = NU+2.+deltastr(cc)
2206 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2207 ani(cc)=((qi(cc,k)*gammnu)/(rhobar(cc)*ni(cc,k)*alphv* &
2208 gamma_tab(gi)))**(1./betam)
2209 cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc)
2211 !.. volume and volume*aspect ratio process rates
2212 ai(cc,k) = ai(cc,k) + cniold(cc)*nibnuc(cc)*2.*aniold(cc)*(ani(cc)-aniold(cc)) + &
2213 nibnuc(cc)*aniold(cc)**2*(cni(cc)-cniold(cc)) + &
2214 aniold(cc)**2*cniold(cc)*(ni(cc,k)-nibnuc(cc))
2216 ci(cc,k) = ci(cc,k) + aniold(cc)*nibnuc(cc)*2.*cniold(cc)*(cni(cc)-cniold(cc)) + &
2217 nibnuc(cc)*cniold(cc)**2*(ani(cc)-aniold(cc)) + &
2218 cniold(cc)**2*aniold(cc)*(ni(cc,k)-nibnuc(cc))
2220 !.. Keep everything in bounds
2221 ni(cc,k)=max(ni(cc,k),QNSMALL)
2222 ai(cc,k)=max(ai(cc,k),QASMALL)
2223 ci(cc,k)=max(ci(cc,k),QASMALL)
2225 if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then
2226 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
2227 else
2228 deltastr(cc)=1.
2229 endif
2231 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
2232 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
2233 alphstr, alphv, betam)
2235 dsold(cc) = deltastr(cc)
2237 !.. Limit the fall speeds here to 25 m s^-1
2238 vtrzi1(cc,k)=min(vtrzi1(cc,k),25.)
2239 vtrmi1(cc,k)=min(vtrmi1(cc,k),25.)
2240 vtrni1(cc,k)=min(vtrni1(cc,k),25.)
2241 else
2243 !.. If there is no ice even after nucleation,
2244 !.. add remainder to vapor and update temperature
2245 qv(k)=qv(k)+qi(cc,k)
2246 theta(k)=theta(k)+theta(k)*i_temp*(qi(cc,k)*xxls)*i_cp
2247 qi(cc,k)=0.
2248 ai(cc,k)=0.
2249 ni(cc,k)=0.
2250 ci(cc,k)=0.
2251 vtrzi1(cc,k)=0.
2252 vtrmi1(cc,k)=0.
2253 vtrni1(cc,k)=0.
2254 endif !.. End check for ice
2255 enddo !.. Loop over all ice species
2257 !--------------------------------------------------------------------------------------------------------------!
2258 !--------------------------------------------------------------------------------------------------------------!
2259 if(temp.le.T0) then
2260 !.. Aggregation occurs similarly to nulceation
2261 !.. We assume that the loss of ice1 and ice2 to aggregates does
2262 !.. not significantly change the shape and density of the
2263 !.. remaining particles. Both ani and cni are updated
2264 !.. with no change to density and delta*.
2265 !.. Aggregation (see Jensen et al. 2017, 2018a, 2018b)
2266 do cc = 1, cat
2267 dsold(cc) = deltastr(cc)
2268 nibnuc(cc) = max(ni(cc,k),QNSMALL)
2269 qi(cc,k) = qi(cc,k) + qagg(cc)
2270 ni(cc,k) = ni(cc,k) + nagg(cc)
2271 if(qi(cc,k).gt.QSMALL) then
2273 !.. Aggregate (unmelted) density and aspect ratio are assumed
2274 !.. Force the shape of aggregates
2275 if(cc.eq.ICE3) then
2276 rhobar(cc)=50.
2277 phiagg3 = 0.2
2278 ani(cc)=.5*dn3
2279 gamma_arg = NU-1.+deltastr(cc)
2280 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2281 cni(cc)=phiagg3*ani(cc)*gammnu/gamma_tab(gi)
2283 ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k)
2284 ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k)
2286 ni(cc,k)=max(ni(cc,k),QNSMALL)
2287 ai(cc,k)=max(ai(cc,k),QASMALL)
2288 ci(cc,k)=max(ci(cc,k),QASMALL)
2290 if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then
2291 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
2292 else
2293 deltastr(cc)=1.
2294 endif
2296 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
2297 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
2298 alphstr, alphv, betam)
2300 else
2302 !.. Variable check on ice-1 and ice-2 after aggregation
2303 ni(cc,k)=max(ni(cc,k),QNSMALL)
2304 ai(cc,k)=max(ai(cc,k),QASMALL)
2305 ci(cc,k)=max(ci(cc,k),QASMALL)
2306 ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*nibnuc(cc)))**0.333333333333),2.e-6)
2307 cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*nibnuc(cc)))**0.333333333333),2.e-6)
2309 aniold(cc) = ani(cc)
2310 cniold(cc) = cni(cc)
2311 deltastr(cc)= dsold(cc)
2313 !.. Update ice1 and ice2 sizes after aggregation from
2314 !.. updated mass and number mixiing ratios (transfer to aggregates) and an
2315 !.. assumed constant deltastar
2316 alphstr=ao**(1.-deltastr(cc))
2317 gamma_arg = NU+2.+deltastr(cc)
2318 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2319 ani(cc) = ((qi(cc,k))/(ni(cc,k)*rhobar(cc)*fourthirdspi*alphstr*gamma_tab(gi)*i_gammnu))**&
2320 (1./(2.+deltastr(cc)))
2321 ani(cc)=max(ani(cc),2.e-6)
2322 cni(cc)=alphstr*ani(cc)**deltastr(cc)
2324 !.. volume and volume*aspect ratio updates
2325 ai(cc,k) = ai(cc,k) + cniold(cc)*nibnuc(cc)*2.*aniold(cc)*(ani(cc)-aniold(cc)) + &
2326 nibnuc(cc)*aniold(cc)**2*(cni(cc)-cniold(cc)) + &
2327 aniold(cc)**2*cniold(cc)*(ni(cc,k)-nibnuc(cc))
2329 ci(cc,k) = ci(cc,k) + aniold(cc)*nibnuc(cc)*2.*cniold(cc)*(cni(cc)-cniold(cc)) + &
2330 nibnuc(cc)*cniold(cc)**2*(ani(cc)-aniold(cc)) + &
2331 cniold(cc)**2*aniold(cc)*(ni(cc,k)-nibnuc(cc))
2333 !.. Keep everything in bounds
2334 ni(cc,k)=max(ni(cc,k),QNSMALL)
2335 ai(cc,k)=max(ai(cc,k),QASMALL)
2336 ci(cc,k)=max(ci(cc,k),QASMALL)
2338 if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then
2339 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
2340 else
2341 deltastr(cc)=1.
2342 endif
2344 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
2345 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
2346 alphstr, alphv, betam)
2348 endif
2349 endif !.. If ice exist
2350 enddo !.. Loop over all species
2351 endif !.. If temperature < 273.15 K
2353 !--------------------------------------------------------------------------------------------------------------!
2354 !--------------------------------------------------------------------------------------------------------------!
2355 !.. Update the thermodynamics
2356 theta(k)=theta(k)+theta(k)*i_temp*((mnuccd)*xxls*i_cp*dt)
2357 do cc = 1, cat
2358 theta(k)=theta(k)+theta(k)*i_temp*((prd(cc))*xxls*i_cp*dt)
2359 theta(k)=theta(k)+theta(k)*i_temp*((qmlt(cc))*xxlf*i_cp*dt)
2360 enddo
2362 theta(k)=theta(k)+theta(k)*i_temp*(pre*xxlv)*i_cp*dt
2364 if(temp.le.T0) then
2365 theta(k)=theta(k)+theta(k)*i_temp*(mim+mimr+mbiggr)*xxlf*i_cp*dt
2367 do cc = 1, cat
2368 theta(k)=theta(k)+theta(k)*i_temp*(prdr(cc)+dQRfzri(cc))*xxlf*i_cp*dt
2369 enddo
2370 endif
2372 temp = theta(k)/(100000./pres_e(k))**(RCP)
2373 i_temp= 1./temp
2374 qvs = 0.622*polysvp(temp,0)/(pres_e(k)-polysvp(temp,0))
2375 xxlv = 3.1484E6-2370.*temp
2376 rhoair(k)= pres_e(k)/(RD*temp)
2377 i_rhoair(k) = 1./rhoair(k)
2379 !--------------------------------------------------------------------------------------------------------------!
2380 !--------------------------------------------------------------------------------------------------------------!
2381 !.. Saturation adjustment for qc
2382 pcc = (qv(k)-qvs) / (1.0 + xxlv**2*qvs/(CP*rv*temp**2))*i_dt
2383 if (pcc*dt+qc(k).lt.0.) then
2384 pcc=-qc(k)/dt
2385 end if
2386 theta(k)=theta(k) + theta(k)*i_temp*pcc*xxlv*i_cp*dt
2387 qv(k) = qv(k) - pcc*dt
2388 qc(k) = qc(k) + pcc*dt
2390 if(qc(k).le.QSMALL) then
2391 qv(k)=qv(k)+qc(k)
2392 theta(k)=theta(k)+theta(k)*i_temp*(qc(k)*xxlv)*i_cp
2393 qc(k)=0.
2394 nc(k)=0.
2395 endif
2397 !--------------------------------------------------------------------------------------------------------------!
2398 !--------------------------------------------------------------------------------------------------------------!
2399 !.. Cloud water effective radius
2400 lrsig = 1.3
2401 sig = log(lrsig) !.. Distribution standard deviation
2402 lwc = qc(k)*rhoair(k) !.. liquid water content
2403 ncm3dum = nc(k)*rhoair(k)
2405 if(ncm3dum .gt. 0.) then
2406 r_n = ((lwc)/(fourthirdspi*ncm3dum*RHOW*exp(4.5*(sig**2))))**(0.33333333333)
2407 else
2408 r_n = 0.
2409 endif
2411 effc1d(k) = r_n * (5.*EXP((sig**2)/2.))
2412 effc1d(k) = MIN(effc1d(k),50.e-6)
2413 effc1d(k) = MAX(effc1d(k),2.51e-6)
2415 ! For now, assume cloud droplet fallspeed is zero
2416 vtrmc(k) = 0.0
2417 vtrnc(k) = 0.0
2419 !--------------------------------------------------------------------------------------------------------------!
2420 !--------------------------------------------------------------------------------------------------------------!
2421 !.. Final check on ice before sedimentation
2422 do cc = 1, cat
2424 if(qi(cc,k).gt.QSMALL) then
2425 ni(cc,k)=max(ni(cc,k),QNSMALL)
2426 ai(cc,k)=max(ai(cc,k),QASMALL)
2427 ci(cc,k)=max(ci(cc,k),QASMALL)
2428 ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
2429 cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
2431 ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k)
2432 ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k)
2434 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
2436 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
2437 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
2438 alphstr, alphv, betam)
2440 !.. Final check on aggregates
2441 if(cc.eq.ICE3) then
2442 if(temp.le.T0) then
2443 rhobar(cc)=50.
2444 phiagg3 = 0.2
2445 gamma_arg = NU-1.+deltastr(cc)
2446 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2447 cni(cc)=phiagg3*ani(cc)*gammnu/gamma_tab(gi)
2449 ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k)
2450 ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k)
2451 endif
2453 if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then
2454 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
2455 else
2456 deltastr(cc)=1.
2457 endif
2458 !.. Keep aggregates spherical or oblate (for now)
2459 deltastr(cc) = min(deltastr(cc),1.)
2461 alphstr=ao**(1.-deltastr(cc))
2462 gamma_arg = NU+2.+deltastr(cc)
2463 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2464 ani(cc) = ((qi(cc,k))/(ni(cc,k)*rhobar(cc)*fourthirdspi*alphstr*&
2465 gamma_tab(gi)/gammnu))**(1./(2.+deltastr(cc)))
2467 ani(cc) = max(ani(cc),2.e-6)
2468 cni(cc)=alphstr*ani(cc)**deltastr(cc)
2470 !.. Final check on aggregate size (implicit breakup)
2471 if(ani(cc).gt.0.5e-3) then
2472 ani(cc)=0.5e-3
2473 cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc)
2474 gamma_arg = NU+2.+deltastr(cc)
2475 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2476 ni(cc,k) = (qi(cc,k))/(rhobar(cc)*fourthirdspi*ao**(1.-deltastr(cc))* &
2477 ani(cc)**(2.+deltastr(cc))*gamma_tab(gi)/gammnu)
2478 ni(cc,k)=max(ni(cc,k),QNSMALL)
2479 endif
2481 ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k)
2482 ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k)
2484 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
2485 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
2486 alphstr, alphv, betam)
2488 endif ! Check on aggregates
2490 !.. Fall speed moved to here
2491 !.. See Harrington et al. 2013 and Mitchell and Heymsfield 2005
2492 fv(cc) = 1.0
2493 fh(cc) = 1.0
2495 gamma_arg = NU-1.+deltastr(cc)
2496 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2497 phiivt = cni(cc)/ani(cc)*gamma_tab(gi)*i_gammnu
2498 alphstr=ao**(1.-deltastr(cc))
2500 !.. Determine coefficeients for the Best number for fall speeds
2501 if(phiivt.lt.1.0)then
2502 bl = 1.
2503 al = 2.
2504 aa = PI
2505 ba = 2.
2506 qe = (1.-phiivt)*(rhobar(cc)/RHOI) + phiivt
2507 else if(phiivt .gt. 1.0)then
2508 al = 2.0
2509 bl = 1.
2510 aa = PI*alphstr
2511 ba = deltastr(cc) + 1.
2512 qe = 1.0
2513 else if(phiivt .eq. 1.0)then
2514 bl = 1.
2515 al = 2.
2516 aa = PI
2517 ba = 2.
2518 qe = 1.0
2519 endif
2520 qe = min(qe, 1.0)
2522 !.. Fall speed and ventilation (Best number formulation)
2523 xn = 2./rhobar(cc)*(rhobar(cc)-rhoair(k))*G_HOME*rhoair(k)/mu**2 * &
2524 (fourthirdspi*rhobar(cc))*alphstr*al**2/(aa*qe) * qe**(3./4.)
2525 bx = deltastr(cc)+2.+2.*bl-ba
2527 !.. Number-average Best Number
2528 xm = xn*ani(cc)**bx * (gamma(NU+bx))*i_gammnu
2530 !.. The following fall speed coefficients are from
2531 !.. Mitchell and Heymsfiled 2005
2532 f_c1 = 4.0 / (5.83 * 5.83 * SQRT(0.6))
2533 f_c2 = (5.83 * 5.83) / 4.0
2534 bm = ((f_c1 * SQRT(xm)) / &
2535 (2.0 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0) * &
2536 (SQRT(1.0 + f_c1 * SQRT(xm))))) - &
2537 (1.0e-5 * xm) / &
2538 (f_c2 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0)**2.0)
2539 am = ((f_c2 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0)**2.0) - &
2540 (1.0e-5 * xm)) / (xm**bm)
2542 if(xm.gt.1.e8)then
2543 am = 1.0865
2544 bm = 0.499
2545 endif
2547 !.. Reynolds Number
2548 Nre = am*xm**bm
2550 !.. Number-averaged ice fall speed
2551 vtrni1(cc,k) = min((mu/rhoair(k)*0.5 * am*(xn)**bm*ani(cc)**(bx*bm-1.) * &
2552 (gamma(NU+bx*bm-1.))*i_gammnu),25.)
2554 !.. Mass-averaged ice fall speed
2555 vtrmi1(cc,k) = min((mu/rhoair(k)*0.5 * am*(xn)**bm*ani(cc)**(bx*bm-1.) * &
2556 (gamma(NU+bx*bm-1.+2.+deltastr(cc)))/(gamma(NU+2.+deltastr(cc)))),25.)
2558 !.. Remove size sorting when melting
2559 if(temp.gt.T0.and.qmlt(cc).lt.0.) then
2560 vtrni1(cc,k) = vtrmi1(cc,k)
2561 endif
2563 !--------------------------------------------------------------------------------------------------------------!
2564 !--------------------------------------------------------------------------------------------------------------!
2565 !.. Determine the ice radiative properties
2566 !.. For simlicity this is
2567 !.. Planar ice: integral of a^2c / integral of a^2
2568 !.. Columnar ice: integral of a^2c / integral of a*c
2570 gamma_arg = NU-1.+deltastr(cc)
2571 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2572 phirad = cni(cc)/ani(cc)*gamma_tab(gi)*i_gammnu
2573 if(phirad.le.1.) then
2574 gamma_arg = NU+2.+deltastr(cc)
2575 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2576 gamma_arg = NU+2.
2577 gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2578 effi(cc,k) = cni(cc)*gamma_tab(gi)/gamma_tab(gi2)
2579 else
2580 gamma_arg = NU+2.+deltastr(cc)
2581 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2582 gamma_arg = NU+1.+deltastr(cc)
2583 gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2584 effi(cc,k) = ani(cc)*gamma_tab(gi)/gamma_tab(gi2)
2585 endif
2587 !.. For advection use volume and volume times aspect ratio mixing ratios
2588 ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k)
2589 ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k)
2591 ai(cc,k)=max(ai(cc,k),QASMALL)
2592 ci(cc,k)=max(ci(cc,k),QASMALL)
2593 ni(cc,k)=max(ni(cc,k),QNSMALL)
2595 else
2596 !.. If not enough ice convert to vapor
2597 qv(k)=qv(k)+qi(cc,k)
2598 theta(k)=theta(k)+theta(k)*i_temp*(qi(cc,k)*xxls)*i_cp
2599 qi(cc,k)=0.
2600 ai(cc,k)=0.
2601 ni(cc,k)=0.
2602 ci(cc,k)=0.
2603 vtrzi1(cc,k)=0.
2604 vtrmi1(cc,k)=0.
2605 vtrni1(cc,k)=0.
2606 endif !.. If ice exists
2608 !--------------------------------------------------------------------------------------------------------------!
2609 !--------------------------------------------------------------------------------------------------------------!
2610 !.. Output diagnostics
2611 gamma_arg = NU-1.+deltastr(cc)
2612 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2614 !.. For aspect ratio diagnostic set lower size limit
2615 if(ai(cc,k).gt.1.e-12.and.ci(cc,k).gt.1.e-12.and. &
2616 ani(cc).gt.2.e-6.and.cni(cc).gt.2.e-6.and.qi(cc,k).gt.1.e-9) then
2617 phii1d(cc,k) = cni(cc)/ani(cc)*gamma_tab(gi)*i_gammnu
2618 rhopo1d(cc,k) = rhobar(cc)
2619 endif
2621 gamma_arg = NU+3.+deltastr(cc)
2622 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2623 gamma_arg = NU+2.+deltastr(cc)
2624 gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2625 gamma_arg = NU+2.+(2.*deltastr(cc))
2626 gi3=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2628 if(ai(cc,k).gt.1.e-12.and.ci(cc,k).gt.1.e-12.and. &
2629 ani(cc).gt.2.e-6.and.cni(cc).gt.2.e-6.and.qi(cc,k).gt.1.e-9) then
2630 di1d(cc,k) = MAX((2.*ani(cc)*gamma_tab(gi)/gamma_tab(gi2)), &
2631 (2.*cni(cc)*gamma_tab(gi3)/gamma_tab(gi2)))
2632 vmi1d(cc,k) = vtrmi1(cc,k)
2633 endif
2635 if(ni(cc,k).gt.0..and.ani(cc).gt.0.) then
2636 vtrzi1(cc,k) = 0.176/0.93*(6./PI)**2*((fourthirdspi*rhobar(cc)* &
2637 (ao)**(1.-deltastr(cc)))/(2.**(2.+deltastr(cc))))**2 * &
2638 rhoair(k)/900./900.*ni(cc,k)*(2.*ani(cc))**(2.*(2.+deltastr(cc))) * &
2639 gamma(NU+2.*(2.+deltastr(cc)))*i_gammnu
2640 vtrzi1(cc,k) = 1.e18*vtrzi1(cc,k)
2641 else
2642 vtrzi1(cc,k) = 0.
2643 endif
2644 enddo !.. Loop over all ice species
2646 dbzsum = dbzr + vtrzi1(ICE1,k) + vtrzi1(ICE2,k) + vtrzi1(ICE3,k)
2647 if(dbzsum.ne.0.) then
2648 dbz1d(k) = 10.*Log10(dbzsum)
2649 dbz1d(k) = MAX(-35.,dbz1d(k))
2650 endif
2652 !--------------------------------------------------------------------------------------------------------------!
2653 !--------------------------------------------------------------------------------------------------------------!
2654 !.. Ice effective radius is a mass-weighting of all ice species
2655 totaliwc = 0.
2656 do cc = 1, cat
2657 totaliwc = totaliwc + qi(cc,k)
2658 enddo
2659 if(totaliwc.gt.0.) then
2660 radw1 = MAX((MIN((qi(ICE1,k)/totaliwc),1.)),0.)
2661 radw2 = MAX((MIN((qi(ICE2,k)/totaliwc),1.)),0.)
2662 radw3 = MAX((MIN((qi(ICE3,k)/totaliwc),1.)),0.)
2663 effi1d(k) = MAX((MIN((effi(ICE1,k)*radw1 + effi(ICE2,k)*radw2 + &
2664 effi(ICE3,k)*radw3),999.e-6)),5.e-6)
2665 else
2666 effi1d(k)=0.
2667 endif
2669 !--------------------------------------------------------------------------------------------------------------!
2670 !--------------------------------------------------------------------------------------------------------------!
2671 !.. Check for rain or ice to do sedimentation
2672 if(qr(k).gt.QSMALL) then
2673 sedi=.true.
2674 endif
2675 do cc = 1, cat
2676 if(qi(cc,k).gt.QSMALL) then
2677 sedi=.true.
2678 endif
2679 enddo
2681 !--------------------------------------------------------------------------------------------------------------!
2682 !--------------------------------------------------------------------------------------------------------------!
2683 !.. else on domicro
2684 else
2686 !.. If we don't do all of the above
2687 !.. Push to vapor and get density for seimentation in the column
2689 xxls = 3.15e6-2370.*temp+0.3337e6
2690 xxlv = 3.1484E6-2370.*temp
2692 i_temp= 1./temp
2693 rhoair(k)= pres_e(k)/(RD*temp)
2694 i_rhoair(k) = 1./rhoair(k)
2696 qv(k)=qv(k)+qc(k)+qr(k)
2697 theta(k)=theta(k)+theta(k)*i_temp*((qc(k)+qr(k))*xxlv)*i_cp
2698 do cc = 1, cat
2699 qv(k)=qv(k)+qi(cc,k)
2700 theta(k)=theta(k)+theta(k)*i_temp*((qi(cc,k))*xxls)*i_cp
2701 qi(cc,k)=0.
2702 ni(cc,k)=0.
2703 ai(cc,k)=0.
2704 ci(cc,k)=0.
2705 vtrmi1(cc,k)=0.
2706 vtrni1(cc,k)=0.
2707 vtrzi1(cc,k)=0.
2708 vmi1d(cc,k)=0.
2709 di1d(cc,k)=0.
2710 phii1d(cc,k)=0.
2711 rhopo1d(cc,k)=0.
2712 qipre1d(cc) =0.
2713 icetype1d(cc,k)=0.
2714 enddo
2716 effc1d(k)=0.
2717 effi1d(k)=0.
2718 dbz1d(k)=-35.
2719 qc(k)=0.
2720 qr(k)=0.
2721 nr(k)=0.
2722 vtrm(k)=0.
2723 vtrn(k)=0.
2724 qrpre1d=0.
2726 endif !.. If domicro
2727 enddo !.. do k = kts, kte
2729 !--------------------------------------------------------------------------------------------------------------!
2730 !--------------------------------------------------------------------------------------------------------------!
2731 !.. Sedimentation (Reisner et al. 1998)
2732 !.. Fallout terms are calculated on split time steps to ensure numerical stability (Courant < 1)
2734 nstep = 1
2735 maxfall=0.
2737 if(sedi) then
2738 do k = kte,kts,-1
2740 !.. Calculate number of split timesteps
2741 maxfall = MAX(vtrmi1(ICE1,k),vtrni1(ICE1,k),vtrmi1(ICE2,k),vtrni1(ICE2,k), &
2742 vtrmi1(ICE3,k),vtrni1(ICE3,k),vtrm(k),vtrn(k))
2743 nstep = MAX(INT(maxfall*dt/dzmic(k)+1.),nstep)
2745 do cc = 1, cat
2746 qi(cc,k) = qi(cc,k)*rhoair(k)
2747 ni(cc,k) = ni(cc,k)*rhoair(k)
2748 ai(cc,k) = ai(cc,k)*rhoair(k)
2749 ci(cc,k) = ci(cc,k)*rhoair(k)
2750 enddo
2752 qr(k) = qr(k)*rhoair(k)
2753 nr(k) = nr(k)*rhoair(k)
2754 enddo
2756 do nn = 1, nstep
2757 do k = kts, kte
2759 do cc = 1, cat
2760 fluxqi(cc,k) = vtrmi1(cc,k)*qi(cc,k)
2761 fluxni(cc,k) = vtrni1(cc,k)*ni(cc,k)
2762 fluxai(cc,k) = vtrmi1(cc,k)*ai(cc,k)
2763 fluxci(cc,k) = vtrmi1(cc,k)*ci(cc,k)
2765 enddo
2766 fluxqr(k) = vtrm(k)*qr(k)
2767 fluxnr(k) = vtrn(k)*nr(k)
2768 enddo
2770 do cc = 1, cat
2771 falltndqi(cc) = fluxqi(cc,KTE)/dzmic(KTE)
2772 falltndni(cc) = fluxni(cc,KTE)/dzmic(KTE)
2773 falltndai(cc) = fluxai(cc,KTE)/dzmic(KTE)
2774 falltndci(cc) = fluxci(cc,KTE)/dzmic(KTE)
2775 enddo
2776 falltndqr = fluxqr(KTE)/dzmic(KTE)
2777 falltndnr = fluxnr(KTE)/dzmic(KTE)
2779 do cc = 1, cat
2780 qi(cc,KTE) = qi(cc,KTE)-falltndqi(cc)*dt/nstep
2781 ni(cc,KTE) = ni(cc,KTE)-falltndni(cc)*dt/nstep
2782 ai(cc,KTE) = ai(cc,KTE)-falltndai(cc)*dt/nstep
2783 ci(cc,KTE) = ci(cc,KTE)-falltndci(cc)*dt/nstep
2784 enddo
2785 qr(KTE) = qr(KTE)-falltndqr*dt/nstep
2786 nr(KTE) = nr(KTE)-falltndnr*dt/nstep
2789 do k = kte-1,kts,-1
2790 do cc = 1, cat
2791 falltndqi(cc) = (fluxqi(cc,k+1)-fluxqi(cc,k))/dzmic(k)
2792 falltndni(cc) = (fluxni(cc,k+1)-fluxni(cc,k))/dzmic(k)
2793 falltndai(cc) = (fluxai(cc,k+1)-fluxai(cc,k))/dzmic(k)
2794 falltndci(cc) = (fluxci(cc,k+1)-fluxci(cc,k))/dzmic(k)
2796 enddo
2797 falltndqr = (fluxqr(k+1)-fluxqr(k))/dzmic(k)
2798 falltndnr = (fluxnr(k+1)-fluxnr(k))/dzmic(k)
2800 do cc = 1, cat
2801 qi(cc,k) = qi(cc,k)+falltndqi(cc)*dt/nstep
2802 ni(cc,k) = ni(cc,k)+falltndni(cc)*dt/nstep
2803 ai(cc,k) = ai(cc,k)+falltndai(cc)*dt/nstep
2804 ci(cc,k) = ci(cc,k)+falltndci(cc)*dt/nstep
2806 !.. GET PRECIPITATION AND SNOWFALL ACCUMULATION DURING THE TIME STEP
2807 !.. FACTOR OF 1000 CONVERTS FROM M TO MM, BUT DIVISION BY DENSITY
2808 !.. OF LIQUID WATER CANCELS THIS FACTOR OF 1000
2810 !.. Ice precipitation rate
2811 if(k.eq.kts.and.qi(cc,k).gt.1.e-9) then
2812 qipre1d(cc) = qipre1d(cc)+fluxqi(cc,k)*dt/nstep
2813 endif
2815 enddo
2816 qr(k) = qr(k)+falltndqr*dt/nstep
2817 nr(k) = nr(k)+falltndnr*dt/nstep
2819 !.. Rain precipitation rate
2820 if(k.eq.kts.and.qr(k).gt.1.e-9) then
2821 qrpre1d = qrpre1d + fluxqr(k)*dt/nstep
2822 endif
2823 enddo
2824 enddo
2826 do k = kts, kte
2827 do cc = 1, cat
2828 qi(cc,k) = qi(cc,k)*i_rhoair(k)
2829 ni(cc,k) = ni(cc,k)*i_rhoair(k)
2830 ai(cc,k) = ai(cc,k)*i_rhoair(k)
2831 ci(cc,k) = ci(cc,k)*i_rhoair(k)
2832 enddo
2833 qr(k) = qr(k)*i_rhoair(k)
2834 nr(k) = nr(k)*i_rhoair(k)
2835 enddo
2836 endif
2838 !--------------------------------------------------------------------------------------------------------------!
2839 !--------------------------------------------------------------------------------------------------------------!
2840 !.. Check after sedimentation
2841 do k = kts, kte
2842 temp = theta(k)/(100000./pres_e(k))**(RCP)
2843 i_temp= 1./temp
2844 xxls = 3.15e6-2370.*temp+0.3337e6
2845 xxlv = 3.1484E6-2370.*temp
2847 do cc = 1, CAT
2848 if(qi(cc,k).gt.QSMALL) then
2850 !.. Keep variables in bounds
2851 ni(cc,k)=max(ni(cc,k),QNSMALL)
2852 ai(cc,k)=max(ai(cc,k),QASMALL)
2853 ci(cc,k)=max(ci(cc,k),QASMALL)
2854 ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
2855 cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6)
2857 !.. Set incoming limit after sedimentation
2858 !.. smallest values for which we care about shape
2859 if(ai(cc,k).lt.1.e-12.or.ci(cc,k).lt.1.e-12) then
2860 ai(cc,k)=min(ai(cc,k),ci(cc,k))
2861 ci(cc,k)=ai(cc,k)
2862 ani(cc)=(ai(cc,k)/ni(cc,k))**0.3333333333
2863 cni(cc)=(ci(cc,k)/ni(cc,k))**0.3333333333
2864 endif
2866 if(ani(cc).lt.2.e-6.or.cni(cc).lt.2.e-6) then
2867 ani(cc)=2.e-6
2868 cni(cc)=2.e-6
2869 endif
2871 ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k)
2872 ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k)
2874 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
2876 !.. Final check on aggregates
2877 if(cc.eq.ICE3) then
2878 if(temp.le.T0) then
2879 rhobar(cc)=50.
2880 phiagg3 = 0.2
2881 gamma_arg = NU-1.+deltastr(cc)
2882 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2883 cni(cc)=phiagg3*ani(cc)*gammnu/gamma_tab(gi)
2884 endif
2886 ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k)
2887 ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k)
2889 if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then
2890 deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao))
2891 else
2892 deltastr(cc)=1.
2893 endif
2894 !.. Keep aggregates spherical or oblate (for now)
2895 deltastr(cc) = min(deltastr(cc),1.)
2897 alphstr=ao**(1.-deltastr(cc))
2898 gamma_arg = NU+2.+deltastr(cc)
2899 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2900 ani(cc) = ((qi(cc,k))/(ni(cc,k)*rhobar(cc)*fourthirdspi*alphstr*&
2901 gamma_tab(gi)/gammnu))**(1./(2.+deltastr(cc)))
2903 ani(cc) = max(ani(cc),2.e-6)
2904 cni(cc)=alphstr*ani(cc)**deltastr(cc)
2906 !.. Final check on aggregate size (implicit breakup)
2907 if(ani(cc).gt.0.5e-3) then
2908 ani(cc)=0.5e-3
2909 cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc)
2910 gamma_arg = NU+2.+deltastr(cc)
2911 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2912 ni(cc,k) = (qi(cc,k))/(rhobar(cc)*fourthirdspi*ao**(1.-deltastr(cc))* &
2913 ani(cc)**(2.+deltastr(cc))*gamma_tab(gi)/gammnu)
2914 ni(cc,k)=max(ni(cc,k),QNSMALL)
2915 endif
2917 ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k)
2918 ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k)
2920 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
2921 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
2922 alphstr, alphv, betam)
2924 endif !.. Check on aggregates
2926 call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), &
2927 ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), &
2928 alphstr, alphv, betam)
2930 else
2931 qv(k)=qv(k)+qi(cc,k)
2932 theta(k)=theta(k)+theta(k)*i_temp*(qi(cc,k)*xxls)*i_cp
2933 qi(cc,k)=0.
2934 ni(cc,k)=0.
2935 ai(cc,k)=0.
2936 ci(cc,k)=0.
2937 endif
2938 enddo
2940 if(qr(k).gt.QSMALL) then
2941 nr(k) = max(nr(k),QNSMALL)
2942 lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333
2943 n0rr = nr(k)*lamr
2944 if(lamr.LT.LAMMINR) then
2945 lamr = LAMMINR
2946 n0rr = lamr**4*qr(k)/(PI*RHOW)
2947 nr(k) = n0rr/lamr
2948 elseif(lamr.gt.LAMMAXR) then
2949 lamr = LAMMAXR
2950 n0rr = lamr**4*qr(k)/(PI*RHOW)
2951 nr(k) = n0rr/lamr
2952 endif
2953 else
2954 qv(k)=qv(k)+qr(k)
2955 theta(k)=theta(k)+theta(k)*i_temp*(qr(k)*xxlv)*i_cp
2956 qr(k) =0.
2957 nr(k) =0.
2958 endif
2960 enddo
2962 end subroutine me_ishmael
2964 !--------------------------------------------------------------------------------------------------------------!
2965 !--------------------------------------------------------------------------------------------------------------!
2966 !.. This subroutine keeps the ice properties in check
2968 subroutine var_check(NU, ao, fourthirdspi, gammnu, qidum, dsdum, ani, &
2969 cni, rni, rbdum, nidum, aidum, cidum, alphstr, alphv, betam)
2971 implicit none
2973 REAL, INTENT(IN) :: NU, ao, fourthirdspi, gammnu, qidum
2974 REAL :: voltmp, maxsize, gamma_arg
2975 INTEGER :: gi
2976 REAL, INTENT(INOUT) :: dsdum, ani, cni, rbdum, nidum, aidum, cidum
2977 REAL, INTENT(OUT) :: rni, alphstr, alphv, betam
2979 !.. Deltastr check
2980 if(dsdum.lt.0.55) then
2981 dsdum=0.55
2982 ani=(cni/(ao**(1.-dsdum)))**(1./dsdum)
2983 aidum=ani**2*cni*nidum
2984 cidum=cni**2*ani*nidum
2985 else if (dsdum.gt.1.3) then
2986 dsdum=1.3
2987 cni=ao**(1.-dsdum)*ani**dsdum
2988 aidum=ani**2*cni*nidum
2989 cidum=cni**2*ani*nidum
2990 endif
2992 alphstr=ao**(1.-dsdum)
2993 alphv = fourthirdspi*alphstr
2994 betam = 2.0 + dsdum
2995 gamma_arg = NU+2.+dsdum
2996 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
2998 !.. Ice density check
2999 !.. Keep ice density between 50 and RHOI=920 kg m^-3
3000 if(ani.gt.2.e-6) then
3001 rbdum = qidum*gammnu/(nidum*alphv* &
3002 ani**betam*gamma_tab(gi))
3003 else
3004 rbdum = RHOI
3005 endif
3007 if(rbdum.gt.RHOI) then
3008 rbdum=RHOI
3009 ani=((qidum*gammnu)/(rbdum*nidum*alphv* &
3010 gamma_tab(gi)))**(1./betam)
3011 cni=ao**(1.-dsdum)*ani**dsdum
3012 aidum=ani**2*cni*nidum
3013 cidum=cni**2*ani*nidum
3014 elseif(rbdum.lt.50.) then
3015 rbdum=50.
3016 ani=((qidum*gammnu)/(rbdum*nidum*alphv* &
3017 gamma_tab(gi)))**(1./betam)
3018 cni=ao**(1.-dsdum)*ani**dsdum
3019 aidum=ani**2*cni*nidum
3020 cidum=cni**2*ani*nidum
3021 endif
3023 !.. Small ice limit
3024 !.. Keep rni > 2 microns
3025 rni= (qidum*3./(nidum*rbdum*4.*PI* &
3026 (gamma_tab(gi)/gammnu)))**0.333333333333
3027 if(rni.lt.2.e-6) then
3028 rni=2.e-6
3029 nidum=3.*qidum*gammnu/(4.*pi*rbdum*rni**3* &
3030 (gamma_tab(gi)))
3031 ani=((qidum*gammnu)/(rbdum*nidum*alphv* &
3032 gamma_tab(gi)))**(1./betam)
3033 cni=ao**(1.-dsdum)*ani**dsdum
3034 aidum=ani**2*cni*nidum
3035 cidum=cni**2*ani*nidum
3036 endif
3038 !.. Large ice limit
3039 !.. This is a number weighted diameter of 8 mm and
3040 !.. a spherical mass weighted diameter of 14 mm
3041 maxsize = max(ani,cni)
3042 if (maxsize.gt.1.e-3) then
3043 if(ani.ge.cni) then
3044 ani = 1.e-3
3045 nidum=qidum*gammnu/(fourthirdspi*rbdum*ao**(1.-dsdum)*ani**(2.+dsdum)* &
3046 (gamma_tab(gi)))
3047 cni=ao**(1.-dsdum)*ani**dsdum
3048 aidum=ani**2*cni*nidum
3049 cidum=cni**2*ani*nidum
3050 else
3051 cni = 1.e-3
3052 ani = (cni/(ao**(1.-dsdum)))**(1./dsdum)
3053 nidum=qidum*gammnu/(fourthirdspi*rbdum*ao**(1.-dsdum)*ani**(2.+dsdum)* &
3054 (gamma_tab(gi)))
3055 aidum=ani**2*cni*nidum
3056 cidum=cni**2*ani*nidum
3057 endif
3058 rni= (qidum/(nidum*rbdum*fourthirdspi* &
3059 (gamma_tab(gi)/gammnu)))**0.333333333333
3060 end if
3062 return
3063 end subroutine var_check
3065 !--------------------------------------------------------------------------------------------------------------!
3066 !--------------------------------------------------------------------------------------------------------------!
3067 !.. This is the ventilation, fall speed, vapor growth routine
3068 !.. See Harrington et al. 2013
3070 subroutine vaporgrow(dt, ani, cni, rni, igr, nidum, temp, rimesum, presdum, &
3071 NU, alphstr, sui, sup, qvs, qvi, mu, iwci, rhodum, qidum, dv, kt, ao, &
3072 nsch, npr, gammnu, i_gammnu, fourthirdspi, svpi, xxls, xxlv, xxlf, &
3073 capgam, vtbarb, vtbarbm, vtbarbz, anf, cnf, rnf, iwcf, fvdum, fhdum, &
3074 rbdum, dsdum, rdout, dsdumout)
3077 implicit none
3079 real, parameter :: QASMALL= 1.e-19 ! Smallest ice size for sublimation (squared)
3080 REAL, INTENT(IN) :: dt, ani, cni, rni, igr, nidum, temp, presdum
3081 REAL, INTENT(IN) :: NU, alphstr, mu, iwci, rhodum, qidum
3082 REAL, INTENT(IN) :: dv, kt, ao, nsch, npr, gammnu, i_gammnu, fourthirdspi, svpi
3083 REAL, INTENT(IN) :: xxls, xxlv, xxlf, capgam, rimesum, sui, sup, qvs, qvi, dsdum
3084 REAL :: gammnubet, phii, fs, alphanr, bl, al, aa, ba
3085 REAL :: xn, bx, xm, am, bm, f_c1, f_c2, del1, del2, alpha, afn, del
3086 REAL :: xvent, ntherm, bv1, bv2, gv, bt1, bt2, gt, vtbranch, rhodep
3087 REAL :: videp, vmin, betavol, vi, vf, cf, maxsui, qe, Nre, phif, gamma_arg, rbdumtmp
3088 INTEGER :: gi, gi2
3089 REAL, INTENT(INOUT) :: vtbarb, vtbarbm, vtbarbz
3090 REAL, INTENT(OUT) :: anf, cnf, rnf, iwcf, rdout
3091 REAL, INTENT(INOUT) :: fvdum, fhdum, dsdumout, rbdum
3093 fvdum = 1.0
3094 fhdum = 1.0
3096 dsdumout = dsdum
3097 gamma_arg = NU-1.+dsdum
3098 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
3099 phii = cni/ani*gamma_tab(gi)*i_gammnu
3100 fs = capgam/rni
3101 alphanr = ani/rni**(3./(2.+igr))
3103 !.. Determine coefficeients for the Best number for fall speeds
3104 if(phii.lt.1.0)then
3105 bl = 1.
3106 al = 2.
3107 aa = PI
3108 ba = 2.
3109 qe = (1.-phii)*(rbdum/RHOI) + phii
3110 else if(phii .gt. 1.0)then
3111 al = 2.0
3112 bl = 1.
3113 aa = PI*alphstr
3114 ba = dsdum + 1.
3115 qe = 1.0
3116 else if(phii .eq. 1.0)then
3117 bl = 1.
3118 al = 2.
3119 aa = pi
3120 ba = 2.
3121 qe = 1.0
3122 endif
3123 qe = min(qe, 1.0)
3125 !.. Fall speed and ventilation (Best number formulation)
3126 xn = 2./rbdum*(rbdum-rhodum)*G_HOME*rhodum/mu**2 * &
3127 (fourthirdspi*rbdum)*alphstr*al**2/(aa*qe) * qe**(3./4.)
3128 bx = dsdum+2.+2.*bl-ba
3130 !.. Number-average Best Number
3131 xm = xn*ani**bx * (gamma(nu+bx))*i_gammnu
3133 !.. The following fall speed coefficients are from
3134 !.. Mitchell and Heymsfiled 2005
3135 f_c1 = 4.0 / (5.83 * 5.83 * SQRT(0.6))
3136 f_c2 = (5.83 * 5.83) / 4.0
3137 bm = ((f_c1 * SQRT(xm)) / &
3138 (2.0 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0) * &
3139 (SQRT(1.0 + f_c1 * SQRT(xm))))) - &
3140 (1.0e-5 * xm) / &
3141 (f_c2 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0)**2.0)
3142 am = ((f_c2 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0)**2.0) - &
3143 (1.0e-5 * xm)) / (xm**bm)
3145 if(xm.gt.1.e8)then
3146 am = 1.0865
3147 bm = 0.499
3148 endif
3150 !.. Reynolds Number
3151 Nre = am*xm**bm
3153 !.. Number-averaged ice fall speed
3154 vtbarb = mu/rhodum*0.5 * am*(xn)**bm*ani**(bx*bm-1.) * &
3155 (gamma(nu+bx*bm-1.))*i_gammnu
3157 !.. Mass-averaged ice fall speed
3158 vtbarbm = mu/rhodum*0.5 * am*(xn)**bm*ani**(bx*bm-1.) * &
3159 (gamma(nu+bx*bm-1.+2.+dsdum))/(gamma(nu+2.+dsdum))
3161 !.. Reflectivity-weighted fall speed (if needed)
3162 vtbarbz = mu/rhodum*0.5 * am*(xn)**bm*ani**(bx*bm-1.) * &
3163 exp(gammln(nu+bx*bm-1.+4.+2.*dsdum))/exp(gammln(nu+4.+2.*dsdum))
3165 !.. Calculate Ventilation
3166 xvent = nsch**0.333333333*Nre**0.5
3167 ntherm = Nre**0.5*npr**0.333333333
3168 if(xvent.le.1.0)then
3169 bv1 = 1.0
3170 bv2 = 0.14
3171 gv = 2.
3172 else
3173 bv1 = 0.86
3174 bv2 = 0.28
3175 gv = 1.
3176 endif
3177 if(ntherm.lt.1.4)then
3178 bt1 = 1.0
3179 bt2 = 0.108
3180 gt = 2.0
3181 else
3182 bt1 = 0.78
3183 bt2 = 0.308
3184 gt = 1.0
3185 endif
3186 fvdum = bv1 + bv2*xvent**gv
3187 fhdum = bt1 + bt2*ntherm**gt
3189 if(temp.le.T0) then !.. If T < T0 do vapor growth/sublimation (otherwise assume water on surface)
3191 !.. Fall speed needed to determine when branching occurs
3192 vtbranch = vtbarbm
3194 if(sup.ge.0.) then
3195 maxsui = 1.
3196 elseif(sui.ge.0..and.qvi.lt.qvs) then
3197 maxsui = ((sui+1.)*qvi)/(qvs-qvi) - qvi/(qvs-qvi)
3198 maxsui = min(maxsui,1.)
3199 maxsui = max(maxsui,0.)
3200 else
3201 maxsui = 0.
3202 endif
3204 !.. Vapor growth density
3205 if(igr.le.1.) then !.. Planar
3206 if(vtbranch.gt.0.) then
3207 if(ani.gt.SQRT((dv*PI*2.*cni)/(vtbranch*nu))) then
3208 rhodep = (RHOI*igr)*maxsui + RHOI*(1.-maxsui)
3209 else
3210 rhodep = RHOI
3211 endif
3212 else
3213 rhodep = RHOI
3214 endif
3215 else !.. Columnar
3216 rhodep = (RHOI/igr)*maxsui + RHOI*(1.-maxsui)
3217 endif
3218 !.. high limit on rhodep for vapor growth = 700 kg m^-3
3219 !.. Cotton et al. 2012 (QJRMS)
3220 rhodep = min(rhodep,700.)
3222 !.. Vapor growth solution which includes heating from rime mass
3223 !.. See J.P. Chen's thesis
3224 alpha = (dv*fvdum*svpi*xxls)/(RV*kt*fhdum*temp)
3225 if(nidum.gt.0.0.and.rimesum.gt.0.0) then
3226 del1 = (xxlf*(rimesum/nidum)/(4.*PI*kt*fhdum*capgam)) * &
3227 ((temp + alpha*((xxls/(RV*temp))-1.))**(-1.0))
3228 else
3229 del1 = 0.0
3230 endif
3231 del2 = sui * &
3232 (((temp/alpha) + ((xxls/(RV*temp))-1.))**(-1.0))
3233 del = del1 + del2
3234 afn = ((dv*fvdum*polysvp(temp,1))/(RV*temp)) * &
3235 (sui - del*((xxls/(RV*temp))-1.))
3237 !.. During sublimation using polynomial removal of density
3238 if(afn.lt.0.0) then
3239 rhodep = rbdum
3240 videp = rni**3
3241 vmin = (10.e-6)**3
3242 if(Vmin.lt.videp)then
3243 betavol = log(RHOI/rbdum)*1./(log(Vmin/videp))
3244 rhodep = rbdum*(1.+betavol)
3245 else
3246 rhodep = rbdum
3247 endif
3248 endif
3249 rhodep=max(rhodep,50.)
3250 rhodep=min(rhodep,RHOI)
3252 gamma_arg = NU+2.+dsdum
3253 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
3254 gammnubet = gamma_tab(gi)
3256 !.. characteristic r-axis and a-axis after growth timestep
3257 rnf = (max((rni**2 + 2.*afn*fs/rhodep*gammnu/gammnubet*dt),QASMALL))**(0.5)
3258 anf = alphanr*rnf**(3./(2.+igr))
3260 !.. Do not sublimation change the shape of ice from
3261 !.. prolate to oblate or vice versa
3262 phif = phii*(rnf**3/rni**3)**((igr-1.)/(igr+2.))
3263 if(sui.lt.0.0.or.afn.lt.0.0) then
3264 if(phii.gt.1.0.and.phif.lt.1.0) then
3265 phif = phii
3266 alphanr = ani/rni
3267 anf = alphanr*rnf
3268 endif
3269 if(phii.lt.1.0.and.phif.gt.1.0) then
3270 phif = phii
3271 alphanr = ani/rni
3272 anf = alphanr*rnf
3273 endif
3275 !.. Do not let sublimation create extreme shapes
3276 if(phii.gt.1.) then
3277 if(phif.gt.phii) then
3278 phif = phii
3279 alphanr = ani/rni
3280 anf = alphanr*rnf
3281 endif
3282 else
3283 if(phif.lt.phii) then
3284 phif = phii
3285 alphanr = ani/rni
3286 anf = alphanr*rnf
3287 endif
3288 endif
3289 endif
3291 vi = fourthirdspi*rni**3*gamma_tab(gi)*i_gammnu
3292 vf = fourthirdspi*rnf**3*gamma_tab(gi)*i_gammnu
3293 rdout = rhodep
3294 rbdumtmp = rbdum*(vi/vf) + rhodep*(1.-vi/vf)
3295 rbdumtmp = min(rbdumtmp,RHOI)
3296 iwcf = nidum*rbdumtmp*vf
3298 !.. Update delta* from vapor growth/sublimation
3299 if(igr.ne.1.0)then
3301 if(anf.gt.(1.1*ao)) then
3302 dsdumout = (3.*log(rnf)-2.*log(anf)-log(ao))/ &
3303 (log(anf)-log(ao))
3304 else
3305 dsdumout=1.
3306 endif
3307 endif
3309 !.. Do not let particles sublimate to sizes that are too small
3310 if(afn.lt.0..and.rnf.lt.1.e-6) then
3311 rbdum =RHOI
3312 rdout =RHOI
3313 dsdumout =1.
3314 phif =1.
3315 alphanr = ani/rni
3316 anf = alphanr*rnf
3317 endif
3319 if(afn.lt.0..and.anf.lt.1.e-6) then
3320 rbdum =RHOI
3321 rdout =RHOI
3322 dsdumout =1.
3323 phif =1.
3324 alphanr = ani/rni
3325 anf = alphanr*rnf
3326 endif
3329 !.. Sublimation check
3330 if(afn.lt.0..and.dsdumout.le.0.) then
3331 dsdumout=1.
3332 anf =rnf
3333 endif
3335 !.. C-axis after vapor growth
3336 gamma_arg = NU-1.+dsdumout
3337 gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001)
3338 cnf = phif*anf*gammnu/gamma_tab(gi)
3340 else !.. T > T0
3342 anf =ani
3343 cnf =cni
3344 rnf =rni
3345 iwcf=iwci
3346 endif
3348 end subroutine vaporgrow
3350 !--------------------------------------------------------------------------------------------------------------!
3351 !--------------------------------------------------------------------------------------------------------------!
3352 !.. This routine is used to calculate the distribution-averaged
3353 !.. capacitance (see Harrington et al. 2013)
3355 real function capacitance_gamma(ani,dsdum,NU,alphstr,i_gammnu)
3357 implicit none
3359 REAL, INTENT(IN) :: ani, dsdum, NU, alphstr, i_gammnu
3360 REAL :: a1, a2, b1, b2, c1, c2, d1, d2, gammad1, gammad2
3362 !.. Oblate Spheroid
3363 if (dsdum.le.1.0) then
3364 a1 = 0.6369427
3365 a2 = 0.57*a1
3366 b1 = 0.0
3367 b2 = 0.95
3368 c1 = a1*alphstr**b1
3369 c2 = a2*alphstr**b2
3370 d1 = b1*(dsdum - 1.0) + 1.0
3371 d2 = b2*(dsdum - 1.0) + 1.0
3372 !.. Prolate Spheroid
3373 elseif (dsdum.gt.1.) then
3374 a1 = 0.5714285
3375 a2 = 0.75*a1
3376 b1 = -1.0
3377 b2 = -0.18
3378 c1 = a1*alphstr**(b1+1.0)
3379 c2 = a2*alphstr**(b2+1.0)
3380 d1 = b1*(dsdum - 1.0) + dsdum
3381 d2 = b2*(dsdum - 1.0) + dsdum
3382 endif
3384 if(dsdum.eq.1..and.nu.eq.1.) then
3385 gammad1 = (gamma(nu+1.0))
3386 capacitance_gamma = ani*gammad1*i_gammnu
3387 elseif(dsdum.le.1.)then
3388 gammad1 = (gamma(nu+d1))
3389 gammad2 = (gamma(nu+d2))
3390 capacitance_gamma = c1*ani**d1 * gammad1*i_gammnu + &
3391 c2*ani**d2 * gammad2*i_gammnu
3392 elseif(dsdum.gt.1.) then
3393 gammad1 = (gamma(nu+d1))
3394 gammad2 = (gamma(nu+d2))
3395 capacitance_gamma = c1*ani**d1 * gammad1*i_gammnu + &
3396 c2*ani**d2 * gammad2*i_gammnu
3397 endif
3399 return
3400 end function capacitance_gamma
3402 !--------------------------------------------------------------------------------------------------------------!
3403 !--------------------------------------------------------------------------------------------------------------!
3404 SUBROUTINE wet_growth_check(NU, temp, rhodum, xxlv, xxlf, qvdum, dv, kt, qs0, &
3405 fvdum, fhdum, rimedum, rni, nidum, dgflag)
3407 IMPLICIT NONE
3409 !.. Check for wet growth conditions
3410 !.. See Lamb and Verlinde (2011)
3411 REAL, INTENT(IN) :: NU, temp, rhodum, xxlv, xxlf, qvdum, dv, kt
3412 REAL, INTENT(IN) :: fvdum, fhdum, rni, nidum
3413 REAL, INTENT(IN) :: qs0, rimedum
3414 REAL :: wetg, dum
3415 LOGICAL, INTENT(INOUT) :: dgflag
3417 dum=nidum*NU*2.*rni
3418 wetg=2.*PI*dum*(kt*fhdum*(T0-temp)+rhodum*xxlv*dv*fvdum*(qs0-qvdum))/(xxlf + (CPW*(temp-T0)))
3419 if(rimedum/rhodum.gt.wetg) then
3420 dgflag = .false.
3421 endif
3423 end subroutine wet_growth_check
3425 !--------------------------------------------------------------------------------------------------------------!
3426 !--------------------------------------------------------------------------------------------------------------!
3427 real function polysvp(T,TYPE)
3429 !--------------------------------------------------------------
3430 !.. Taken from 'module_mp_morr_two_moment.F' (WRFV3.4)
3432 !.. COMPUTE SATURATION VAPOR PRESSURE
3434 !.. POLYSVP RETURNED IN UNITS OF PA.
3435 !.. T IS INPUT IN UNITS OF K.
3436 !.. TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1)
3438 !.. REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM FLATAU ET AL. 1992,
3439 !.. TABLE 4 (RIGHT-HAND COLUMN)
3440 !--------------------------------------------------------------
3442 IMPLICIT NONE
3444 REAL DUM
3445 REAL T
3446 INTEGER TYPE
3447 !.. ice
3448 real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i
3449 data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /&
3450 6.11147274, 0.503160820, 0.188439774e-1, &
3451 0.420895665e-3, 0.615021634e-5,0.602588177e-7, &
3452 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/
3454 !.. liquid
3455 real a0,a1,a2,a3,a4,a5,a6,a7,a8
3457 !.. V1.7
3458 data a0,a1,a2,a3,a4,a5,a6,a7,a8 /&
3459 6.11239921, 0.443987641, 0.142986287e-1, &
3460 0.264847430e-3, 0.302950461e-5, 0.206739458e-7, &
3461 0.640689451e-10,-0.952447341e-13,-0.976195544e-15/
3462 real dt
3464 !.. ICE
3465 IF (TYPE.EQ.1) THEN
3467 dt = max(-80.,t-273.16)
3468 polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt)))))))
3469 polysvp = polysvp*100.
3471 END IF
3473 !.. LIQUID
3474 IF (TYPE.EQ.0) THEN
3476 dt = max(-80.,t-273.16)
3477 polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt)))))))
3478 polysvp = polysvp*100.
3480 END IF
3482 end function polysvp
3484 !--------------------------------------------------------------------------------------------------------------!
3485 !--------------------------------------------------------------------------------------------------------------!
3486 !.. Access and interpolate lookup table values (from Hugh Morrison)
3487 subroutine access_lookup_table(itabdum,dumjj,dumii,dumi,dumk,index, &
3488 dum1,dum2,dum4,dum5,proc)
3490 implicit none
3492 real :: itabdum(:,:,:,:,:), dum1,dum2,dum4,dum5,proc
3493 real :: dproc1,dproc2,iproc1,gproc1,tmp1,tmp2
3494 integer :: dumjj,dumii,dumi,dumk,index
3496 !.. get value at current density index
3497 !.. first interpolate for current rimed fraction index
3498 dproc1 =itabdum(dumjj,dumii,dumi,dumk,index)+(dum1-real(dumi))* &
3499 (itabdum(dumjj,dumii,dumi+1,dumk,index)-itabdum(dumjj,dumii,dumi,dumk,index))
3501 dproc2 =itabdum(dumjj,dumii,dumi,dumk+1,index)+(dum1-real(dumi))* &
3502 (itabdum(dumjj,dumii,dumi+1,dumk+1,index)-itabdum(dumjj,dumii,dumi,dumk+1,index))
3504 iproc1=dproc1+(dum2-real(dumk))*(dproc2-dproc1)
3506 !.. linearly interpolate to get process rates for rimed fraction index + 1
3507 dproc1 =itabdum(dumjj,dumii+1,dumi,dumk,index)+(dum1-real(dumi))* &
3508 (itabdum(dumjj,dumii+1,dumi+1,dumk,index)- &
3509 itabdum(dumjj,dumii+1,dumi,dumk,index))
3511 dproc2 =itabdum(dumjj,dumii+1,dumi,dumk+1,index)+(dum1-real(dumi))* &
3512 (itabdum(dumjj,dumii+1,dumi+1,dumk+1,index)- &
3513 itabdum(dumjj,dumii+1,dumi,dumk+1,index))
3515 gproc1=dproc1+(dum2-real(dumk))*(dproc2-dproc1)
3517 tmp1=iproc1+(dum4-real(dumii))*(gproc1-iproc1)
3519 !.. get value at density index + 1
3520 !.. first interpolate for current rimed fraction index
3521 dproc1 =itabdum(dumjj+1,dumii,dumi,dumk,index)+(dum1-real(dumi))* &
3522 (itabdum(dumjj+1,dumii,dumi+1,dumk,index)- &
3523 itabdum(dumjj+1,dumii,dumi,dumk,index))
3525 dproc2 =itabdum(dumjj+1,dumii,dumi,dumk+1,index)+ &
3526 (dum1-real(dumi))*(itabdum(dumjj+1,dumii,dumi+1,dumk+1,index)- &
3527 itabdum(dumjj+1,dumii,dumi,dumk+1,index))
3529 iproc1=dproc1+(dum2-real(dumk))*(dproc2-dproc1)
3531 !.. linearly interpolate to get process rates for rimed fraction index + 1
3532 dproc1 =itabdum(dumjj+1,dumii+1,dumi,dumk,index)+(dum1-real(dumi))* &
3533 (itabdum(dumjj+1,dumii+1,dumi+1,dumk,index)- &
3534 itabdum(dumjj+1,dumii+1,dumi,dumk,index))
3536 dproc2 =itabdum(dumjj+1,dumii+1,dumi,dumk+1,index)+(dum1-real(dumi))* &
3537 (itabdum(dumjj+1,dumii+1,dumi+1,dumk+1,index)- &
3538 itabdum(dumjj+1,dumii+1,dumi,dumk+1,index))
3540 gproc1=dproc1+(dum2-real(dumk))*(dproc2-dproc1)
3542 tmp2=iproc1+(dum4-real(dumii))*(gproc1-iproc1)
3544 !.. get final process rate
3545 proc=tmp1+(dum5-real(dumjj))*(tmp2-tmp1)
3547 end subroutine access_lookup_table
3549 !--------------------------------------------------------------------------------------------------------------!
3550 !--------------------------------------------------------------------------------------------------------------!
3551 !.. Computes ln(gamma(xx)) using Lanczos
3552 real function gammln(xx)
3554 implicit none
3556 REAL :: xx
3557 REAL*8 :: cof(6), stp, x, y, tmp, ser
3558 INTEGER :: j
3559 DATA cof, stp/76.18009172947146d0, &
3560 -86.50532032941677d0, &
3561 24.01409824083091d0, &
3562 -1.231739572450155d0, &
3563 0.1208650973866179d-2, &
3564 -0.5395239384953d-5, &
3565 2.5066282746310005d0/
3566 x = xx
3567 y = x
3568 tmp = x + 5.5d0
3569 tmp = (x + 0.5d0) * log(tmp) - tmp
3570 ser = 1.000000000190015d0
3571 do j = 1, 6
3572 y = y + 1.0d0
3573 ser = ser + cof(j) / y
3574 enddo
3576 gammln = tmp + log(stp * ser / x)
3578 return
3579 end function gammln
3581 !--------------------------------------------------------------------------------------------------------------!
3582 !--------------------------------------------------------------------------------------------------------------!
3583 REAL FUNCTION GAMMA(X)
3584 !----------------------------------------------------------------------
3585 !
3586 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
3587 ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
3588 ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
3589 ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
3590 ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
3591 ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
3592 ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
3593 ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
3594 ! MACHINE-DEPENDENT CONSTANTS.
3595 !
3596 !
3597 !*******************************************************************
3598 !*******************************************************************
3599 !
3600 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
3601 !
3602 ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
3603 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
3604 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
3605 ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
3606 ! GAMMA(XBIG) = BETA**MAXEXP
3607 ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
3608 ! APPROXIMATELY BETA**MAXEXP
3609 ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3610 ! 1.0+EPS .GT. 1.0
3611 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
3612 ! 1/XMININ IS MACHINE REPRESENTABLE
3613 !
3614 ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
3615 !
3616 ! BETA MAXEXP XBIG
3617 !
3618 ! CRAY-1 (S.P.) 2 8191 966.961
3619 ! CYBER 180/855
3620 ! UNDER NOS (S.P.) 2 1070 177.803
3621 ! IEEE (IBM/XT,
3622 ! SUN, ETC.) (S.P.) 2 128 35.040
3623 ! IEEE (IBM/XT,
3624 ! SUN, ETC.) (D.P.) 2 1024 171.624
3625 ! IBM 3033 (D.P.) 16 63 57.574
3626 ! VAX D-FORMAT (D.P.) 2 127 34.844
3627 ! VAX G-FORMAT (D.P.) 2 1023 171.489
3628 !
3629 ! XINF EPS XMININ
3630 !
3631 ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
3632 ! CYBER 180/855
3633 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
3634 ! IEEE (IBM/XT,
3635 ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
3636 ! IEEE (IBM/XT,
3637 ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
3638 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
3639 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
3640 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
3641 !
3642 !*******************************************************************
3643 !*******************************************************************
3644 !
3645 ! ERROR RETURNS
3646 !
3647 ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
3648 ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
3649 ! TO BE FREE OF UNDERFLOW AND OVERFLOW.
3650 !
3651 !
3652 ! INTRINSIC FUNCTIONS REQUIRED ARE:
3653 !
3654 ! INT, DBLE, EXP, LOG, REAL, SIN
3655 !
3656 !
3657 ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
3658 ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
3659 ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
3660 ! (ED.), SPRINGER VERLAG, BERLIN, 1976.
3661 !
3662 ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
3663 ! SONS, NEW YORK, 1968.
3664 !
3665 ! LATEST MODIFICATION: OCTOBER 12, 1989
3666 !
3667 ! AUTHORS: W. J. CODY AND L. STOLTZ
3668 ! APPLIED MATHEMATICS DIVISION
3669 ! ARGONNE NATIONAL LABORATORY
3670 ! ARGONNE, IL 60439
3671 !
3672 !----------------------------------------------------------------------
3673 implicit none
3674 INTEGER I,N
3675 LOGICAL PARITY
3676 REAL &
3677 CONV,EPS,FACT,HALF,ONE,RES,SUM,TWELVE, &
3678 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
3679 REAL, DIMENSION(7) :: C
3680 REAL, DIMENSION(8) :: P
3681 REAL, DIMENSION(8) :: Q
3682 REAL, PARAMETER :: xxx = 0.9189385332046727417803297
3683 !----------------------------------------------------------------------
3684 ! MATHEMATICAL CONSTANTS
3685 !----------------------------------------------------------------------
3686 DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/
3689 !----------------------------------------------------------------------
3690 ! MACHINE DEPENDENT PARAMETERS
3691 !----------------------------------------------------------------------
3692 DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,XINF/3.4E38/
3693 !----------------------------------------------------------------------
3694 ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
3695 ! APPROXIMATION OVER (1,2).
3696 !----------------------------------------------------------------------
3697 DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
3698 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
3699 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
3700 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
3701 DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
3702 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
3703 2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
3704 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
3705 !----------------------------------------------------------------------
3706 ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
3707 !----------------------------------------------------------------------
3708 DATA C/-1.910444077728E-03,8.4171387781295E-04, &
3709 -5.952379913043012E-04,7.93650793500350248E-04, &
3710 -2.777777777777681622553E-03,8.333333333333333331554247E-02, &
3711 5.7083835261E-03/
3712 !----------------------------------------------------------------------
3713 ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
3714 !----------------------------------------------------------------------
3715 CONV(I) = REAL(I)
3716 PARITY=.FALSE.
3717 FACT=ONE
3718 N=0
3719 Y=X
3720 IF(Y.LE.ZERO)THEN
3721 !----------------------------------------------------------------------
3722 ! ARGUMENT IS NEGATIVE
3723 !----------------------------------------------------------------------
3724 Y=-X
3725 Y1=AINT(Y)
3726 RES=Y-Y1
3727 IF(RES.NE.ZERO)THEN
3728 IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
3729 FACT=-PI/SIN(PI*RES)
3730 Y=Y+ONE
3731 ELSE
3732 RES=XINF
3733 GOTO 900
3734 ENDIF
3735 ENDIF
3736 !----------------------------------------------------------------------
3737 ! ARGUMENT IS POSITIVE
3738 !----------------------------------------------------------------------
3739 IF(Y.LT.EPS)THEN
3740 !----------------------------------------------------------------------
3741 ! ARGUMENT .LT. EPS
3742 !----------------------------------------------------------------------
3743 IF(Y.GE.XMININ)THEN
3744 RES=ONE/Y
3745 ELSE
3746 RES=XINF
3747 GOTO 900
3748 ENDIF
3749 ELSEIF(Y.LT.TWELVE)THEN
3750 Y1=Y
3751 IF(Y.LT.ONE)THEN
3752 !----------------------------------------------------------------------
3753 ! 0.0 .LT. ARGUMENT .LT. 1.0
3754 !----------------------------------------------------------------------
3755 Z=Y
3756 Y=Y+ONE
3757 ELSE
3758 !----------------------------------------------------------------------
3759 ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
3760 !----------------------------------------------------------------------
3761 N=INT(Y)-1
3762 Y=Y-CONV(N)
3763 Z=Y-ONE
3764 ENDIF
3765 !----------------------------------------------------------------------
3766 ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
3767 !----------------------------------------------------------------------
3768 XNUM=ZERO
3769 XDEN=ONE
3770 DO I=1,8
3771 XNUM=(XNUM+P(I))*Z
3772 XDEN=XDEN*Z+Q(I)
3773 END DO
3774 RES=XNUM/XDEN+ONE
3775 IF(Y1.LT.Y)THEN
3776 !----------------------------------------------------------------------
3777 ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
3778 !----------------------------------------------------------------------
3779 RES=RES/Y1
3780 ELSEIF(Y1.GT.Y)THEN
3781 !----------------------------------------------------------------------
3782 ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
3783 !----------------------------------------------------------------------
3784 DO I=1,N
3785 RES=RES*Y
3786 Y=Y+ONE
3787 END DO
3788 ENDIF
3789 ELSE
3790 !----------------------------------------------------------------------
3791 ! EVALUATE FOR ARGUMENT .GE. 12.0,
3792 !----------------------------------------------------------------------
3793 IF(Y.LE.XBIG)THEN
3794 YSQ=Y*Y
3795 SUM=C(7)
3796 DO I=1,6
3797 SUM=SUM/YSQ+C(I)
3798 END DO
3799 SUM=SUM/Y-Y+xxx
3800 SUM=SUM+(Y-HALF)*LOG(Y)
3801 RES=EXP(SUM)
3802 ELSE
3803 RES=XINF
3804 GOTO 900
3805 ENDIF
3806 ENDIF
3807 !----------------------------------------------------------------------
3808 ! FINAL ADJUSTMENTS AND RETURN
3809 !----------------------------------------------------------------------
3810 IF(PARITY)RES=-RES
3811 IF(FACT.NE.ONE)RES=FACT/RES
3812 900 GAMMA=RES
3813 RETURN
3814 ! ---------- LAST LINE OF GAMMA ----------
3815 END FUNCTION GAMMA
3817 !--------------------------------------------------------------------------------------------------------------!
3818 !--------------------------------------------------------------------------------------------------------------!
3819 !.. Get the inherent growth ratio data as a function of temperature
3820 real function get_igr(igrdatadum, temp)
3822 implicit none
3824 real :: igrdatadum(:), dum, igr1, igr2, temp
3826 !.. Inherent growth ratio (IGR) from data
3827 !.. See Lamb and Scott (1972) and Chen and Lamb (1994)
3828 if((temp-T0).ge.-59..and.(temp-T0).le.-1.) then
3829 dum = (abs(real(int(temp-T0))) + 1.) - abs(temp-T0)
3830 igr1 = igrdatadum(max((int(temp-T0)*(-1)),1))
3831 igr2 = igrdatadum(min(((int(temp-T0)*(-1))+1),60))
3832 get_igr = dum*igr1 + (1.-dum)*igr2
3833 elseif((temp-T0).gt.-1..and.(temp-T0).le.0.) then
3834 dum = 1.0 - abs(temp-T0)
3835 igr1 = 1.0
3836 igr2 = igrdatadum(1)
3837 get_igr = dum*igr1 + (1.-dum)*igr2
3838 elseif((temp-T0).ge.-60..and.(temp-T0).lt.-59.) then
3839 dum = 60.0 - abs(temp-T0)
3840 igr1 = igrdatadum(59)
3841 igr2 = igrdatadum(60)
3842 get_igr = dum*igr1 + (1.-dum)*igr2
3843 elseif((temp-T0).lt.-60.) then
3844 get_igr = igrdatadum(60)
3845 else
3846 get_igr = 1.
3847 endif
3849 end function get_igr
3851 !--------------------------------------------------------------------------------------------------------------!
3852 !--------------------------------------------------------------------------------------------------------------!
3853 !.. Aggregation routine (see Jensen et al. 2017, implentation by JYH and AAJ)
3854 !.. See Walko et al. 1995 and Meyers et al. 1997 for more details
3855 subroutine aggregation(dtdum,rhoairdum,tempdum,qdum1,ndum1,ddum1,qdum2,ndum2,ddum2, &
3856 qdum3,ndum3,ddum3,rhodum1,rhodum2,phidum1,phidum2,coltab,coltabn, &
3857 qagg1,qagg2,qagg3,nagg1,nagg2,nagg3)
3859 implicit none
3861 real, parameter :: QNSMALL= 1.25e-7 ! Smallest ice number
3862 integer icat,k,kp1,kp2,isnow,ipris,it,ihab,jhab
3863 integer, parameter :: ncat=7,npair=35,ndn=60,nftab=1000,ngam=5000,nz=1
3864 real, parameter :: rsmall=1.e-12
3865 real gnu(ncat),shape(ncat),cfmas(ncat),pwmas(ncat),&
3866 cfden(ncat),pwden(ncat)&
3867 ,cfvt(ncat),pwvt(ncat),tablo(ncat),tabhi(ncat),tabloln(ncat),&
3868 tabhiln(ncat)&
3869 ,dnmin(ncat),dnmax(ncat),vtfac(ncat),pow1(ncat),pow2(ncat),&
3870 frefac1(ncat)&
3871 ,frefac2(ncat),cfmasft(ncat),dict(ncat),dbmi(ncat),&
3872 gamm(ncat),gamn1(ncat)&
3873 ,gam(ngam,ncat),ncrossloss,rcrossloss,nselfloss,&
3874 rselfloss,ncrossgain,rcrossgain,nselfgain,rselfgain
3875 real rictmin,rictmax,rict,rictmm,dtlt,totalc(nz),totaldc(nz)
3876 real dn(nz,ncat),en(nz,ncat),qr(nz,ncat),r(nz,ncat),qq(nz,ncat)
3877 real wct1(nz,ncat),wct2(nz,ncat),rxfer(nz,ncat,ncat),&
3878 qrxfer(nz,ncat,ncat),enxfer(nz,ncat,ncat),&
3879 parxfer(nz,ncat,ncat), paqrxfer(nz,ncat,ncat),&
3880 paenxfer(nz,ncat,ncat), pprxfer(nz,ncat,ncat),&
3881 ppqrxfer(nz,ncat,ncat), ppenxfer(nz,ncat,ncat)
3882 real qs(nz),ns(nz),dns(nz),qa(nz),na(nz),dna(nz)
3883 real t(nz,ncat),tc(nz),rhoa(nz),tk(nz),denfac(nz),eff(nz)
3884 real ftable(nftab,ncat)
3885 real, intent(IN) :: coltab(ndn,ndn,npair),coltabn(ndn,ndn,npair)
3886 integer ipair(ncat,ncat),ict(nz,ncat)
3887 real, dimension(ncat,9) :: dstprms
3888 real, intent(IN) :: qdum1,ndum1,qdum2,ndum2,qdum3,ndum3,ddum1,ddum2
3889 real, intent(INOUT) :: qagg1,qagg2,qagg3,nagg1,nagg2,nagg3,ddum3
3890 real, intent(IN) :: rhodum1, rhodum2, phidum1, phidum2
3891 real :: rhoeff, phieff, rhoeffmax, phieffmax, efffact
3892 real :: source3, sink3, source4, sink4, ratioagg
3893 real, intent(IN) :: dtdum, rhoairdum, tempdum
3895 dstprms=&
3896 !-----------------------------------------------------------------------
3897 !-----------------------------------------------------------------------
3898 ! cloud rain planar columnar aggreg graup hail
3899 ! 1 2 3 4 5 6 7 ---> ncat
3900 !-----------------------------------------------------------------------
3901 reshape((/.5, .5, .318, .318, .5, .5, .5, &!shape (capacitance shape)
3902 524., 524., .333, .333, .496, 157., 471., &!cfmas (m=cfmas*D**pwmas)
3903 3., 3., 2.4, 2.4, 2.4, 3., 3., &!pwmas
3904 3173., 149., 4.836, 4.836, 3.084, 93.3, 161., &!cfvt (v = cfvt*D**pwvt)
3905 2., .5, 0.25, .25, .2, .5, .5, &!pwvt
3906 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, &!tablo (size limits for table)
3907 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, &!tabhi
3908 .1e-5, .1e-4, .1e-5, .1e-5, .1e-4, .1e-4, .1e-4, &!dnmin (characteristic diam
3909 .1e-2, .1e-1, .1e-1, .1e-1, .1e-1, .1e-1, .1e-1/), &!dnmax limits)
3910 (/ncat,9/))
3911 !-----------------------------------------------------------------------
3913 ipair=&
3914 reshape((/0,0,3,7,11,15,19 ,1,0,5,9,13,17,21 & ! collection pairs
3915 ,2,4,22,24,0,0,0 ,6,8,23,28,0,0,0 &
3916 ,10,12,25,29,35,0,0 ,14,16,26,30,32,0,0 &
3917 ,18,20,27,31,33,34,0/), &
3918 (/ncat,ncat/))
3920 do icat=1,ncat
3921 shape(icat)=dstprms(icat,1)
3922 cfmas(icat)=dstprms(icat,2)
3923 pwmas(icat)=dstprms(icat,3)
3924 cfvt (icat)=dstprms(icat,4)
3925 pwvt (icat)=dstprms(icat,5)
3926 tablo(icat)=dstprms(icat,6)
3927 tabhi(icat)=dstprms(icat,7)
3928 dnmin(icat)=dstprms(icat,8)
3929 dnmax(icat)=dstprms(icat,9)
3930 gnu(icat) = 4.0 ! shape distribution
3932 tabloln(icat)=log(tablo(icat))
3933 tabhiln(icat)=log(tabhi(icat))
3934 dict(icat) = float(ndn - 1) / (tabhiln(icat) - tabloln(icat))
3935 enddo
3936 rictmin=1.0001
3937 rictmax=0.9999*float(ndn)
3939 !.. Collection integral multipliers for number and mass mixing ratio loss
3940 !.. due to self collection from a cateogry. Multipliers for the gain term
3941 !.. for aggregates given as well. Below this are the multipliers for
3942 !.. cross collection. These should be correct but a more detailed
3943 !.. mathematical treatment is warranted. JYH Sept 2016.
3944 nselfloss = 1.
3945 rselfloss = 1.
3946 nselfgain = 1./2.
3947 rselfgain = 1.0
3949 ncrossloss = 1.
3950 rcrossloss = 1.0
3951 ncrossgain = 1./2.
3952 rcrossgain = 1.0
3954 !.. Assumed mass-size coefficients
3955 pwmas(3) = 3.0
3956 cfmas(3) = 3.14159/6.*250.*0.05
3957 pwmas(4) = 3.0
3958 cfmas(4) = 3.14159/6.*250.*0.05
3959 pwmas(5) = 3.0
3960 cfmas(5) = 3.14159/6.*50.*0.2
3962 dtlt = dtdum ! time step
3963 kp1 = 1
3964 kp2 = nz
3965 ipris= 5 ! double moment
3967 !.. collection routine
3968 !.. nz = # vertical levels
3969 !.. kp1 = lowest level with cloud
3970 !.. kp2 = higest level with cloud
3971 !.. 3 => the first is the collector
3972 !.. 3 => the second is the collectee
3973 !.. 5 => where the collection goes (aggregates)
3974 !.. dn(nz,ncat) = characteristic diameter
3975 !.. en(nz,ncat) = number concentration
3976 !.. r(nz,ncat) = mass mixing ratio
3977 !.. qr(nz,ncat) = hydrometeor internal enthalpy
3978 !.. qq(nz,ncat) = hydrometeor temperature (internal energy)
3979 !.. t(nz,ncat) = hydrometeor temperature (air temperature)
3980 !.. ict,wct1,wct2 -> table interpolation
3981 !.. rhoa(nz) = air density
3982 !.. denfac(nz) = sqrt(1/rhoa)
3983 !.. eff(nz) = collection efficiency
3984 !.. ipris = 5 double moment
3985 !.. rxfer(nz,ncat,ncat) = change in mass mixing ratio due to collection
3986 !.. enxfer(nz,ncat,ncat) = change in number concentration due to collection
3987 !.. qrxfer(nz,ncat,ncat) = change in internal energy due to collection
3988 !.. ipair,2 -> 2 is for self-collection, 1 if cross collection
3990 do k = 1,nz !.. nz = 1 currently (the k loop is done outsize of this routine)
3991 rhoa(k) = rhoairdum
3992 tk(k) = tempdum
3994 denfac(k) = sqrt(1./rhoa(k))
3995 dn(k,3) = ddum1
3996 en(k,3) = ndum1*rhoa(k)
3997 r(k,3) = qdum1
3999 dn(k,4) = ddum2
4000 en(k,4) = ndum2*rhoa(k)
4001 r(k,4) = qdum2
4003 qr(k,3) = 0.0 !.. internal energy... not needed so pass in zero
4004 qq(k,3) = t(k,3)
4005 t(k,3) = tk(k) - 273.15
4006 qr(k,4) = 0.0 !.. internal energy... not needed so pass in zero
4007 qq(k,4) = t(k,3)
4008 t(k,4) = t(k,3)
4009 qr(k,5) = 0.0
4010 qq(k,5) = t(k,5)
4011 t(k,5) = t(k,3)
4013 dn(k,5) = ddum3
4014 en(k,5) = ndum3*rhoa(k)
4015 r(k,5) = qdum3
4017 totaldc(k)=0.
4018 enddo
4020 do it = 1,1 !.. One timestep at a time
4021 do k = 1,nz !.. One vertical grid at a time
4023 do ihab = 1,7
4024 do jhab = 1,7
4026 !.. planar = qice, columnar=qice2, aggregates = qice3
4027 !.. planar + columnar = aggregates (qice + qice2 -> qice3)
4028 rxfer(k,ihab,jhab) = 0.0
4029 enxfer(k,ihab,jhab) = 0.0
4030 qrxfer(k,ihab,jhab) = 0.0
4031 rxfer(k,ihab,jhab) = 0.0
4032 enxfer(k,ihab,jhab) = 0.0
4033 qrxfer(k,ihab,jhab) = 0.0
4035 !.. planar + aggregates = aggregates
4036 parxfer(k,ihab,jhab) = 0.0
4037 paenxfer(k,ihab,jhab) = 0.0
4038 paqrxfer(k,ihab,jhab) = 0.0
4040 !.. columnar + aggregates = aggregates
4041 parxfer(k,ihab,jhab) = 0.0
4042 paenxfer(k,ihab,jhab) = 0.0
4043 paqrxfer(k,ihab,jhab) = 0.0
4045 !.. Self collection (aggregates)
4046 pprxfer(k,ihab,jhab) = 0.0
4047 ppenxfer(k,ihab,jhab) = 0.0
4048 ppqrxfer(k,ihab,jhab) = 0.0
4049 pprxfer(k,ihab,jhab) = 0.0
4050 ppenxfer(k,ihab,jhab) = 0.0
4051 ppqrxfer(k,ihab,jhab) = 0.0
4052 enddo
4053 enddo
4055 !.. needed for table interpolation
4056 rict=dict(3)*(log(max(1.e-10,dn(k,3)))-tabloln(3))+1.
4057 rictmm=max(rictmin,min(rictmax,rict))
4058 ict(k,3)=int(rictmm)
4059 wct2(k,3)=rictmm-float(ict(k,3))
4060 wct1(k,3)=1.0-wct2(k,3)
4062 !.. needed for table interpolation
4063 rict=dict(4)*(log(max(1.e-10,dn(k,4)))-tabloln(4))+1.
4064 rictmm=max(rictmin,min(rictmax,rict))
4065 ict(k,4)=int(rictmm)
4066 wct2(k,4)=rictmm-float(ict(k,4))
4067 wct1(k,4)=1.0-wct2(k,4)
4069 !.. needed for table interpolation
4070 rict=dict(5)*(log(max(1.e-10,dn(k,5)))-tabloln(5))+1.
4071 rictmm=max(rictmin,min(rictmax,rict))
4072 ict(k,5)=int(rictmm)
4073 wct2(k,5)=rictmm-float(ict(k,5))
4074 wct1(k,5)=1.0-wct2(k,5)
4076 totalc(k)=0.
4077 enddo
4079 !.. planar + columnar = aggregates
4080 !.. let high density reduce the collection efficiency
4081 rhoeffmax = max(rhodum1,rhodum2)
4082 phieffmax = max(min(phidum1,1./phidum1),min(phidum2,1./phidum2))
4083 if(rhoeffmax.le.400.) then
4084 rhoeff = 1.
4085 else
4086 rhoeff = -0.001923*rhoeffmax + 1.76916
4087 endif
4088 rhoeff = max(rhoeff,0.)
4089 rhoeff = min(rhoeff,1.)
4091 !.. Spherical particles collect less efficiently than more extreme
4092 !.. particle shapes
4093 !.. Based on Connolly et al 2012
4094 !.. See Jensen et al. 2018b
4095 if(phieffmax.le.0.03) then
4096 phieff = 1.
4097 elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then
4098 phieff = 0.0001*((phieffmax+0.07)**(-4.))
4099 else
4100 phieff = 0.
4101 endif
4103 if(phieff.le.0.001) then
4104 phieff = 0.
4105 endif
4107 phieff = max(phieff,0.01) !.. Low limit currently set to 0.01 AAJ
4108 phieff = min(phieff,1.)
4109 efffact = phieff * rhoeff
4111 call col1 &
4112 (nz,kp1,kp2,ndn,ncat,3,4,5 &
4113 ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab &
4114 ,coltabn,ipair,1.,eff,efffact,ipris,rxfer,qrxfer,enxfer)
4115 call col1 &
4116 (nz,kp1,kp2,ndn,ncat,4,3,5 &
4117 ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab &
4118 ,coltabn,ipair,1.,eff,efffact,ipris,rxfer,qrxfer,enxfer)
4120 !.. plates + aggregates = aggregates
4121 rhoeffmax = rhodum1
4122 phieffmax = min(phidum1,1./phidum1)
4124 if(rhoeffmax.le.400.) then
4125 rhoeff = 1.
4126 else
4127 rhoeff = -0.001923*rhoeffmax + 1.76916
4128 endif
4129 rhoeff = max(rhoeff,0.)
4130 rhoeff = min(rhoeff,1.)
4132 !.. Based on Connolly et al 2012
4133 !.. See Jensen et al. 2018b
4134 if(phieffmax.le.0.03) then
4135 phieff = 1.
4136 elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then
4137 phieff = 0.0001*((phieffmax+0.07)**(-4.))
4138 else
4139 phieff = 0.
4140 endif
4142 if(phieff.le.0.001) then
4143 phieff = 0.
4144 endif
4146 phieff = max(phieff,0.01)
4147 phieff = min(phieff,1.)
4148 efffact = phieff * rhoeff
4150 call col1 &
4151 (nz,kp1,kp2,ndn,ncat,3,5,5 &
4152 ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab &
4153 ,coltabn,ipair,1.,eff,efffact,ipris,parxfer,paqrxfer,paenxfer)
4155 !.. Columnar + Aggregates = Aggregates
4156 rhoeffmax = rhodum2
4157 phieffmax = min(phidum2,1./phidum2)
4159 if(rhoeffmax.le.400.) then
4160 rhoeff = 1.
4161 else
4162 rhoeff = -0.001923*rhoeffmax + 1.76916
4163 endif
4164 rhoeff = max(rhoeff,0.)
4165 rhoeff = min(rhoeff,1.)
4167 !.. Based on Connolly et al 2012
4168 !.. See Jensen et al. 2018b
4169 if(phieffmax.le.0.03) then
4170 phieff = 1.
4171 elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then
4172 phieff = 0.0001*((phieffmax+0.07)**(-4.))
4173 else
4174 phieff = 0.
4175 endif
4177 if(phieff.le.0.001) then
4178 phieff = 0.
4179 endif
4181 phieff = max(phieff,0.01)
4182 phieff = min(phieff,1.)
4183 efffact = phieff * rhoeff
4185 call col1 &
4186 (nz,kp1,kp2,ndn,ncat,4,5,5 &
4187 ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab &
4188 ,coltabn,ipair,1.,eff,efffact,ipris,parxfer,paqrxfer,paenxfer)
4190 !.. Plates + Plates = Aggregates
4191 rhoeffmax = rhodum1
4192 phieffmax = min(phidum1,1./phidum1)
4194 if(rhoeffmax.le.400.) then
4195 rhoeff = 1.
4196 else
4197 rhoeff = -0.001923*rhoeffmax + 1.76916
4198 endif
4199 rhoeff = max(rhoeff,0.)
4200 rhoeff = min(rhoeff,1.)
4202 !.. Based on Connolly et al 2012
4203 !.. See Jensen et al. 2018b
4204 if(phieffmax.le.0.03) then
4205 phieff = 1.
4206 elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then
4207 phieff = 0.0001*((phieffmax+0.07)**(-4.))
4208 else
4209 phieff = 0.
4210 endif
4212 if(phieff.le.0.001) then
4213 phieff = 0.
4214 endif
4216 phieff = max(phieff,0.01)
4217 phieff = min(phieff,1.)
4218 efffact = phieff * rhoeff
4220 call col1 &
4221 (nz,kp1,kp2,ndn,ncat,3,3,5 &
4222 ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab &
4223 ,coltabn,ipair,1.,eff,efffact,ipris,pprxfer,ppqrxfer,ppenxfer)
4225 !.. Columns + Columns = Aggregates
4226 rhoeffmax = rhodum2
4227 phieffmax = min(phidum2,1./phidum2)
4229 if(rhoeffmax.le.400.) then
4230 rhoeff = 1.
4231 else
4232 rhoeff = -0.001923*rhoeffmax + 1.76916
4233 endif
4234 rhoeff = max(rhoeff,0.)
4235 rhoeff = min(rhoeff,1.)
4237 !.. Based on Connolly et al 2012
4238 !.. See Jensen et al. 2018b
4239 if(phieffmax.le.0.03) then
4240 phieff = 1.
4241 elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then
4242 phieff = 0.0001*((phieffmax+0.07)**(-4.))
4243 else
4244 phieff = 0.
4245 endif
4247 if(phieff.le.0.001) then
4248 phieff = 0.
4249 endif
4251 phieff = max(phieff,0.01)
4252 phieff = min(phieff,1.)
4253 efffact = phieff * rhoeff
4255 call col1 &
4256 (nz,kp1,kp2,ndn,ncat,4,4,5 &
4257 ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab &
4258 ,coltabn,ipair,1.,eff,efffact,ipris,pprxfer,ppqrxfer,ppenxfer)
4260 efffact=1.
4261 !.. Aggregate self collection
4262 call col1 &
4263 (nz,kp1,kp2,ndn,ncat,5,5,5 &
4264 ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab &
4265 ,coltabn,ipair,1.,eff,efffact,ipris,pprxfer,ppqrxfer,ppenxfer)
4267 do k = 1,nz
4269 !.. Do not over-deplete from aggregation
4270 sink3 = rcrossloss*rxfer(k,3,5) + &
4271 rcrossloss*parxfer(k,3,5) + &
4272 rselfloss*pprxfer(k,3,5)
4273 source3 = r(k,3)
4275 sink4 = rcrossloss*rxfer(k,4,5) + &
4276 rcrossloss*parxfer(k,4,5) + &
4277 rselfloss*pprxfer(k,4,5)
4278 source4 = r(k,4)
4280 if(sink3.gt.source3.and.source3.gt.1.e-8) then
4281 ratioagg= source3/sink3
4282 rxfer(k,3,5) = rxfer(k,3,5)*ratioagg
4283 parxfer(k,3,5) = parxfer(k,3,5)*ratioagg
4284 pprxfer(k,3,5) = pprxfer(k,3,5)*ratioagg
4285 endif
4287 if(sink4.gt.source4.and.source4.gt.1.e-8) then
4288 ratioagg= source4/sink4
4289 rxfer(k,4,5) = rxfer(k,4,5)*ratioagg
4290 parxfer(k,4,5) = parxfer(k,4,5)*ratioagg
4291 pprxfer(k,4,5) = pprxfer(k,4,5)*ratioagg
4292 endif
4294 totaldc(k)=totaldc(k) + &
4295 ncrossloss*enxfer(k,3,5)+ncrossloss*paenxfer(k,3,5)+ &
4296 nselfloss*ppenxfer(k,3,5)+ &
4297 ncrossloss*enxfer(k,4,5)+ncrossloss*paenxfer(k,4,5)+ &
4298 nselfloss*ppenxfer(k,4,5)
4300 !.. qagg is the mass tendency, nagg is the number tendency
4301 qagg1 = -rcrossloss*rxfer(k,3,5) - rcrossloss*parxfer(k,3,5) - rselfloss*pprxfer(k,3,5)
4302 nagg1 = -ncrossloss*enxfer(k,3,5) - ncrossloss*paenxfer(k,3,5) - nselfloss*ppenxfer(k,3,5)
4304 qagg2 = -rcrossloss*rxfer(k,4,5) - rcrossloss*parxfer(k,4,5) - rselfloss*pprxfer(k,4,5)
4305 nagg2 = -ncrossloss*enxfer(k,4,5) - ncrossloss*paenxfer(k,4,5) - nselfloss*ppenxfer(k,4,5)
4307 qagg3 = rcrossgain*rxfer(k,3,5) + rcrossgain*rxfer(k,4,5) + rcrossgain*parxfer(k,3,5) + &
4308 rcrossgain*parxfer(k,4,5) + rselfgain*pprxfer(k,3,5) + rselfgain*pprxfer(k,4,5)
4309 nagg3 = ncrossgain*enxfer(k,3,5) + ncrossgain*enxfer(k,4,5) + nselfgain*ppenxfer(k,3,5) + &
4310 nselfgain*ppenxfer(k,4,5) - 0.5*ppenxfer(k,5,5)
4312 qagg1 = min(qagg1,0.)
4313 qagg2 = min(qagg2,0.)
4314 qagg3 = max(qagg3,0.)
4316 !.. Update to characteristic diameters after aggregation
4317 !.. Note that for ice-one and ice-two, the axes will be re-diagnosed
4318 !.. after aggregation based on delta*, so only output aggregation size
4319 dn(k,5) = ((r(k,5)+qagg3)*rhoa(k)/(MAX((en(k,5)+nagg3),QNSMALL)*cfmas(5) &
4320 *(gamma(gnu(5)+ &
4321 pwmas(5)))/(gamma(gnu(5))) ) )**(1./pwmas(5))
4322 ddum3=dn(k,5)
4324 nagg1 = nagg1/rhoa(k)
4325 nagg2 = nagg2/rhoa(k)
4326 nagg3 = nagg3/rhoa(k)
4328 totalc(k)=0.
4329 totalc(k) = en(k,3)+en(k,4) + 0.
4331 rxfer(k,3,5) = 0.0
4332 enxfer(k,3,5) = 0.0
4333 qrxfer(k,3,5) = 0.0
4335 rxfer(k,4,5) = 0.0
4336 enxfer(k,4,5) = 0.0
4337 qrxfer(k,4,5) = 0.0
4338 enddo
4339 enddo
4341 end subroutine aggregation
4343 !--------------------------------------------------------------------------------------------------------------!
4344 !--------------------------------------------------------------------------------------------------------------!
4345 !.. Needed for the aggregation routine
4346 subroutine col1(n1,k1,k2,ndn,ncat,mx,my,mz,dtlt,dn,en,r,qr,qq,t &
4347 ,ict,wct1,wct2,dn0,denfac,coltab,coltabn,ipair,cr,eff,efdum,icf &
4348 ,rxfer,qrxfer,enxfer)
4349 integer ict(n1,*),ipair(ncat,ncat), ndn, ncat,n1
4350 real dn(n1,*),en(n1,*),r(n1,*),qr(n1,*),qq(n1,*),t(n1,*) &
4351 ,wct1(n1,*),wct2(n1,*),dn0(*),denfac(*) &
4352 ,coltab(ndn,ndn,*),eff(*),coltabn(ndn,ndn,*) &
4353 ,rxfer(n1,ncat,ncat),qrxfer(n1,ncat,ncat),enxfer(n1,ncat,ncat)
4354 integer :: k1, k2, mx, my, mz, icf, k, ip
4355 real :: dtlt, cr, csizei, tmax, colamt, colamtn, wght, deltan,efdum
4357 ip=ipair(mx,my)
4358 csizei=7.e4
4360 do k=k1,k2
4362 tmax=max(t(k,mx),t(k,my))
4364 !.. efdum added by AAJ from ramsagg
4365 !.. if takes into account the shape and density
4366 !.. dependence on the collection efficiency
4367 eff(k)=min(0.2,10.**(0.035*tmax-0.7))*efdum
4369 if(my.eq.6.and.qr(k,6).gt.0)eff(k)=1.
4370 if(my.eq.7.and.qr(k,7).gt.0)eff(k)=1.
4372 !.. see note above on efdum
4373 if(mz.eq.5.and.abs(t(k,mx)+14.).le.2.)eff(k)=1.4*efdum
4375 if(mx.eq.1)eff(k)=min(1.,dn(k,1)*csizei)
4377 colamt=dtlt*cr*0.785*eff(k)*denfac(k)/dn0(k) &
4378 *en(k,mx)*en(k,my)*( &
4379 wct1(k,mx)*wct1(k,my)*coltab(ict(k,mx) ,ict(k,my) ,ip) &
4380 +wct2(k,mx)*wct1(k,my)*coltab(ict(k,mx)+1,ict(k,my) ,ip) &
4381 +wct1(k,mx)*wct2(k,my)*coltab(ict(k,mx) ,ict(k,my)+1,ip) &
4382 +wct2(k,mx)*wct2(k,my)*coltab(ict(k,mx)+1,ict(k,my)+1,ip))
4384 colamt=min(colamt,r(k,mx))
4385 rxfer(k,mx,mz) = rxfer(k,mx,mz) + colamt
4386 qrxfer(k,mx,mz) = qrxfer(k,mx,mz) + colamt * qq(k,mx)
4388 colamtn=dtlt*cr*0.785*eff(k)*denfac(k)/dn0(k) &
4389 *en(k,mx)*en(k,my)*( &
4390 wct1(k,mx)*wct1(k,my)*coltabn(ict(k,mx) ,ict(k,my) ,ip) &
4391 +wct2(k,mx)*wct1(k,my)*coltabn(ict(k,mx)+1,ict(k,my) ,ip) &
4392 +wct1(k,mx)*wct2(k,my)*coltabn(ict(k,mx) ,ict(k,my)+1,ip) &
4393 +wct2(k,mx)*wct2(k,my)*coltabn(ict(k,mx)+1,ict(k,my)+1,ip))
4395 colamtn=min(colamtn,en(k,mx))
4397 if(icf.eq.5)then
4399 wght = min(max(0.,100.*abs(colamt)/max(1.e-20,r(k,mx))),1.)
4400 deltan = colamtn*(1.-wght) + &
4401 wght*en(k,mx) * colamt / max(1.e-20,r(k,mx))
4402 enxfer(k,mx,mz) = enxfer(k,mx,mz) + deltan
4404 endif
4406 enddo
4407 return
4408 end subroutine col1
4410 !--------------------------------------------------------------------------------------------------------------!
4411 !--------------------------------------------------------------------------------------------------------------!
4412 !.. Needed for the aggregation routine
4413 subroutine mkcoltb(ndn,ncat,coltab,coltabn,ipair,gnu,tablo,tabhi &
4414 ,cfmas,pwmas,cfvt,pwvt)
4415 integer :: ncat
4416 real :: coltab(ndn,ndn,35),gnu(ncat) &
4417 ,tablo(ncat), tabhi(ncat), cfmas(ncat),pwmas(ncat) &
4418 ,cfvt(ncat), pwvt(ncat), coltabn(ndn,ndn,35)
4419 integer :: ipair(ncat,ncat), ix, iy, idny, idnx, idx, ndn
4420 integer, parameter :: ndx=20
4421 real :: dx(ndx),gx(ndx),fx(20), gyn,gyn1,gyn2,gynp1,gynp2,gynp
4422 real :: dxhi, dny, vny, dxlo, dnx, ans
4424 do ix=1,ncat
4425 do iy=1,ncat
4426 if(ipair(ix,iy).gt.0)then
4427 gyn=(gamma(gnu(iy)))
4428 gyn1=(gamma(gnu(iy)+1.))/gyn
4429 gyn2=(gamma(gnu(iy)+2.))/gyn
4430 gynp=(gamma(gnu(iy)+pwvt(iy)))/gyn
4431 gynp1=(gamma(gnu(iy)+pwvt(iy)+1.))/gyn
4432 gynp2=(gamma(gnu(iy)+pwvt(iy)+2.))/gyn
4433 dxlo=tablo(ix)*0.01
4434 dxhi=tabhi(ix)*10.0
4435 do idny=1,ndn
4436 dny=tablo(iy)*(tabhi(iy)/tablo(iy))**(float(idny-1) &
4437 /float(ndn-1))
4438 vny=cfvt(iy)*dny**pwvt(iy)
4439 do idnx=1,ndn
4440 dnx=tablo(ix)*(tabhi(ix)/tablo(ix))**(float(idnx-1) &
4441 /float(ndn-1))
4442 do idx=1,ndx
4443 dx(idx)=dxlo*(dxhi/dxlo) &
4444 **(float(idx-1)/float(ndx-1))
4445 fx(idx) = xjnum(dx(idx),cfvt(ix),pwvt(ix),cfvt(iy) &
4446 ,pwvt(iy),vny,dnx,dny,gnu(ix),gnu(iy) &
4447 ,gyn1,gyn2,gynp,gynp1,gynp2)
4448 if(fx(idx).lt.1.e-15) then !.. low limit added by AAJ
4449 fx(idx) = 0.
4450 endif
4451 gx(idx) = fx(idx) * cfmas(ix) &
4452 * dx(idx) ** pwmas(ix)
4453 enddo
4454 call avint(dx,gx,ndx,dxlo,dxhi,ans)
4455 coltab(idnx,idny,ipair(ix,iy))=max(0.,ans)
4456 call avint(dx,fx,ndx,dxlo,dxhi,ans)
4457 coltabn(idnx,idny,ipair(ix,iy))=max(0.,ans)
4458 enddo
4459 enddo
4460 endif
4461 enddo
4462 enddo
4463 return
4464 end subroutine mkcoltb
4466 !--------------------------------------------------------------------------------------------------------------!
4467 !--------------------------------------------------------------------------------------------------------------!
4468 !.. Needed for the aggregation routine
4469 function xjnum(dx,cvx,pvx,cvy,pvy,vny,dnx,dny,xnu,ynu &
4470 ,gyn1,gyn2,gynp,gynp1,gynp2)
4471 implicit none
4473 real :: dx,cvx,pvx,cvy,pvy,vny,dnx,dny,xnu,ynu,gyn1,gyn2,gynp, &
4474 gynp1,gynp2,xjnum &
4475 ,dnxi,rdx,vx,dxy,ynup
4477 dnxi = 1. / dnx
4478 rdx = dx * dnxi
4479 vx = max((cvx * dx ** pvx),1.e-6) !.. low limit added by AAJ
4480 dxy = (vx / cvy) ** (1. / pvy) / dny
4481 dxy = max(dxy,1.e-5)
4482 dxy = min(dxy,70.)
4483 ynup = ynu + pvy
4485 if (rdx .lt. 38.) then
4486 xjnum=exp(-rdx-gammln(xnu)-gammln(ynu))*rdx**(xnu-1.)*dnxi*( &
4487 vx*(dx*dx*(gammap(ynu,dxy)-gammq(ynu,dxy)) &
4488 +2.*dx*dny*gyn1*(gammap(ynu+1.,dxy)-gammq(ynu+1.,dxy)) &
4489 +dny*dny*gyn2*(gammap(ynu+2.,dxy)-gammq(ynu+2.,dxy))) &
4490 -vny*(dx*dx*gynp*(gammap(ynup,dxy)-gammq(ynup,dxy)) &
4491 +2.*dx*dny*gynp1*(gammap(ynup+1.,dxy)-gammq(ynup+1.,dxy)) &
4492 +dny*dny*gynp2*(gammap(ynup+2.,dxy)-gammq(ynup+2.,dxy))))
4493 else
4494 xjnum=0.
4495 endif
4496 return
4497 end function xjnum
4499 !--------------------------------------------------------------------------------------------------------------!
4500 !--------------------------------------------------------------------------------------------------------------!
4501 !.. Used by aggregation
4502 REAL FUNCTION GAMMQ(A,X)
4504 real :: x, a, gln, gamser
4506 IF (X .LT. A+1.) THEN
4507 CALL LOWGSERIES(GAMSER,A,X,GLN)
4508 GAMMQ = 1. - GAMSER
4509 ELSE
4510 CALL HIGHGCONTFRAC(GAMMQ,A,X,GLN)
4511 END IF
4512 RETURN
4513 END FUNCTION GAMMQ
4515 !--------------------------------------------------------------------------------------------------------------!
4516 !--------------------------------------------------------------------------------------------------------------!
4517 !.. Computes P(a, x) = lower_gamma(a, x) / gamma(a)
4518 !.. Using either series solution or continued fractions
4519 real function gammap(a, x)
4521 implicit none
4523 REAL :: a, x, gammaser, gammacf, gln
4525 if(x.lt.0.0.or.a.lt.0.0) then
4526 write(*,*) 'either x or a is less than zero'
4527 write(*,*) x, a
4528 return
4529 endif
4530 if(x.lt.(a + 1.0)) then
4531 call lowgseries(gammaser, a, x, gln)
4532 gammap = gammaser
4533 else
4534 call highgcontfrac(gammacf, a, x, gln)
4535 gammap = 1.0 - gammacf
4536 endif
4538 return
4539 end function gammap
4541 !--------------------------------------------------------------------------------------------------------------!
4542 !--------------------------------------------------------------------------------------------------------------!
4543 !.. Returns: gln = ln(gamma(a))
4544 !.. Returns: gammaser = lower_gamma(a, x) as a series solution
4545 subroutine lowgseries(gammaser, a, x, gln)
4547 implicit none
4549 INTEGER, PARAMETER :: ITMAX = 100
4550 INTEGER :: n
4551 REAL :: gammaser, a, x, gln, sum, del, ap
4552 REAL, PARAMETER :: EPS = 3.0e-7
4554 gln = gammln(a)
4555 if(x.le.0) then
4556 gammaser = 0.0
4557 return
4558 endif
4560 ap=a
4561 sum=1./a
4562 del=sum
4563 do n = 1, ITMAX
4564 ap=ap+1.
4565 del=del*x/ap
4566 sum=sum+del
4567 if(abs(del).lt.abs(sum)*EPS) exit
4568 enddo
4570 gammaser = sum*exp(-x+a*log(x)-gln)
4571 return
4572 end subroutine lowgseries
4574 !--------------------------------------------------------------------------------------------------------------!
4575 !--------------------------------------------------------------------------------------------------------------!
4576 !.. Returns: gln = ln(gamma(a))
4577 !.. Returns: gammacf = upper_gamma(a, x) as a continued fractions solution
4578 subroutine highgcontfrac(gammacf, a, x, gln)
4580 implicit none
4582 INTEGER, PARAMETER :: ITMAX = 100
4583 INTEGER :: n
4584 REAL :: gammacf, a, x, gln, b, c, d, h, del, an
4585 REAL, PARAMETER :: EPS = 3.0e-7, FPMIN = 1.0e-30
4587 gln = gammln(a)
4588 b=x+1.-a
4589 c=1./FPMIN
4590 d=1./b
4591 h=d
4592 do n = 1, ITMAX
4593 an=-n*(n-a)
4594 b=b+2.
4595 d=an*d+b
4596 if(abs(d).lt.FPMIN) d = FPMIN
4597 c=b+an/c
4598 if(abs(c).lt.FPMIN) c = FPMIN
4599 d = 1./d
4600 del=d*c
4601 h = h*del
4602 if(abs(del-1.).lt.EPS) exit
4603 enddo
4605 gammacf = exp(-x+a*log(x)-gln)*h
4606 return
4607 end subroutine highgcontfrac
4609 !--------------------------------------------------------------------------------------------------------------!
4610 !--------------------------------------------------------------------------------------------------------------!
4611 !.. Used by aggregation
4612 SUBROUTINE AVINT(X,Y,N,XLO,XUP,ANS)
4613 DOUBLE PRECISION :: R3,RP5,SUM,SYL,SYL2,SYL3,SYU,SYU2,SYU3,X1,X2,X3 &
4614 ,X12,X13,X23,TERM1,TERM2,TERM3,A,B,C,CA,CB,CC
4615 REAL :: XLO,XUP,ANS,FL,FR,SLOPE
4616 INTEGER :: N, I, INLFT,INRT,ISTART,ISTOP
4617 REAL, DIMENSION(N) :: X,Y
4619 ANS =0.0
4620 IF (XLO-XUP) 3,100,200
4621 3 IF (N.LT.2) GO TO 215
4622 DO 5 I=2,N
4623 IF (X(I).LE.X(I-1)) GO TO 210
4624 IF (X(I).GT.XUP) GO TO 6
4625 5 CONTINUE
4626 6 CONTINUE
4627 IF (N.GE.3) GO TO 9
4628 !.. special n=2 case
4629 SLOPE = (Y(2)-Y(1))/(X(2)-X(1))
4630 FL = Y(1) + SLOPE*(XLO-X(1))
4631 FR = Y(2) + SLOPE*(XUP-X(2))
4632 ANS = 0.5*(FL+FR)*(XUP-XLO)
4633 RETURN
4634 9 CONTINUE
4635 IF (X(N-2).LT.XLO) GO TO 205
4636 IF (X(3).GT.XUP) GO TO 205
4637 I = 1
4638 10 IF (X(I).GE.XLO) GO TO 15
4639 I = I+1
4640 GO TO 10
4641 15 INLFT = I
4642 I = N
4643 20 IF (X(I).LE.XUP) GO TO 25
4644 I = I-1
4645 GO TO 20
4646 25 INRT = I
4647 IF ((INRT-INLFT).LT.2) GO TO 205
4648 ISTART = INLFT
4649 IF (INLFT.EQ.1) ISTART = 2
4650 ISTOP = INRT
4651 IF (INRT.EQ.N) ISTOP = N-1
4653 R3 = 3.0D0
4654 RP5= 0.5D0
4655 SUM = 0.0
4656 SYL = XLO
4657 SYL2= SYL*SYL
4658 SYL3= SYL2*SYL
4660 DO 50 I=ISTART,ISTOP
4661 X1 = X(I-1)
4662 X2 = X(I)
4663 X3 = X(I+1)
4664 X12 = X1-X2
4665 X13 = X1-X3
4666 X23 = X2-X3
4667 TERM1 = DBLE(Y(I-1))/(X12*X13)
4668 TERM2 =-DBLE(Y(I)) /(X12*X23)
4669 TERM3 = DBLE(Y(I+1))/(X13*X23)
4670 A = TERM1+TERM2+TERM3
4671 B = -(X2+X3)*TERM1 - (X1+X3)*TERM2 - (X1+X2)*TERM3
4672 C = X2*X3*TERM1 + X1*X3*TERM2 + X1*X2*TERM3
4673 IF (I-ISTART) 30,30,35
4674 30 CA = A
4675 CB = B
4676 CC = C
4677 GO TO 40
4678 35 CA = 0.5*(A+CA)
4679 CB = 0.5*(B+CB)
4680 CC = 0.5*(C+CC)
4681 40 SYU = X2
4682 SYU2= SYU*SYU
4683 SYU3= SYU2*SYU
4684 SUM = SUM + CA*(SYU3-SYL3)/R3 + CB*RP5*(SYU2-SYL2) + CC*(SYU-SYL)
4685 CA = A
4686 CB = B
4687 CC = C
4688 SYL = SYU
4689 SYL2= SYU2
4690 SYL3= SYU3
4691 50 CONTINUE
4692 SYU = XUP
4693 ANS = SUM + CA*(SYU**3-SYL3)/R3 + CB*RP5*(SYU**2-SYL2) &
4694 + CC*(SYU-SYL)
4695 100 RETURN
4696 200 PRINT*, 'Upper limit of integration not greater than lower limit.'
4697 STOP 'AVINT2'
4698 205 PRINT*, 'Less than 3 function values between integration limits.'
4699 STOP 'AVINT3'
4700 210 PRINT*, 'Abscissas not strictly increasing.'
4701 STOP 'AVINT4'
4702 215 PRINT*, 'Less than 2 function values were supplied.'
4703 STOP 'AVINT5'
4704 END SUBROUTINE
4706 !--------------------------------------------------------------------------------------------------------------!
4707 !--------------------------------------------------------------------------------------------------------------!
4709 end module module_mp_jensen_ishmael
4711 !.. Full stop