| blob | 82194c6ef037b3723975ca0f448880be57ff889b |
1 PROGRAM create_p3_lookuptable_1
3 !______________________________________________________________________________________
4 !
5 ! This program creates the lookup tables (that are combined into the single file
6 ! 'p3_lookupTable_1.dat') for the microphysical processes, used by the P3 microphysics scheme.
7 ! Note, 'p3_lookupTable_2.dat', used in ice-ice interactions for the multi-ice-category
8 ! configuration, are created by a separate program, 'create_lookupTable_2.f90').
9 !
10 ! This code is used to generate lookupTable_1 for either 2momI or 3momI, depending on
11 ! the specified value of the parameter 'log_3momI'.
12 !
13 ! All other parameter settings are linked uniquely to the version number.
14 !
15 !--------------------------------------------------------------------------------------
16 ! Version: 5.4
17 ! Last modified: 2022 June
18 !______________________________________________________________________________________
20 !______________________________________________________________________________________
21 !
22 ! To generate 'p3_lookupTable_1.dat' using this code, do the following steps :
23 !
24 ! 1. Break up this code into two parts (-top.f90 and -bottom.90). Search for the string
25 ! 'RUNNING IN PARALLEL MODE' and follow instructions.
26 ! For 2momI, need to comment do/enddo for 'i_rhor_loop' for parallelization (search for
27 ! 'i_rhor_loop' label); for 3momI, the 'i_rhor_loop' needs to be uncomment.
28 ! (In the future, a script will be written to automate this.)
29 !
30 ! 2. Copy the 3 pieces of text below to create indivudual executable shell scripts.
32 ! 3. Run script 1 (./go_1-compile.ksh). This script will recreate several
33 ! versions of the full code, concatenating the -top.f90 and the -bottom.90 with
34 ! the following code (e.g.) in between (looping through values of i_rhor (for 2momI) or
35 ! i_Znorm (for 3momI):
36 !
37 ! i_rhor = 1 ! 2-moment-ice
38 ! i_Znorm = 1 ! 3-moment-ice
39 !
40 ! Each version of full_code.f90 is then compiled, with a unique executable name.
41 ! Note, this is done is place the outer i1 loop in order to parallelized .
42 !
43 ! 4. Run script 2 (./go_2-submit.csh) This create temporary work directories,
44 ! moves each executable to each work directory, and runs the executables
45 ! simultaneously. (Note, it is assumed that the machine on which this is done
46 ! has multiple processors, though it is not necessary.)
47 !
48 ! 5. Run script 3 (./go_3-concatenate.ksh). This concatenates the individual
49 ! partial tables (the output in each working directory) into a single, final
50 ! file 'p3_lookupTable_1.dat' Once it is confirmed that the table is correct,
51 ! the work directories and their contents can be removed.
52 !
53 ! Note: For testing or to run in serial, ompile with double-precision
54 ! (e.g. ifort -r8 create_p3_lookupTable_1.f90)
55 ! gfortran -fdefault-real-8 create_p3_lookupTable_1.f90)
56 !______________________________________________________________________________________
58 !--------------------------------------------------------------------------------------------
59 ! Parallel script 1 (of 3): [copy text below (uncommented) to file 'go_1-compile.ksh']
60 ! - creates individual parallel codes, compiles each
61 !
62 !-----------------------------
63 ! For 2-MOMENT-ICE:
65 ! #!/bin/ksh
66 !
67 ! for i_rhor in 01 02 03 04 05
68 ! do
69 !
70 ! rm cfg_input full_code.f90
71 ! cat > cfg_input << EOF
72 ! i_rhor = ${i_rhor}
73 ! EOF
74 ! cat create_p3_lookupTable_1-top.f90 cfg_input create_p3_lookupTable_1-bottom.f90 > full_code.f90
75 ! echo 'Compiling 'exec_${i_rhor}
76 ! ifort -r8 full_code.f90
77 ! mv a.out exec_${i_rhor}
78 !
79 ! done
80 !
81 ! rm cfg_input full_code.f90
83 !-----------------------------
84 ! For 3-MOMENT-ICE:
86 !#!/bin/ksh
87 !
88 ! for i_Znorm in 01 02 03 04 05 06 07 08 09 10 11
89 ! do
90 !
91 ! rm cfg_input full_code.f90
92 ! cat > cfg_input << EOF
93 ! i_Znorm = ${i_Znorm}
94 ! EOF
95 ! cat create_p3_lookupTable_1-top.f90 cfg_input create_p3_lookupTable_1-bottom.f90 > full_code.f90
96 ! echo 'Compiling 'exec_${i_Znorm}
97 ! ifort -r8 full_code.f90
98 ! mv a.out exec_${i_Znorm}
99 !
100 ! done
101 !
102 ! rm cfg_input full_code.f90
104 !--------------------------------------------------------------------------------------------
105 !# Parallel script 2 (of 3): [copy text below (uncommented) to file 'go_2-submit.ksh']
106 !# - creates individual work directories, launches each executable
108 !#!/bin/ksh
109 !
110 ! for exec in `ls exec_*`
111 ! do
112 ! echo Submitting: ${exec}
113 ! mkdir ${exec}-workdir
114 ! mv ${exec} ${exec}-workdir
115 ! cd ${exec}-workdir
116 ! ./${exec} > log &
117 ! cd ..
118 ! done
120 !--------------------------------------------------------------------------------------------
121 !# Parallel script 3 (of 3): [copy text below (uncommented) to file 'go_3-concatenate.ksh]
122 !# - concatenates the output of each parallel job into a single output file.
124 !#!/bin/ksh
125 !
126 ! rm lt_total
127 !
128 ! for i in `ls exec*/*dat`
129 ! do
130 ! echo $i
131 ! cat lt_total $i > lt_total_tmp
132 ! mv lt_total_tmp lt_total
133 ! done
134 !
135 ! mv lt_total p3_lookupTable_1.dat
136 !
137 ! echo 'Done. Work directories and contents can now be removed.'
138 ! echo 'Be sure to re-name the file with the appropriate extension, with the version number'
139 ! echo 'corresponding to that in the header. (e.g. 'p3_lookupTable_1.dat-v5.3-3momI')'
141 !--------------------------------------------------------------------------------------------
143 implicit none
145 !-----
146 character(len=20), parameter :: version = '5.4'
147 logical, parameter :: log_3momI = .true. !switch to create table for 2momI (.false.) or 3momI (.true.)
148 !-----
150 integer :: i_Znorm ! index for normalized (by Q) Z (passed in through script; [1 .. n_Znorm])
151 integer :: i_rhor ! index for rho_rime (passed in through script; [1 .. n_rhor])
152 integer :: i_Fr ! index for rime-mass-fraction loop [1 .. n_Fr]
153 integer :: i_Qnorm ! index for normalized (by N) Q loop [1 .. n_Qnorm]
154 integer :: i_Drscale ! index for scaled mean rain size loop [1 .. n_Drscale]
156 ! NOTE: n_Znorm (number of i_Znorm values) is currently equal to 80. It is not actually used herein and therefore not declared.
157 ! Rather, the outer (i_Znorm) "loop" is treated by individual complilations/execuations with specified values of i_Znorm,
158 ! with resulting sub-tables subsequently concatenated. The same is true for the second (i_rhor) "loop"; however, n_rhor
159 ! is used to decare the ranges of other arrays, hence it is declared/initialized here
161 integer, parameter :: n_Znorm = 10 ! number of indices for i_Znorm loop (1nd "loop") [not used in parallelized version]
162 integer, parameter :: n_rhor = 5 ! number of indices for i_rhor loop (2nd "loop")
163 integer, parameter :: n_Fr = 4 ! number of indices for i_Fr loop (3rd loop)
164 integer, parameter :: n_Qnorm = 50 ! number of indices for i_Qnorm loop (4th loop)
165 integer, parameter :: n_Drscale = 30 ! number of indices for scaled mean rain size (5th [inner] loop)
167 integer, parameter :: num_int_bins = 40000 ! number of bins for numerical integration of ice processes
168 integer, parameter :: num_int_coll_bins = 1500 ! number of bins for numerical integration of ice-ice and ice-rain collection
170 real, parameter :: mu_i_min = 0.
171 real, parameter :: mu_i_max = 20.
173 integer :: i,ii,iii,jj,kk,kkk,dumii,i_iter,n_iter_psdSolve
175 real :: N,q,qdum,dum1,dum2,cs1,ds1,lam,n0,lamf,qerror,del0,c0,c1,c2,dd,ddd,sum1,sum2, &
176 sum3,sum4,xx,a0,b0,a1,b1,dum,bas1,aas1,aas2,bas2,gammq,gamma,d1,d2,delu,lamold, &
177 cap,lamr,dia,amg,dv,n0dum,sum5,sum6,sum7,sum8,dg,cg,bag,aag,dcritg,dcrits, &
178 dcritr,Fr,csr,dsr,duml,dum3,rhodep,cgpold,m1,m2,m3,dt,mu_r,initlamr,lamv, &
179 rdumii,lammin,lammax,pi,g,p,t,rho,mu,mu_i,ds,cs,bas,aas,dcrit,mu_dum,gdum, &
180 Z_value,sum9,mom3,mom6,intgrR1,intgrR2,intgrR3,intgrR4,dum4,cs2,ds2
182 ! function to compute mu for triple moment
183 real :: compute_mu_3moment
185 ! function to return diagnostic value of shape paramter, mu_i
186 real :: diagnostic_mui
188 ! outputs from lookup table (i.e. "f1prxx" read by access_lookup_table in s/r p3_main)
189 real, dimension(n_Qnorm,n_Fr) :: uns,ums,refl,dmm,rhomm,nagg,nrwat,qsave,nsave,vdep, &
190 eff,lsave,a_100,n_100,vdep1,i_qsmall,i_qlarge,nrwats,lambda_i,mu_i_save
192 ! outputs for triple moment
193 ! HM zsmall, zlarge no longer needed
194 ! real, dimension(n_Qnorm,n_Fr) :: uzs,zlarge,zsmall
195 real, dimension(n_Qnorm,n_Fr) :: uzs
197 real, dimension(n_Qnorm,n_Drscale,n_Fr) :: qrrain,nrrain,nsrain,qsrain,ngrain
199 real, dimension(n_Drscale) :: lamrs
200 real, dimension(num_int_bins) :: fall1
201 real, dimension(num_int_coll_bins) :: fall2,fallr,num,numi
202 real, dimension(n_rhor) :: cgp,crp
203 real, dimension(150) :: mu_r_table
205 real, parameter :: Dm_max = 40000.e-6 ! max. mean ice [m] size for lambda limiter
206 !real, parameter :: Dm_max = 2000.e-6 ! max. mean ice [m] size for lambda limiter
207 real, parameter :: Dm_min = 2.e-6 ! min. mean ice [m] size for lambda limiter
209 real, parameter :: thrd = 1./3.
210 real, parameter :: sxth = 1./6.
211 real, parameter :: cutoff = 1.e-90
213 character(len=2014) :: filename
215 !=== end of variable declaration ===
217 if (log_3momI) then
218 ! HM, no longer need iteration loop
219 ! n_iter_psdSolve = 3 ! 3 iterations found to be sufficient (trial-and-error)
220 n_iter_psdSolve = 1
221 else
222 n_iter_psdSolve = 1
223 endif
226 ! RUNNING IN PARALLEL MODE:
227 !
228 !------------------------------------------------------------------------------------
229 ! CODE ABOVE HERE IS FOR THE "TOP" OF THE BROKEN UP CODE (for running in parallel)
230 !
231 ! Before running ./go_1-compile.ksh, delete all lines below this point and
232 ! and save as 'create_p3_lookupTable_1-top.f90'
233 !------------------------------------------------------------------------------------
235 ! For testing single values, uncomment the following:
236 ! i_Znorm = 1
237 ! i_rhor = 1
239 !------------------------------------------------------------------------------------
240 ! CODE BELOW HERE IS FOR THE "BOTTOM" OF THE BROKEN UP CODE (for running in parallel)
241 !
242 ! Before running ./go_1-compile.ksh, delete all lines below this point and
243 ! and save as 'create_p3_lookupTable_1-bottom.f90'
244 !------------------------------------------------------------------------------------
246 ! set constants and parameters
248 ! assume 600 hPa, 253 K for p and T for fallspeed calcs (for reference air density)
249 pi = acos(-1.)
250 g = 9.861 ! gravity
251 p = 60000. ! air pressure (pa)
252 t = 253.15 ! temp (K)
253 rho = p/(287.15*t) ! air density (kg m-3)
254 mu = 1.496E-6*t**1.5/(t+120.)/rho ! viscosity of air
255 dv = 8.794E-5*t**1.81/p ! diffusivity of water vapor in air
256 dt = 10. ! time step for collection (s)
258 ! parameters for surface roughness of ice particle used for fallspeed
259 ! see mitchell and heymsfield 2005
260 del0 = 5.83
261 c0 = 0.6
262 c1 = 4./(del0**2*c0**0.5)
263 c2 = del0**2/4.
265 dd = 2.e-6 ! bin width for numerical integration of ice processes (units of m)
266 ddd = 50.e-6 ! bin width for numerical integration for ice-ice and ice-rain collection (units of m)
268 !--- specified mass-dimension relationship (cgs units) for unrimed crystals:
270 ! ms = cs*D^ds
271 !
272 ! for graupel:
273 ! mg = cg*D^dg no longer used, since bulk volume is predicted
274 !===
276 !---- Choice of m-D parameters for large unrimed ice:
278 ! Heymsfield et al. 2006
279 ! ds=1.75
280 ! cs=0.0040157+6.06e-5*(-20.)
282 ! sector-like branches (P1b)
283 ! ds=2.02
284 ! cs=0.00142
286 ! bullet-rosette
287 ! ds=2.45
288 ! cs=0.00739
290 ! side planes
291 ! ds=2.3
292 ! cs=0.00419
294 ! radiating assemblages of plates (mitchell et al. 1990)
295 ! ds=2.1
296 ! cs=0.00239
298 ! aggreagtes of side planes, bullets, etc. (Mitchell 1996)
299 ! ds=2.1
300 ! cs=0.0028
302 !-- ** note: if using one of the above (.i.e. not brown and francis, which is already in mks units),
303 ! then uncomment line below to convert from cgs units to mks
304 ! cs=cs*100.**ds/1000.
305 !==
307 ! Brown and Francis (1995)
308 ds = 1.9
309 !cs = 0.01855 ! original (pre v2.3), based on assumption of Dmax
310 cs = 0.0121 ! scaled value based on assumtion of Dmean from Hogan et al. 2012, JAMC
312 !====
314 ! specify m-D parameter for fully rimed ice
315 ! note: cg is not constant, due to variable density
316 dg = 3.
319 !--- projected area-diam relationship (mks units) for unrimed crystals:
320 ! note: projected area = aas*D^bas
322 ! sector-like branches (P1b)
323 ! bas = 1.97
324 ! aas = 0.55*100.**bas/(100.**2)
326 ! bullet-rosettes
327 ! bas = 1.57
328 ! aas = 0.0869*100.**bas/(100.**2)
330 ! aggreagtes of side planes, bullets, etc.
331 bas = 1.88
332 aas = 0.2285*100.**bas/(100.**2)
334 !===
336 !--- projected area-diam relationship (mks units) for fully rimed ice:
337 ! note: projected area = aag*D^bag
339 ! assumed non-spherical
340 ! bag = 2.0
341 ! aag = 0.625*100.**bag/(100.**2)
343 ! assumed spherical:
344 bag = 2.
345 aag = pi*0.25
346 !===
348 dcrit = (pi/(6.*cs)*900.)**(1./(ds-3.))
350 ! check to make sure projected area at dcrit not greater than than of solid sphere
351 ! stop and give warning message if that is the case
353 if (pi/4.*dcrit**2.lt.aas*dcrit**bas) then
354 print*,'STOP, area > area of solid ice sphere, unrimed'
355 stop
356 endif
358 !.........................................................
359 ! generate lookup table for mu (for rain)
360 !
361 ! space of a scaled q/N -- initlamr
363 !Compute mu_r using diagnostic relation:
364 ! ! do i = 1,150 ! loop over lookup table values
365 ! ! initlamr = (real(i)*2.)*1.e-6 + 250.e-6
366 ! ! initlamr = 1./initlamr
367 ! ! ! iterate to get mu_r
368 ! ! ! mu_r-lambda relationship is from Cao et al. (2008), eq. (7)
369 ! ! mu_r = 0. ! first guess
370 ! ! do ii = 1,50
371 ! ! lamr = initlamr*(gamma(mu_r+4.)/(6.*gamma(mu_r+1.)))**thrd
372 ! ! ! new estimate for mu_r based on lambda:
373 ! ! ! set max lambda in formula for mu to 20 mm-1, so Cao et al.
374 ! ! ! formula is not extrapolated beyond Cao et al. data range
375 ! ! dum = min(20.,lamr*1.e-3)
376 ! ! mu_r = max(0.,-0.0201*dum**2+0.902*dum-1.718)
377 ! ! ! if lambda is converged within 0.1%, then exit loop
378 ! ! if (ii.ge.2) then
379 ! ! if (abs((lamold-lamr)/lamr).lt.0.001) goto 111
380 ! ! endif
381 ! ! lamold = lamr
382 ! ! enddo !ii-loop
383 ! ! 111 continue
384 ! ! mu_r_table(i) = mu_r
385 ! ! enddo !i-loop
387 !Precribe a constant mu_r:
388 mu_r_table(:) = 0.
390 !.........................................................
392 ! 3-moment-ice only:
393 ! compute Z value from input Z index whose value is "passed in" through the script
394 ! Z_value = 2.1**(i_Znorm)*1.e-23 ! range from 2x10^(-23) to 600 using 80 values
396 ! alpha parameter of m-D for rimed ice
397 crp(1) = 50.*pi*sxth
398 crp(2) = 250.*pi*sxth
399 crp(3) = 450.*pi*sxth
400 crp(4) = 650.*pi*sxth
401 crp(5) = 900.*pi*sxth
403 !------------------------------------------------------------------------
405 ! open file to write to lookup table:
406 if (log_3momI) then
407 ! write (filename, "(A12,I0.2,A1,I0.2,A4)") "lookupTable_1-",i_Znorm,"_",i_rhor,".dat" !if parallelized over both i_Znorm and i_rhor
408 write (filename, "(A12,I0.2,A4)") "lookupTable_1-LT1.dat"
409 filename = trim(filename)
410 open(unit=1, file=filename, status='unknown')
411 else
412 write (filename, "(A12,I0.2,A4)") "lookupTable_1-",i_rhor,".dat"
413 filename = trim(filename)
414 open(unit=1, file=filename, status='unknown')
415 endif
417 !--
418 ! The values of i_Znorm (and possibly i_rhor) are "passed in" for parallelized version of code for 3-moment.
419 ! The values of i_rhor are "passed in" for parallelized version of code for 2-moment.
420 ! Thus, the loops 'i_Znorm_loop' and 'i_rhor_loop' are commented out accordingingly.
422 if (log_3momI) then
423 Z_value =2.*(i_Znorm-1) !mu_i values set to 0., 2., 4., ... 20.
424 else
425 Z_value = -99. !this gets overwritten below
426 endif
428 !i_Znorm_loop: do i_Znorm = 1,n_Znorm !normally commented (kept to illustrate the structure (and to run in serial)
429 ! i_rhor_loop: do i_rhor = 1,n_rhor !COMMENT OUT FOR PARALLELIZATION (2-MOMENT ONLY)
430 i_Fr_loop_1: do i_Fr = 1,n_Fr !COMMENT OUT FOR PARALLELIZATION (2-MOMENT ONLY)
432 ! write header to first file:
433 if (log_3momI .and. i_Znorm==1 .and. i_rhor==1 .and. i_Fr==1) then
434 write(1,*) 'LOOKUP_TABLE_1-version: ',trim(version),'-3momI'
435 write(1,*)
436 ! elseif (i_rhor==1) then
437 elseif (.not.log_3momI .and. i_rhor==1 .and. i_Fr==1) then
438 write(1,*) 'LOOKUP_TABLE_1-version: ',trim(version),'-2momI'
439 write(1,*)
440 endif
442 !-- these lines to be replaced by Fr(i_Fr) initialization outside of loops
443 ! OR: replace with: Fr = 1./float(n_Fr-1)
444 if (i_Fr.eq.1) Fr = 0.
445 if (i_Fr.eq.2) Fr = 0.333
446 if (i_Fr.eq.3) Fr = 0.667
447 if (i_Fr.eq.4) Fr = 1.
448 !==
450 ! calculate mass-dimension relationship for partially-rimed crystals
451 ! msr = csr*D^dsr
452 ! formula from P3 Part 1 (JAS)
454 ! dcritr is critical size separating fully-rimed from partially-rime ice
456 cgp(i_rhor) = crp(i_rhor) ! first guess
458 if (i_Fr.eq.1) then ! case of no riming (Fr = 0), then we need to set dcrits and dcritr to arbitrary large values
460 dcrits = 1.e+6
461 dcritr = dcrits
462 csr = cs
463 dsr = ds
465 elseif (i_Fr.eq.2.or.i_Fr.eq.3) then ! case of partial riming (Fr between 0 and 1)
467 do
468 dcrits = (cs/cgp(i_rhor))**(1./(dg-ds))
469 dcritr = ((1.+Fr/(1.-Fr))*cs/cgp(i_rhor))**(1./(dg-ds))
470 csr = cs*(1.+Fr/(1.-Fr))
471 dsr = ds
472 ! get mean density of vapor deposition/aggregation grown ice
473 rhodep = 1./(dcritr-dcrits)*6.*cs/(pi*(ds-2.))*(dcritr**(ds-2.)-dcrits**(ds-2.))
474 ! get density of fully-rimed ice as rime mass fraction weighted rime density plus
475 ! density of vapor deposition/aggregation grown ice
476 cgpold = cgp(i_rhor)
477 cgp(i_rhor) = crp(i_rhor)*Fr+rhodep*(1.-Fr)*pi*sxth
478 if (abs((cgp(i_rhor)-cgpold)/cgp(i_rhor)).lt.0.01) goto 115
479 enddo
480 115 continue
482 else ! case of complete riming (Fr=1.0)
484 ! set threshold size between partially-rimed and fully-rimed ice as arbitrary large
485 dcrits = (cs/cgp(i_rhor))**(1./(dg-ds))
486 dcritr = 1.e+6 ! here is the "arbitrary large"
487 csr = cgp(i_rhor)
488 dsr = dg
490 endif
492 !---------------------------------------------------------------------------------------
493 ! set up particle fallspeed arrays
494 ! fallspeed is a function of mass dimension and projected area dimension relationships
495 ! following mitchell and heymsfield (2005), jas
497 ! set up array of particle fallspeed to make computationally efficient
499 ! for high-resolution (in diameter space), ice fallspeed is stored in 'fall1' array (m/s)
500 ! for lower-resolution (in diameter space), ice fallspeed is stored in 'fall2' array (m/s)
501 ! rain fallspeed is stored in 'fallr' (m/s)
503 ! loop over particle size
505 jj_loop_1: do jj = 1,num_int_bins
507 ! particle size (m)
508 d1 = real(jj)*dd - 0.5*dd
510 !----- get mass-size and projected area-size relationships for given size (d1)
511 ! call get_mass_size
513 if (d1.le.dcrit) then
514 cs1 = pi*sxth*900.
515 ds1 = 3.
516 bas1 = 2.
517 aas1 = pi/4.
518 else if (d1.gt.dcrit.and.d1.le.dcrits) then
519 cs1 = cs
520 ds1 = ds
521 bas1 = bas
522 aas1 = aas
523 else if (d1.gt.dcrits.and.d1.le.dcritr) then
524 cs1 = cgp(i_rhor)
525 ds1 = dg
526 bas1 = bag
527 aas1 = aag
528 else if (d1.gt.dcritr) then
529 cs1 = csr
530 ds1 = dsr
531 if (i_Fr.eq.1) then
532 aas1 = aas
533 bas1 = bas
534 else
535 ! for projected area, keep bas1 constant, but modify aas1 according to rimed fraction
536 bas1 = bas
537 dum1 = aas*d1**bas
538 dum2 = aag*d1**bag
539 m1 = cs1*d1**ds1
540 m2 = cs*d1**ds
541 m3 = cgp(i_rhor)*d1**dg
542 ! linearly interpolate based on particle mass
543 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
544 ! dum3 = (1.-Fr)*dum1+Fr*dum2 !DELETE?
545 aas1 = dum3/(d1**bas)
546 endif
547 endif
548 !=====
550 ! correction for turbulence
551 ! if (d1.lt.500.e-6) then
552 ! a0 = 0.
553 ! b0 = 0.
554 ! else
555 ! a0=1.7e-3
556 ! b0=0.8
557 ! endif
558 ! neglect turbulent correction for aggregates for now (no condition)
559 a0 = 0.
560 b0 = 0.
562 ! fall speed for particle
563 ! Best number:
564 xx = 2.*cs1*g*rho*d1**(ds1+2.-bas1)/(aas1*(mu*rho)**2)
565 ! drag terms:
566 b1 = c1*xx**0.5/(2.*((1.+c1*xx**0.5)**0.5-1.)*(1.+c1*xx**0.5)**0.5)-a0*b0*xx** &
567 b0/(c2*((1.+c1*xx**0.5)**0.5-1.)**2)
568 a1 = (c2*((1.+c1*xx**0.5)**0.5-1.)**2-a0*xx**b0)/xx**b1
569 ! velocity in terms of drag terms
570 fall1(jj) = a1*mu**(1.-2.*b1)*(2.*cs1*g/(rho*aas1))**b1*d1**(b1*(ds1-bas1+2.)-1.)
572 enddo jj_loop_1
574 !................................................................
575 ! fallspeed array for ice-ice and ice-rain collision calculations
577 jj_loop_2: do jj = 1,num_int_coll_bins
579 ! particle size:
580 d1 = real(jj)*ddd - 0.5*ddd
582 if (d1.le.dcrit) then
583 cs1 = pi*sxth*900.
584 ds1 = 3.
585 bas1 = 2.
586 aas1 = pi/4.
587 else if (d1.gt.dcrit.and.d1.le.dcrits) then
588 cs1 = cs
589 ds1 = ds
590 bas1 = bas
591 aas1 = aas
592 else if (d1.gt.dcrits.and.d1.le.dcritr) then
593 cs1 = cgp(i_rhor)
594 ds1 = dg
595 bas1 = bag
596 aas1 = aag
597 else if (d1.gt.dcritr) then
598 cs1 = csr
599 ds1 = dsr
600 if (i_Fr.eq.1) then
601 aas1 = aas
602 bas1 = bas
603 else
604 ! for area, keep bas1 constant, but modify aas1 according to rimed fraction
605 bas1 = bas
606 dum1 = aas*d1**bas
607 dum2 = aag*d1**bag
608 ! dum3 = (1.-Fr)*dum1+Fr*dum2
609 m1 = cs1*d1**ds1
610 m2 = cs*d1**ds
611 m3 = cgp(i_rhor)*d1**dg
612 ! linearly interpolate based on particle mass:
613 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
614 aas1 = dum3/(d1**bas)
615 endif
616 endif
618 ! correction for turbulence
619 ! if (d1.lt.500.e-6) then
620 ! a0 = 0.
621 ! b0 = 0.
622 ! else
623 ! a0=1.7e-3
624 ! b0=0.8
625 ! endif
626 ! neglect turbulent correction for aggregates for now (no condition)
627 a0 = 0.
628 b0 = 0.
630 ! fall speed for ice
631 ! Best number:
632 xx = 2.*cs1*g*rho*d1**(ds1+2.-bas1)/(aas1*(mu*rho)**2)
633 ! drag terms
634 b1 = c1*xx**0.5/(2.*((1.+c1*xx**0.5)**0.5-1.)*(1.+c1*xx**0.5)**0.5)-a0*b0*xx** &
635 b0/(c2*((1.+c1*xx**0.5)**0.5-1.)**2)
636 a1 = (c2*((1.+c1*xx**0.5)**0.5-1.)**2-a0*xx**b0)/xx**b1
637 ! velocity in terms of drag terms
638 fall2(jj) = a1*mu**(1.-2.*b1)*(2.*cs1*g/(rho*aas1))**b1*d1**(b1*(ds1-bas1+2.)-1.)
640 !------------------------------------
641 ! fall speed for rain particle
643 dia = d1 ! diameter m
644 amg = pi*sxth*997.*dia**3 ! mass [kg]
645 amg = amg*1000. ! convert kg to g
647 if (dia.le.134.43e-6) then
648 dum2 = 4.5795e5*amg**(2.*thrd)
649 goto 101
650 endif
652 if(dia.lt.1511.64e-6) then
653 dum2 = 4.962e3*amg**thrd
654 goto 101
655 endif
657 if(dia.lt.3477.84e-6) then
658 dum2 = 1.732e3*amg**sxth
659 goto 101
660 endif
662 dum2 = 917.
664 101 continue
666 fallr(jj) = dum2*1.e-2 ! convert (cm s-1) to (m s-1)
668 enddo jj_loop_2
670 !---------------------------------------------------------------------------------
672 ! loop around normalized Q (Qnorm)
673 i_Qnorm_loop: do i_Qnorm = 1,n_Qnorm
675 ! lookup table values of normalized Qnorm
676 ! (range of mean mass diameter from ~ 1 micron to x cm)
678 !q = 261.7**((i_Qnorm+10)*0.1)*1.e-18 ! old (strict) lambda limiter
679 q = 800.**((i_Qnorm+10)*0.1)*1.e-18 ! new lambda limiter
681 !--uncomment to test and print proposed values of qovn
682 ! print*,i_Qnorm,(6./(pi*500.)*q)**0.3333
683 ! enddo
684 ! stop
685 !==
687 ! test values
688 ! N = 5.e+3
689 ! q = 0.01e-3
691 ! initialize qerror to arbitrarily large value:
692 qerror = 1.e+20
694 !.....................................................................................
695 ! Find parameters for gamma distribution
697 ! size distribution for ice is assumed to be
698 ! N(D) = n0 * D^mu_i * exp(-lam*D)
700 ! for the given q and N, we need to find n0, mu_i, and lam
702 ! approach for finding lambda:
703 ! cycle through a range of lambda, find closest lambda to produce correct q
705 ! start with lam, range of lam from 100 to 1 x 10^7 is large enough to
706 ! cover full range over mean size from approximately 1 micron to x cm
708 ! HM, no longer needed
709 ! if (log_3momI) then
710 ! ! assign provisional values for mom3 (first guess for mom3)
711 ! ! NOTE: these are normalized: mom3 = M3/M0, mom6 = M6/M3 (M3 = 3rd moment, etc.)
712 ! mom3 = q/cgp(i_rhor) !note: cgp is pi/6*(mean_density), computed above
713 ! ! update normalized mom6 based on the updated ratio of normalized mom3 and normalized Q
714 ! ! (we want mom6 normalized by mom3 not q)
715 ! dum = mom3/q
716 ! mom6 = Z_value/dum
717 ! endif !log_3momI
718 !==
720 iteration_loop1: do i_iter = 1,n_iter_psdSolve
722 if (log_3momI) then
723 ! compute mu_i from normalized mom3 and mom6:
724 ! HM set to loop value of mu (temporarily called Z_value)
725 ! mu_i = compute_mu_3moment(mom3,mom6,mu_i_max)
726 ! mu_i = max(mu_i,mu_i_min) ! make sure mu_i >= 0 (otherwise size dist is infinity at D = 0)
727 ! mu_i = min(mu_i,mu_i_max) ! set upper limit
728 mu_i = Z_value
729 endif
731 ii_loop_1: do ii = 1,11000 ! this range of ii to calculate lambda chosen by trial and error for the given lambda limiter values
733 ! lam = 1.0013**ii*100. ! old (strict) lambda_i limiter
734 lam = 1.0013**ii*10. ! new lambda_i limiter
736 ! solve for mu_i for 2-moment-ice:
737 if (.not. log_3momI) mu_i = diagnostic_mui(mu_i_min,mu_i_max,lam,q,cgp(i_rhor),Fr,pi)
739 ! for lambda limiter:
740 !dum = Dm_max+Fr*(3000.e-6)
741 dum = Dm_max
742 lam = max(lam,(mu_i+1.)/dum) ! set min lam corresponding to mean size of x
743 lam = min(lam,(mu_i+1.)/Dm_min) ! set max lam corresponding to mean size of Dm_min (2 micron)
745 ! normalized n0:
746 n0 = lam**(mu_i+1.)/(gamma(mu_i+1.))
748 ! calculate integral for each of the 4 parts of the size distribution
749 ! check difference with respect to Qnorm
751 ! set up m-D relationship for solid ice with D < Dcrit:
752 cs1 = pi*sxth*900.
753 ds1 = 3.
755 call intgrl_section(lam,mu_i, ds1,ds,dg,dsr, dcrit,dcrits,dcritr,intgrR1,intgrR2,intgrR3,intgrR4)
756 ! intgrR1 is integral from 0 to dcrit (solid ice)
757 ! intgrR2 is integral from dcrit to dcrits (unrimed large ice)
758 ! intgrR3 is integral from dcrits to dcritr (fully rimed ice)
759 ! intgrR4 is integral from dcritr to inf (partially rimed)
761 ! sum of the integrals from the 4 regions of the size distribution:
762 qdum = n0*(cs1*intgrR1 + cs*intgrR2 + cgp(i_rhor)*intgrR3 + csr*intgrR4)
764 !--- numerical integration for test to make sure incomplete gamma function is working
765 ! sum1 = 0.
766 ! dd = 1.e-6
767 ! do iii=1,50000
768 ! dum=real(iii)*dd
769 ! if (dum.lt.dcrit) then
770 ! sum1 = sum1+n0*dum**mu_i*cs1*dum**ds1*exp(-lam*dum)*dd
771 ! else
772 ! sum1 = sum1+n0*dum**mu_i*cs*dum**ds*exp(-lam*dum)*dd
773 ! end if
774 ! enddo
775 ! print*,'sum1=',sum1
776 ! stop
777 !===
779 if (ii.eq.1) then
780 qerror = abs(q-qdum)
781 lamf = lam
782 endif
784 ! find lam with smallest difference between Qnorm and estimate of Qnorm, assign to lamf
785 if (abs(q-qdum).lt.qerror) then
786 lamf = lam
787 qerror = abs(q-qdum)
788 endif
790 enddo ii_loop_1
792 ! check and print relative error in q to make sure it is not too large
793 ! note: large error is possible if size bounds are exceeded!!!!!!!!!!
794 ! print*,'qerror (%)',qerror/q*100.
796 ! find n0 based on final lam value
797 ! set final lamf to 'lam' variable
798 ! this is the value of lam with the smallest qerror
799 lam = lamf
801 ! recalculate mu_i based on final lam (for 2-moment-ice only; not needed for 3-moment-ice)
802 if (.not. log_3momI) mu_i = diagnostic_mui(mu_i_min,mu_i_max,lam,q,cgp(i_rhor),Fr,pi)
804 ! n0 = N*lam**(mu_i+1.)/(gamma(mu_i+1.))
806 ! find n0 from lam and Qnorm:
807 ! (this is done instead of finding n0 from lam and N, since N;
808 ! may need to be adjusted to constrain mean size within reasonable bounds)
810 call intgrl_section(lam,mu_i, ds1,ds,dg,dsr, dcrit,dcrits,dcritr,intgrR1,intgrR2,intgrR3,intgrR4)
811 n0 = q/(cs1*intgrR1 + cs*intgrR2 + cgp(i_rhor)*intgrR3 + csr*intgrR4)
813 ! print*,'lam,N0,mu:',lam,n0,mu_i
815 !--- test final lam, N0 values
816 ! sum1 = 0.
817 ! dd = 1.e-6
818 ! do iii=1,50000
819 ! dum=real(iii)*dd
820 ! if (dum.lt.dcrit) then
821 ! sum1 = sum1+n0*dum**mu_i*cs1*dum**ds1*exp(-lam*dum)*dd
822 ! elseif (dum.ge.dcrit.and.dum.lt.dcrits) then
823 ! sum1 = sum1+n0*dum**mu_i*cs*dum**ds*exp(-lam*dum)*dd
824 ! elseif (dum.ge.dcrits.and.dum.lt.dcritr) then
825 ! sum1 = sum1+n0*dum**mu_i*cg*dum**dg*exp(-lam*dum)*dd
826 ! elseif (dum.ge.dcritr) then
827 ! sum1 = sum1+n0*dum**mu_i*csr*dum**dsr*exp(-lam*dum)*dd
828 ! endif
829 ! enddo
830 ! print*,'sum1=',sum1
831 ! stop
832 !===
834 ! calculate normalized mom3 directly from PSD parameters (3-moment-ice only)
835 ! HM no longer needed
836 ! if (log_3momI) then
837 ! mom3 = n0*gamma(4.+mu_i)/lam**(4.+mu_i)
838 ! ! update normalized mom6 based on the updated ratio of normalized mom3 and normalized Q
839 ! ! (we want mom6 normalized by mom3 not q)
840 ! dum = mom3/q
841 ! mom6 = Z_value/dum
842 ! endif !log_3momI
844 enddo iteration_loop1
846 lambda_i(i_Qnorm,i_Fr) = lam
847 mu_i_save(i_Qnorm,i_Fr) = mu_i
849 !.....................................................................................
850 ! At this point, we have solved for all of the ice size distribution parameters (n0, lam, mu_i)
851 !.....................................................................................
853 !.....................................................................................
854 ! find max/min Q* to constrain mean size (i.e. lambda limiter), this is stored and passed to
855 ! lookup table, so that N (nitot) can be adjusted during the simulation to constrain mean size
856 ! (computed and written as the inverses (i_qsmall,i_qlarge) to avoid run-time division in p3_main)
858 ! limit based on min size, Dm_min (2 micron):
859 duml = (mu_i+1.)/Dm_min
860 call intgrl_section(duml,mu_i, ds1,ds,dg,dsr, dcrit,dcrits,dcritr,intgrR1,intgrR2,intgrR3,intgrR4)
861 n0dum = q/(cs1*intgrR1 + cs*intgrR2 + cgp(i_rhor)*intgrR3 + csr*intgrR4)
863 !find maximum N applying the lambda limiter (lower size limit)
864 dum =n0dum/(duml**(mu_i+1.)/(gamma(mu_i+1.)))
866 !calculate the lower limit of normalized Q to use in P3 main
867 !(this is based on the lower limit of mean size so we call this 'qsmall')
868 !qsmall(i_Qnorm,i_Fr) = q/dum
869 i_qsmall(i_Qnorm,i_Fr) = dum/q
872 ! limit based on max size, Dm_max:
873 duml = (mu_i+1.)/Dm_max
875 call intgrl_section(duml,mu_i, ds1,ds,dg,dsr, dcrit,dcrits,dcritr,intgrR1,intgrR2,intgrR3,intgrR4)
876 n0dum = q/(cs1*intgrR1 + cs*intgrR2 + cgp(i_rhor)*intgrR3 + csr*intgrR4)
878 ! find minium N applying the lambda limiter (lower size limit)
879 dum = n0dum/(duml**(mu_i+1.)/(gamma(mu_i+1.)))
881 ! calculate the upper limit of normalized Q to use in P3 main
882 ! (this is based on the upper limit of mean size so we call this 'qlarge')
883 !qlarge(i_Qnorm,i_Fr) = q/dum
884 i_qlarge(i_Qnorm,i_Fr) = dum/q
886 ! calculate bounds for normalized Z based on min/max allowed mu: (3-moment-ice only)
887 ! HM no longer needed, don't need to calculate or output zlarge and zsmall
888 ! if (log_3momI) then
889 ! mu_dum = mu_i_min
890 ! gdum = (6.+mu_dum)*(5.+mu_dum)*(4.+mu_dum)/((3.+mu_dum)*(2.+mu_dum)*(1.+mu_dum))
891 ! dum = mom3/q
892 ! zlarge(i_Qnorm,i_Fr) = gdum*mom3*dum
893 ! mu_dum = mu_i_max
894 ! gdum = (6.+mu_dum)*(5.+mu_dum)*(4.+mu_dum)/((3.+mu_dum)*(2.+mu_dum)*(1.+mu_dum))
895 ! zsmall(i_Qnorm,i_Fr) = gdum*mom3*dum
896 ! endif !if (log_3momI)
898 !.....................................................................................
899 ! begin moment and microphysical process calculations for the lookup table
901 !.....................................................................................
902 ! mass- and number-weighted mean fallspeed (m/s)
903 ! add reflectivity
904 !.....................................................................................
906 ! assume conditions for t and p as assumed above (giving rhos), then in microphysics scheme
907 ! multiply by density correction factor (rhos/rho)^0.54, from Heymsfield et al. 2006
909 ! fallspeed formulation from Mitchell and Heymsfield 2005
911 ! initialize for numerical integration
912 sum1 = 0.
913 sum2 = 0.
914 sum3 = 0.
915 sum4 = 0.
916 sum5 = 0. ! reflectivity
917 sum6 = 0. ! mass mean size
918 sum7 = 0. ! mass-weighted mean density
919 sum8 = 0. ! 6th moment * velocity [3momI only]
920 sum9 = 0. ! 6th moment [3momI only]
922 ! numerically integrate over size distribution
923 ii_loop_2: do ii = 1,num_int_bins
925 dum = real(ii)*dd - 0.5*dd ! particle size
927 !assign mass-size parameters (depending on size at ii)
928 if (dum.le.dcrit) then
929 ds1 = 3.
930 cs1 = pi*sxth*900.
931 else if (dum.gt.dcrit.and.dum.le.dcrits) then
932 ds1 = ds
933 cs1 = cs
934 elseif (dum.gt.dcrits.and.dum.le.dcritr) then
935 ds1 = dg
936 cs1 = cgp(i_rhor)
937 elseif (dum.gt.dcritr) then
938 ds1 = dsr
939 cs1 = csr
940 endif
942 ! numerator of number-weighted velocity - sum1:
943 sum1 = sum1+fall1(ii)*dum**mu_i*exp(-lam*dum)*dd
945 ! numerator of mass-weighted velocity - sum2:
946 sum2 = sum2+fall1(ii)*cs1*dum**(ds1+mu_i)*exp(-lam*dum)*dd
948 ! total number and mass for weighting above fallspeeds:
949 ! (note: do not need to include n0 and cs since these parameters are in both numerator and denominator
950 !denominator of number-weighted V:
951 sum3 = sum3+dum**mu_i*exp(-lam*dum)*dd
953 !denominator of mass-weighted V:
954 sum4 = sum4+cs1*dum**(ds1+mu_i)*exp(-lam*dum)*dd
956 !reflectivity (integral of mass moment squared):
957 sum5 = sum5+n0*(cs1/917.)**2*(6./pi)**2*dum**(2.*ds1+mu_i)*exp(-lam*dum)*dd
959 !numerator of mass-weighted mean size
960 sum6 = sum6+cs1*dum**(ds1+mu_i+1.)*exp(-lam*dum)*dd
962 !numerator of mass-weighted density:
963 ! particle density is defined as mass divided by volume of sphere with same D
964 sum7 = sum7+(cs1*dum**ds1)**2/(pi*sxth*dum**3)*dum**mu_i*exp(-lam*dum)*dd
966 !numerator in 6th-moment-weight fall speed [3momI only]
967 sum8 = sum8 + fall1(ii)*dum**(mu_i+6.)*exp(-lam*dum)*dd
969 !denominator in 6th-moment-weight fall speed [3momI only]
970 sum9 = sum9 + dum**(mu_i+6.)*exp(-lam*dum)*dd
972 enddo ii_loop_2
974 ! save mean fallspeeds for lookup table:
975 uns(i_Qnorm,i_Fr) = sum1/sum3
976 ums(i_Qnorm,i_Fr) = sum2/sum4
977 refl(i_Qnorm,i_Fr) = sum5
978 dmm(i_Qnorm,i_Fr) = sum6/sum4
979 rhomm(i_Qnorm,i_Fr) = sum7/sum4
980 if (log_3momI) then
981 uzs(i_Qnorm,i_Fr) = sum8/sum9
982 ! write(6,'(a12,3e15.5)') 'uzs,ums,uns',uzs(i_Qnorm,i_Fr),ums(i_Qnorm,i_Fr),uns(i_Qnorm,i_Fr)
983 endif
985 !.....................................................................................
986 ! self-aggregation
987 !.....................................................................................
989 sum1 = 0.
991 ! set up binned distribution of ice
992 do jj = num_int_coll_bins,1,-1
993 d1 = real(jj)*ddd - 0.5*ddd
994 num(jj) = n0*d1**mu_i*exp(-lam*d1)*ddd
995 enddo !jj-loop
997 ! loop over exponential size distribution
998 ! ! note: collection of ice within the same bin is neglected
1000 jj_loop_3: do jj = num_int_coll_bins,1,-1
1001 kk_loop_1: do kk = 1,jj-1
1003 ! particle size:
1004 d1 = real(jj)*ddd - 0.5*ddd
1005 d2 = real(kk)*ddd - 0.5*ddd
1007 if (d1.le.dcrit) then
1008 bas1 = 2.
1009 aas1 = pi*0.25
1010 elseif (d1.gt.dcrit.and.d1.le.dcrits) then
1011 bas1 = bas
1012 aas1 = aas
1013 else if (d1.gt.dcrits.and.d1.le.dcritr) then
1014 bas1 = bag
1015 aas1 = aag
1016 else if (d1.gt.dcritr) then
1017 cs1 = csr
1018 ds1 = dsr
1019 if (i_Fr.eq.1) then
1020 aas1 = aas
1021 bas1 = bas
1022 else
1023 ! for area, keep bas1 constant, but modify aas1 according to rimed fraction
1024 bas1 = bas
1025 dum1 = aas*d1**bas
1026 dum2 = aag*d1**bag
1027 m1 = cs1*d1**ds1
1028 m2 = cs*d1**ds
1029 m3 = cgp(i_rhor)*d1**dg
1030 ! linearly interpolate based on particle mass
1031 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
1032 aas1 = dum3/(d1**bas)
1033 endif
1034 endif
1036 ! parameters for particle 2
1037 if (d2.le.dcrit) then
1038 bas2 = 2.
1039 aas2 = pi/4.
1040 elseif (d2.gt.dcrit.and.d2.le.dcrits) then
1041 bas2 = bas
1042 aas2 = aas
1043 elseif (d2.gt.dcrits.and.d2.le.dcritr) then
1044 bas2 = bag
1045 aas2 = aag
1046 elseif (d2.gt.dcritr) then
1047 cs2 = csr
1048 ds2 = dsr
1049 if (i_Fr.eq.1) then
1050 aas2 = aas
1051 bas2 = bas
1052 else
1053 ! for area, keep bas1 constant, but modify aas1 according to rime fraction
1054 bas2 = bas
1055 dum1 = aas*d2**bas
1056 dum2 = aag*d2**bag
1057 m1 = cs2*d2**ds2
1058 m2 = cs*d2**ds
1059 m3 = cgp(i_rhor)*d2**dg
1060 ! linearly interpolate based on particle mass
1061 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
1062 aas2 = dum3/(d2**bas)
1063 endif
1064 endif
1066 ! differential fallspeed:
1067 ! (note: in P3_MAIN must multiply by air density correction factor, and collection efficiency
1068 delu = abs(fall2(jj)-fall2(kk))
1070 ! sum for integral
1072 ! sum1 = # of collision pairs
1073 ! the assumption is that each collision pair reduces crystal
1074 ! number mixing ratio by 1 kg^-1 s^-1 per kg/m^3 of air (this is
1075 ! why we need to multiply by air density, to get units of 1/kg^-1 s^-1)
1077 sum1 = sum1+(aas1*d1**bas1+aas2*d2**bas2)*delu*num(jj)*num(kk)
1079 ! remove collected particles from distribution over time period dt, update num
1080 ! note -- dt is time scale for removal, not model time step
1081 ! num(kk) = num(kk)-(aas1*d1**bas1+aas2*d2**bas2)*delu*num(jj)*num(kk)*dt
1082 ! num(kk) = max(num(kk),0.)
1084 ! write(6,'(2i5,8e15.5)')jj,kk,sum1,num(jj),num(kk),delu,aas1,d1,aas2,d2
1085 ! num(kk)=num(kk)-(aas1*d1**bas1+aas2*d2**bas2)*delu*num(jj)*num(kk)*0.1*0.5
1086 ! num(kk)=max(num(kk),0.)
1087 ! sum1 = sum1+0.5*(aas1*d1**bas1+aas2*d2**bas2)*delu*n0*n0*(d1+d2)**mu_i*exp(-lam*(d1+d2))*dd**2
1089 enddo kk_loop_1
1090 enddo jj_loop_3
1092 nagg(i_Qnorm,i_Fr) = sum1 ! save to write to output
1094 ! print*,'nagg',nagg(i_Qnorm,i_Fr)
1096 !.....................................................................................
1097 ! collection of cloud droplets
1098 !.....................................................................................
1099 ! note: In P3_MAIN, needs to be multiplied by collection efficiency Eci
1100 ! Also needs to be multiplied by air density correction factor for fallspeed,
1101 ! ! air density, and cloud water mixing ratio or number concentration
1103 ! initialize sum for integral
1104 sum1 = 0.
1105 sum2 = 0.
1107 ! loop over exponential size distribution (from 1 micron to 2 cm)
1108 jj_loop_4: do jj = 1,num_int_bins
1110 d1 = real(jj)*dd - 0.5*dd ! particle size or dimension (m) for numerical integration
1112 ! get mass-dimension and projected area-dimension relationships
1113 ! for different ice types across the size distribution based on critical dimensions
1114 ! separating these ice types (see Fig. 2, morrison and grabowski 2008)
1116 ! mass = cs1*D^ds1
1117 ! projected area = bas1*D^bas1
1118 if (d1.le.dcrit) then
1119 cs1 = pi*sxth*900.
1120 ds1 = 3.
1121 bas1 = 2.
1122 aas1 = pi*0.25
1123 elseif (d1.gt.dcrit.and.d1.le.dcrits) then
1124 cs1 = cs
1125 ds1 = ds
1126 bas1 = bas
1127 aas1 = aas
1128 elseif (d1.gt.dcrits.and.d1.le.dcritr) then
1129 cs1 = cgp(i_rhor)
1130 ds1 = dg
1131 bas1 = bag
1132 aas1 = aag
1133 elseif (d1.gt.dcritr) then
1134 cs1 = csr
1135 ds1 = dsr
1136 if (i_Fr.eq.1) then
1137 aas1 = aas
1138 bas1 = bas
1139 else
1140 ! for area, ! keep bas1 constant, but modify aas1 according to rimed fraction
1141 bas1 = bas
1142 dum1 = aas*d1**bas
1143 dum2 = aag*d1**bag
1144 m1 = cs1*d1**ds1
1145 m2 = cs*d1**ds
1146 m3 = cgp(i_rhor)*d1**dg
1147 ! linearly interpolate based on particle mass
1148 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
1149 aas1 = dum3/(d1**bas)
1150 endif
1151 endif
1153 ! sum for integral
1154 ! include assumed ice particle size threshold for riming of 100 micron
1155 ! note: sum1 (nrwat) is the scaled collection rate, units of m^3 kg^-1 s^-1
1156 ! sum2 (nrwats) is mass of snow times scaled collection rate, units of m^3 s^-1
1158 if (d1.ge.100.e-6) then
1159 sum1 = sum1+aas1*d1**bas1*fall1(jj)*n0*d1**mu_i*exp(-lam*d1)*dd
1160 !sum2 = sum2+aas1*d1**bas1*fall1(jj)*n0*d1**mu_i*exp(-lam*d1)*dd*cs1*d1**ds1
1161 endif
1163 enddo jj_loop_4
1165 ! save for output
1166 nrwat(i_Qnorm,i_Fr) = sum1 ! note: read in as 'f1pr4' in P3_MAIN
1167 !nrwats(i_Qnorm,i_Fr) = sum2
1169 !.....................................................................................
1170 ! collection of rain
1171 !.....................................................................................
1173 ! note: In P3_MAIN, we need to multiply rate by n0r, collection efficiency,
1174 ! air density, and air density correction factor
1176 ! This approach implicitly assumes that the PSDs are constant during the microphysics
1177 ! time step, this could produce errors if the time step is large. In particular,
1178 ! more mass or number could be removed than is available. This will be taken care
1179 ! of by conservation checks in the microphysics code.
1181 ! loop around lambda for rain
1182 i_Drscale_loop: do i_Drscale = 1,n_Drscale
1184 print*,'** STATUS: ',i_Znorm, i_rhor, i_Fr, i_Qnorm, i_Drscale
1186 dum = 1.24**i_Drscale*10.e-6
1188 ! assumed lamv for tests
1189 ! dum = 7.16e-5
1190 ! note: lookup table for rain is based on lamv, i.e.,inverse volume mean diameter
1191 lamv = 1./dum
1192 lamrs(i_Drscale) = lamv
1194 ! get mu_r from lamr:
1195 dum = 1./lamv
1197 if (dum.lt.282.e-6) then
1198 mu_r = 8.282
1199 elseif (dum.ge.282.e-6 .and. dum.lt.502.e-6) then
1200 ! interpolate:
1201 rdumii = (dum-250.e-6)*1.e6*0.5
1202 rdumii = max(rdumii,1.)
1203 rdumii = min(rdumii,150.)
1204 dumii = int(rdumii)
1205 dumii = min(149,dumii)
1206 mu_r = mu_r_table(dumii)+(mu_r_table(dumii+1)-mu_r_table(dumii))*(rdumii-real(dumii))
1207 elseif (dum.ge.502.e-6) then
1208 mu_r = 0.
1209 endif
1210 ! recalculate slope based on mu_r
1211 !LAMR = (pi*sxth*rhow*nr(i_Qnorm,k)*gamma(mu_r+4.)/(qr(i_Qnorm,k)*gamma(mu_r+1.)))**thrd
1213 ! this is done by re-scaling lamv to account for DSD shape (mu_r)
1214 lamr = lamv*(gamma(mu_r+4.)/(6.*gamma(mu_r+1.)))**thrd
1216 ! set maximum value for rain lambda
1217 !lammax = (mu_r+1.)/10.e-6
1218 lammax = (mu_r+1.)*1.e+5
1220 ! set to small value since breakup is explicitly included (mean size 5 mm)
1221 !lammin = (mu_r+1.)/5000.e-6
1222 lammin = (mu_r+1.)*200.
1223 lamr = min(lamr,lammax)
1224 lamr = max(lamr,lammin)
1226 ! initialize sum
1227 sum1 = 0.
1228 sum2 = 0.
1229 sum6 = 0.
1230 ! sum8 = 0. ! total rain
1232 do jj = 1,num_int_coll_bins
1233 ! particle size:
1234 d1 = real(jj)*ddd - 0.5*ddd
1235 ! num is the scaled binned rain size distribution;
1236 ! need to multiply by n0r to get unscaled distribution
1237 num(jj) = d1**mu_r*exp(-lamr*d1)*ddd
1238 ! get (unscaled) binned ice size distribution
1239 numi(jj) = n0*d1**mu_i*exp(-lam*d1)*ddd
1240 enddo !jj-loop
1242 ! loop over rain and ice size distributions
1243 jj_loop_5: do jj = 1,num_int_coll_bins
1244 kk_loop_2: do kk = 1,num_int_coll_bins
1246 ! particle size:
1247 d1 = real(jj)*ddd - 0.5*ddd ! ice
1248 d2 = real(kk)*ddd - 0.5*ddd ! rain
1249 if (d1.le.dcrit) then
1250 cs1 = pi*sxth*900.
1251 ds1 = 3.
1252 bas1 = 2.
1253 aas1 = pi*0.25
1254 elseif (d1.gt.dcrit.and.d1.le.dcrits) then
1255 cs1 = cs
1256 ds1 = ds
1257 bas1 = bas
1258 aas1 = aas
1259 elseif (d1.gt.dcrits.and.d1.le.dcritr) then
1260 cs1 = cgp(i_rhor)
1261 ds1 = dg
1262 bas1 = bag
1263 aas1 = aag
1264 else if (d1.gt.dcritr) then
1265 cs1 = csr
1266 ds1 = dsr
1267 if (i_Fr.eq.1) then
1268 aas1 = aas
1269 bas1 = bas
1270 else
1271 ! for area, keep bas1 constant, but modify aas1 according to rime fraction
1272 bas1 = bas
1273 dum1 = aas*d1**bas
1274 dum2 = aag*d1**bag
1275 m1 = cs1*d1**ds1
1276 m2 = cs*d1**ds
1277 m3 = cgp(i_rhor)*d1**dg
1278 ! linearly interpolate based on particle mass
1279 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
1280 aas1 = dum3/(d1**bas)
1281 endif
1282 endif
1284 delu = abs(fall2(jj)-fallr(kk)) ! differential fallspeed
1286 !......................................................
1287 ! collection of rain mass and number
1289 ! allow collection of rain both when rain fallspeed > ice fallspeed and ice fallspeed > rain fallspeed
1290 ! this is applied below freezing to calculate total amount of rain mass and number that collides with ice and freezes
1292 ! if (fall2(jj).ge.fallr(kk)) then
1294 ! sum for integral:
1296 ! change in rain N (units of m^4 s^-1 kg^-1), thus need to multiply
1297 ! by air density (units kg m^-3) and n0r (units kg^-1 m^-1) in P3_MAIN
1299 !sum1 = sum1+(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* &
1300 ! exp(-lam*d1)* &dd*num(kk)
1301 sum1 = sum1+(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* &
1302 exp(-lam*d1)*ddd*num(kk)
1303 !sum1 = sum1+min((aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* &
1304 ! exp(-lam*d1)*dd*num(kk),num(kk))
1306 ! change in rain q (units of m^4 s^-1), again need to multiply by air density and n0r in P3_MAIN
1308 !sum2 = sum2+(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* &
1309 ! exp(-lam*d1)*dd*num(kk)*pi*sxth*997.*d2**3
1310 sum2 = sum2+(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* &
1311 exp(-lam*d1)*ddd*num(kk)*pi*sxth*997.*d2**3
1313 ! remove collected rain drops from distribution:
1314 !num(kk) = num(kk)-(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* &
1315 ! exp(-lam*d1)*dd*num(kk)*dt
1316 !num(kk) = max(num(kk),0.)
1318 !......................................................
1319 ! now calculate collection of ice mass by rain
1321 ! ice collecting rain
1323 ! again, allow collection both when ice fallspeed > rain fallspeed
1324 ! and when rain fallspeed > ice fallspeed
1325 ! this is applied to conditions above freezing to calculate
1326 ! acceleration of melting due to collisions with liquid (rain)
1328 ! if (fall2(jj).ge.fallr(kk)) then
1330 ! collection of ice number
1332 ! sum5 = sum5+(aas1*d1**bas1+pi/4.*d2**2)*delu*exp(-lamr*d2)*dd*numi(jj)
1334 ! collection of ice mass (units of m^4 s^-1)
1335 ! note: need to multiply by air density and n0r in microphysics code
1336 sum6 = sum6+(aas1*d1**bas1+pi*0.25*d2**2)*delu*d2**mu_r* &
1337 exp(-lamr*d2)*ddd*numi(jj)*cs1*d1**ds1
1339 ! remove collected snow from distribution:
1340 !numi(jj) = numi(jj)-(aas1*d1**bas1+pi*0.25*d2**2)*delu*d2**mu_r* &
1341 ! exp(-lamr*d2)*dd*numi(jj)*dt
1342 !numi(jj) = max(numi(jj),0.)
1344 enddo kk_loop_2
1345 enddo jj_loop_5
1347 ! save for output:
1348 nrrain(i_Qnorm,i_Drscale,i_Fr) = sum1
1349 qrrain(i_Qnorm,i_Drscale,i_Fr) = sum2
1350 qsrain(i_Qnorm,i_Drscale,i_Fr) = sum6
1352 enddo i_Drscale_loop !(loop around lambda for rain)
1354 !.....................................................................................
1355 ! vapor deposition/melting/wet growth
1356 !.....................................................................................
1358 ! vapor deposition including ventilation effects
1359 ! note: in microphysics code we need to multiply by air density and
1360 ! (mu/dv)^0.3333*(rhofac/mu)^0.5, where rhofac is air density correction factor
1362 sum1 = 0.
1363 sum2 = 0.
1365 ! loop over exponential size distribution:
1366 jj_loop_6: do jj = 1,num_int_bins
1368 d1 = real(jj)*dd - 0.5*dd ! particle size
1370 ! get capacitance for different ice regimes:
1371 if (d1.le.dcrit) then
1372 cap = 1. ! for small spherical crystal use sphere
1373 elseif (d1.gt.dcrit.and.d1.le.dcrits) then
1374 cap = 0.48 ! field et al. 2006
1375 elseif (d1.gt.dcrits.and.d1.le.dcritr) then
1376 cap = 1. ! for graupel assume sphere
1377 elseif (d1.gt.dcritr) then
1378 cs1 = csr
1379 ds1 = dsr
1380 if (i_Fr.eq.1) then
1381 cap = 0.48
1382 else
1383 dum1 = 0.48
1384 dum2 = 1.
1385 m1 = cs1*d1**ds1
1386 m2 = cs*d1**ds
1387 m3 = cgp(i_rhor)*d1**dg
1388 ! linearly interpolate to get capacitance based on particle mass
1389 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
1390 cap = dum3
1391 endif
1392 endif
1394 ! for ventilation, only include fallspeed and size effects, the rest of
1395 ! term Sc^1/3 x Re^1/2 is multiplied in-line in the model code to allow
1396 ! effects of atmospheric conditions on ventilation
1397 !dum = (mu/dv)**0.333333*(fall1(jj)*d1/mu)**0.5
1398 dum = (fall1(jj)*d1)**0.5
1400 ! ventilation from Hall and Pruppacher (1976)
1401 ! only include ventilation for super-100 micron particles
1402 ! if (dum.lt.1.) then
1404 ! units are m^3 kg^-1 s^-1, thus multiplication by air density in P3_MAIN
1406 if (d1.lt.100.e-6) then
1407 sum1 = sum1+cap*n0*d1**(mu_i+1.)*exp(-lam*d1)*dd
1408 else
1409 !sum1 = sum1+cap*n0*(0.65+0.44*dum)*d1**(mu_i+1.)*exp(-lam*d1)*dd
1410 sum1 = sum1+cap*n0*0.65*d1**(mu_i+1.)*exp(-lam*d1)*dd
1411 sum2 = sum2+cap*n0*0.44*dum*d1**(mu_i+1.)*exp(-lam*d1)*dd
1412 endif
1414 enddo jj_loop_6
1416 vdep(i_Qnorm,i_Fr) = sum1
1417 vdep1(i_Qnorm,i_Fr) = sum2
1419 !.....................................................................................
1420 ! ice effective radius
1421 ! use definition of Francis et al. (1994), e.g., Eq. 3.11 in Fu (1996) J. Climate
1422 !.....................................................................................
1423 sum1 = 0.
1424 sum2 = 0.
1426 ! loop over exponential size distribution:
1427 jj_loop_7: do jj = 1,num_int_bins
1429 d1 = real(jj)*dd - 0.5*dd ! particle size
1431 if (d1.le.dcrit) then
1432 cs1 = pi*sxth*900.
1433 ds1 = 3.
1434 bas1 = 2.
1435 aas1 = pi*0.25
1436 elseif (d1.gt.dcrit.and.d1.le.dcrits) then
1437 cs1 = cs
1438 ds1 = ds
1439 bas1 = bas
1440 aas1 = aas
1441 elseif (d1.gt.dcrits.and.d1.le.dcritr) then
1442 cs1 = cgp(i_rhor)
1443 ds1 = dg
1444 bas1 = bag
1445 aas1 = aag
1446 elseif (d1.gt.dcritr) then
1447 cs1 = csr
1448 ds1 = dsr
1449 if (i_Fr.eq.1) then
1450 bas1 = bas
1451 aas1 = aas
1452 else
1453 ! for area, keep bas1 constant, but modify aas1 according to rime fraction
1454 bas1 = bas
1455 dum1 = aas*d1**bas
1456 dum2 = aag*d1**bag
1457 m1 = cs1*d1**ds1
1458 m2 = cs*d1**ds
1459 m3 = cgp(i_rhor)*d1**dg
1460 ! linearly interpolate based on particle mass:
1461 dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2)
1462 aas1 = dum3/(d1**bas)
1463 endif
1464 endif
1466 ! n0 not included below becuase it is in both numerator and denominator
1467 !cs1 = pi*sxth*917.
1468 !ds1 = 3.
1469 !aas1 = pi/4.*2.
1470 !bas1 = 2.
1472 sum1 = sum1 + cs1*d1**ds1*d1**mu_i*exp(-lam*d1)*dd
1473 sum2 = sum2 + aas1*d1**bas1*d1**mu_i*exp(-lam*d1)*dd
1475 ! if (d1.ge.100.e-6) then
1476 ! sum3 = sum3+n0*d1**mu_i*exp(-lam*d1)*dd
1477 ! sum4 = sum4+n0*aas1*d1**bas1*d1**mu_i*exp(-lam*d1)*dd
1478 ! endif
1480 enddo jj_loop_7
1482 ! calculate eff radius:
1484 !eff(i_Qnorm,i_Fr) = sum1/(1.7321*916.7*sum2)
1485 !calculate effective size following Fu (1996)
1486 !eff(i_Qnorm,i_Fr) = sum1/(1.1547*916.7*sum2)
1487 ! calculate for eff rad for twp ice:
1488 eff(i_Qnorm,i_Fr) = 3.*sum1/(4.*sum2*916.7)
1490 !.....................................................................................
1492 522 continue
1494 enddo i_Qnorm_loop
1497 !-- ice table
1498 i_Qnorm_loop_2: do i_Qnorm = 1,n_Qnorm
1501 ! Set values less than cutoff (1.e-99) to 0.
1502 ! note: dim(x,cutoff) actually returns x-cutoff (if x>cutoff; else 0.), but this difference will
1503 ! have no effect since the values will be read in single precision in P3_INIT. The purppse
1504 ! here is to avoid problems trying to writevalues with 3-digit exponents (e.g. 0.123456E-100)
1505 uns(i_Qnorm,i_Fr) = dim( uns(i_Qnorm,i_Fr), cutoff)
1506 ums(i_Qnorm,i_Fr) = dim( ums(i_Qnorm,i_Fr), cutoff)
1507 nagg(i_Qnorm,i_Fr) = dim( nagg(i_Qnorm,i_Fr), cutoff)
1508 nrwat(i_Qnorm,i_Fr) = dim( nrwat(i_Qnorm,i_Fr), cutoff)
1509 vdep(i_Qnorm,i_Fr) = dim( vdep(i_Qnorm,i_Fr), cutoff)
1510 eff(i_Qnorm,i_Fr) = dim( eff(i_Qnorm,i_Fr), cutoff)
1511 i_qsmall(i_Qnorm,i_Fr) = dim( i_qsmall(i_Qnorm,i_Fr), cutoff)
1512 i_qlarge(i_Qnorm,i_Fr) = dim( i_qlarge(i_Qnorm,i_Fr), cutoff)
1513 refl(i_Qnorm,i_Fr) = dim( refl(i_Qnorm,i_Fr), cutoff)
1514 vdep1(i_Qnorm,i_Fr) = dim( vdep1(i_Qnorm,i_Fr), cutoff)
1515 dmm(i_Qnorm,i_Fr) = dim( dmm(i_Qnorm,i_Fr), cutoff)
1516 rhomm(i_Qnorm,i_Fr) = dim( rhomm(i_Qnorm,i_Fr), cutoff)
1517 uzs(i_Qnorm,i_Fr) = dim( uzs(i_Qnorm,i_Fr), cutoff)
1518 ! HM no longer needed
1519 ! zlarge(i_Qnorm,i_Fr) = dim( zlarge(i_Qnorm,i_Fr), cutoff)
1520 ! zsmall(i_Qnorm,i_Fr) = dim( zsmall(i_Qnorm,i_Fr), cutoff)
1521 lambda_i(i_Qnorm,i_Fr) = dim( lambda_i(i_Qnorm,i_Fr), cutoff)
1522 mu_i_save(i_Qnorm,i_Fr) = dim( mu_i_save(i_Qnorm,i_Fr), cutoff)
1524 if (log_3momI) then
1525 write(1,'(4i5,20e15.5)') &
1526 i_Znorm,i_rhor,i_Fr,i_Qnorm, &
1527 uns(i_Qnorm,i_Fr), &
1528 ums(i_Qnorm,i_Fr), &
1529 nagg(i_Qnorm,i_Fr), &
1530 nrwat(i_Qnorm,i_Fr), &
1531 vdep(i_Qnorm,i_Fr), &
1532 eff(i_Qnorm,i_Fr), &
1533 i_qsmall(i_Qnorm,i_Fr), &
1534 i_qlarge(i_Qnorm,i_Fr), &
1535 refl(i_Qnorm,i_Fr), &
1536 vdep1(i_Qnorm,i_Fr), &
1537 dmm(i_Qnorm,i_Fr), &
1538 rhomm(i_Qnorm,i_Fr), &
1539 uzs(i_Qnorm,i_Fr), &
1540 ! HM no longer needed
1541 ! zlarge(i_Qnorm,i_Fr), &
1542 ! zsmall(i_Qnorm,i_Fr), &
1543 lambda_i(i_Qnorm,i_Fr), &
1544 mu_i_save(i_Qnorm,i_Fr)
1545 else
1546 write(1,'(3i5,20e15.5)') &
1547 i_rhor,i_Fr,i_Qnorm, &
1548 uns(i_Qnorm,i_Fr), &
1549 ums(i_Qnorm,i_Fr), &
1550 nagg(i_Qnorm,i_Fr), &
1551 nrwat(i_Qnorm,i_Fr), &
1552 vdep(i_Qnorm,i_Fr), &
1553 eff(i_Qnorm,i_Fr), &
1554 i_qsmall(i_Qnorm,i_Fr), &
1555 i_qlarge(i_Qnorm,i_Fr), &
1556 refl(i_Qnorm,i_Fr), &
1557 vdep1(i_Qnorm,i_Fr), &
1558 dmm(i_Qnorm,i_Fr), &
1559 rhomm(i_Qnorm,i_Fr), &
1560 lambda_i(i_Qnorm,i_Fr), &
1561 mu_i_save(i_Qnorm,i_Fr)
1562 endif
1564 enddo i_Qnorm_loop_2
1566 !-- ice-rain collection table:
1567 do i_Qnorm = 1,n_Qnorm
1568 do i_Drscale = 1,n_Drscale
1570 ! ! !Set values less than cutoff (1.e-99) to 0.
1571 ! ! nrrain(i_Qnorm,i_Drscale,i_Fr) = dim(nrrain(i_Qnorm,i_Drscale,i_Fr), cutoff)
1572 ! ! qrrain(i_Qnorm,i_Drscale,i_Fr) = dim(qrrain(i_Qnorm,i_Drscale,i_Fr), cutoff)
1573 nrrain(i_Qnorm,i_Drscale,i_Fr) = log10(max(nrrain(i_Qnorm,i_Drscale,i_Fr), 1.e-99))
1574 qrrain(i_Qnorm,i_Drscale,i_Fr) = log10(max(qrrain(i_Qnorm,i_Drscale,i_Fr), 1.e-99))
1576 write(1,'(3i5,2e15.5)') &
1577 i_Qnorm,i_Drscale,i_Fr, &
1578 nrrain(i_Qnorm,i_Drscale,i_Fr), &
1579 qrrain(i_Qnorm,i_Drscale,i_Fr)
1581 enddo !i_Drscale-loop
1582 enddo !i_Qnorm-loop
1584 !--
1585 ! The values of i_Znorm (3-momI) and i_rhor/i_Fr (2-momI) are "passed in" for parallelized
1586 ! version of code, thus the loops are commented out.
1587 enddo i_Fr_loop_1
1588 ! enddo i_rhor_loop
1589 ! enddo i_Znorm_loop
1590 !==
1592 close(1)
1594 END PROGRAM create_p3_lookuptable_1
1595 !______________________________________________________________________________________
1597 ! Incomplete gamma function
1598 ! from Numerical Recipes in Fortran 77: The Art of
1599 ! Scientific Computing
1601 function gammq(a,x)
1603 real a,gammq,x
1605 ! USES gcf,gser
1606 ! Returns the incomplete gamma function Q(a,x) = 1-P(a,x)
1608 real gammcf,gammser,gln
1609 ! if (x.lt.0..or.a.le.0) pause 'bad argument in gammq'
1610 if (x.lt.0..or.a.le.0) print*, 'bad argument in gammq'
1611 if (x.lt.a+1.) then
1612 call gser(gamser,a,x,gln)
1613 gammq=1.-gamser
1614 else
1615 call gcf(gammcf,a,x,gln)
1616 gammq=gammcf
1617 end if
1618 return
1619 end
1621 !______________________________________________________________________________________
1623 subroutine gser(gamser,a,x,gln)
1624 integer itmax
1625 real a,gamser,gln,x,eps
1626 parameter(itmax=100,eps=3.e-7)
1627 integer n
1628 real ap,del,sum,gamma
1629 gln = log(gamma(a))
1630 if (x.le.0.) then
1631 ! if (x.lt.0.) pause 'x < 0 in gser'
1632 if (x.lt.0.) print*, 'x < 0 in gser'
1633 gamser = 0.
1634 return
1635 end if
1636 ap=a
1637 sum=1./a
1638 del=sum
1639 do n=1,itmax
1640 ap=ap+1.
1641 del=del*x/ap
1642 sum=sum+del
1643 if (abs(del).lt.abs(sum)*eps) goto 1
1644 end do
1645 ! pause 'a too large, itmax too small in gser'
1646 print*, 'a too large, itmax too small in gser'
1647 1 gamser=sum*exp(-x+a*log(x)-gln)
1648 return
1649 end
1651 !______________________________________________________________________________________
1653 subroutine gcf(gammcf,a,x,gln)
1654 integer itmax
1655 real a,gammcf,gln,x,eps,fpmin
1656 parameter(itmax=100,eps=3.e-7,fpmin=1.e-30)
1657 integer i
1658 real an,b,c,d,del,h,gamma
1659 gln=log(gamma(a))
1660 b=x+1.-a
1661 c=1./fpmin
1662 d=1./b
1663 h=d
1664 do i=1,itmax
1665 an=-i*(i-a)
1666 b=b+2.
1667 d=an*d+b
1668 if(abs(d).lt.fpmin) d=fpmin
1669 c=b+an/c
1670 if(abs(c).lt.fpmin) c=fpmin
1671 d=1./d
1672 del=d*c
1673 h = h*del
1674 if(abs(del-1.).lt.eps)goto 1
1675 end do
1676 ! pause 'a too large, itmax too small in gcf'
1677 print*, 'a too large, itmax too small in gcf'
1678 1 gammcf=exp(-x+a*log(x)-gln)*h
1679 return
1680 end
1682 !______________________________________________________________________________________
1684 real function compute_mu_3moment(mom3,mom6,mu_max)
1686 !--------------------------------------------------------------------------
1687 ! Computes mu as a function of G(mu), where
1688 !
1689 ! G(mu)= N*Z/Q^2 = [(6+mu)(5+mu)(4+mu)]/[(3+mu)(2+mu)(1+mu)]
1690 !
1691 ! 2018-08-08
1692 !--------------------------------------------------------------------------
1694 implicit none
1696 ! Arguments passed:
1697 real, intent(in) :: mom3,mom6 !normalized moments
1698 real, intent(in) :: mu_max !maximum allowable value of mu
1700 ! Local variables:
1701 real :: mu ! shape parameter in gamma distribution
1702 real :: a1,g1,g2,G
1704 ! calculate G from normalized moments
1705 G = mom6/mom3
1707 !----------------------------------------------------------!
1708 ! !Solve alpha numerically: (brute-force)
1709 ! mu= 0.
1710 ! g2= 999.
1711 ! do i=0,4000
1712 ! a1= i*0.01
1713 ! g1= (6.+a1)*(5.+a1)*(4.+a1)/((3.+a1)*(2.+a1)*(1.+a1))
1714 ! if(abs(g-g1)<abs(g-g2)) then
1715 ! mu = a1
1716 ! g2= g1
1717 ! endif
1718 ! enddo
1719 !----------------------------------------------------------!
1721 !Piecewise-polynomial approximation of G(mu) to solve for mu:
1722 if (G >= 20.) then
1723 mu = 0.
1724 else
1725 g2 = G**2
1726 if (G<20. .and.G>=13.31) mu = 3.3638e-3*g2 - 1.7152e-1*G + 2.0857e+0
1727 if (G<13.31.and.G>=7.123) mu = 1.5900e-2*g2 - 4.8202e-1*G + 4.0108e+0
1728 if (G<7.123.and.G>=4.200) mu = 1.0730e-1*g2 - 1.7481e+0*G + 8.4246e+0
1729 if (G<4.200.and.G>=2.946) mu = 5.9070e-1*g2 - 5.7918e+0*G + 1.6919e+1
1730 if (G<2.946.and.G>=1.793) mu = 4.3966e+0*g2 - 2.6659e+1*G + 4.5477e+1
1731 if (G<1.793.and.G>=1.405) mu = 4.7552e+1*g2 - 1.7958e+2*G + 1.8126e+2
1732 if (G<1.405.and.G>=1.230) mu = 3.0889e+2*g2 - 9.0854e+2*G + 6.8995e+2
1733 if (G<1.230) mu = mu_max
1734 endif
1736 compute_mu_3moment = mu
1738 end function compute_mu_3moment
1740 !______________________________________________________________________________________
1742 real function diagnostic_mui(mu_i_min,mu_i_max,lam,q,cgp,Fr,pi)
1744 !----------------------------------------------------------!
1745 ! Compute mu_i diagnostically.
1746 !----------------------------------------------------------!
1748 implicit none
1750 !Arguments:
1751 real :: mu_i_min,mu_i_max,lam,q,cgp,Fr,pi
1753 ! Local variables:
1754 real, parameter :: Di_thres = 0.2 !diameter threshold [mm]
1755 !real, parameter :: Di_thres = 0.6 !diameter threshold [mm]
1756 real :: mu_i,dum1,dum2,dum3
1759 !-- diagnostic mu_i, original formulation: (from Heymsfield, 2003)
1760 ! mu_i = 0.076*(lam/100.)**0.8-2. ! /100 is to convert m-1 to cm-1
1761 ! mu_i = max(mu_i,0.)
1762 ! mu_i = min(mu_i,6.)
1764 !-- diagnostic mu_i, 3-moment-based formulation:
1765 ! dum1 = (q/cgp)**(1./3)*1000. ! estimated Dmvd [mm], assuming spherical
1766 ! if (dum1<=Di_thres) then
1767 ! !diagnostic mu_i, original formulation: (from Heymsfield, 2003)
1768 ! mu_i = 0.076*(lam*0.01)**0.8-2. ! /100 is to convert m-1 to cm-1
1769 ! mu_i = min(mu_i,6.)
1770 ! else
1771 ! dum2 = (6./pi)*cgp ! mean density (total)
1772 ! dum3 = max(1., 1.+0.00842*(dum2-400.)) ! adjustment factor for density
1773 ! mu_i = 4.*(dum1-Di_thres)*dum3*Fr
1774 ! endif
1775 ! mu_i = max(mu_i,mu_i_min) ! make sure mu_i >= 0, otherwise size dist is infinity at D = 0
1776 ! mu_i = min(mu_i,mu_i_max)
1778 dum1 = (q/cgp)**(1./3)*1000. ! estimated Dmvd [mm], assuming spherical
1779 if (dum1<=Di_thres) then
1780 !diagnostic mu_i, original formulation: (from Heymsfield, 2003)
1781 mu_i = 0.076*(lam*0.01)**0.8-2. ! /100 is to convert m-1 to cm-1
1782 mu_i = min(mu_i,6.)
1783 else
1784 dum2 = (6./pi)*cgp ! mean density (total)
1785 ! dum3 = max(1., 1.+0.00842*(dum2-400.)) ! adjustment factor for density
1786 dum3 = max(1., 1.+0.00842*(dum2-400.)) ! adjustment factor for density
1787 mu_i = 0.25*(dum1-Di_thres)*dum3*Fr
1788 endif
1789 mu_i = max(mu_i,mu_i_min) ! make sure mu_i >= 0, otherwise size dist is infinity at D = 0
1790 mu_i = min(mu_i,mu_i_max)
1792 diagnostic_mui = mu_i
1794 end function diagnostic_mui
1796 !______________________________________________________________________________________
1798 subroutine intgrl_section(lam,mu, d1,d2,d3,d4, Dcrit1,Dcrit2,Dcrit3, &
1799 intsec_1,intsec_2,intsec_3,intsec_4)
1800 !-----------------
1801 ! Computes and returns partial integrals (partial moments) of ice PSD.
1802 !-----------------
1804 implicit none
1806 !Arguments:
1807 real, intent(in) :: lam,mu, d1,d2,d3,d4, Dcrit1,Dcrit2,Dcrit3
1808 real, intent(out) :: intsec_1,intsec_2,intsec_3,intsec_4
1810 !Local:
1811 real :: dum,gammq
1812 !-----------------
1814 !Region I -- integral from 0 to Dcrit1 (small spherical ice)
1815 intsec_1 = lam**(-d1-mu-1.)*gamma(mu+d1+1.)*(1.-gammq(mu+d1+1.,Dcrit1*lam))
1817 !Region II -- integral from Dcrit1 to Dcrit2 (non-spherical unrimed ice)
1818 intsec_2 = lam**(-d2-mu-1.)*gamma(mu+d2+1.)*(gammq(mu+d2+1.,Dcrit1*lam))
1819 dum = lam**(-d2-mu-1.)*gamma(mu+d2+1.)*(gammq(mu+d2+1.,Dcrit2*lam))
1820 intsec_2 = intsec_2-dum
1822 !Region III -- integral from Dcrit2 to Dcrit3 (fully rimed spherical ice)
1823 intsec_3 = lam**(-d3-mu-1.)*gamma(mu+d3+1.)*(gammq(mu+d3+1.,Dcrit2*lam))
1824 dum = lam**(-d3-mu-1.)*gamma(mu+d3+1.)*(gammq(mu+d3+1.,Dcrit3*lam))
1825 intsec_3 = intsec_3-dum
1827 !Region IV -- integral from Dcrit3 to infinity (partially rimed ice)
1828 intsec_4 = lam**(-d4-mu-1.)*gamma(mu+d4+1.)*(gammq(mu+d4+1.,Dcrit3*lam))
1830 return
1832 end subroutine intgrl_section