Radiation Scheme in CCPP
Functions/Subroutines | Variables
module_microphysics Module Reference

Functions/Subroutines

subroutine gsmconst (DTPG, mype, first)
 
subroutine my_growth_rates (DTPH)
 
subroutine ice_lookup
 
subroutine rain_lookup
 
subroutine gsmcolumn (ARAING, ASNOWG, DTPG, I_index, J_index,
 
subroutine rsipath (im, ix, ix2, levs, prsl, prsi, t, q, clw
 
subroutine rsipath2
 This subroutine is a modified version of ferrier's original "rsipath" subprogram. It computes layer's cloud liquid, ice, rain, and snow water condensate path and the partical effective radius for liquid droplet, rain drop, and snow flake. More...
 

Variables

real, private abfr
 
real, private cbfr
 
real, private ciacw
 
real, private ciacr
 
real, private c_n0r0
 
integer, private mic_step
 
integer, parameter, private my_t1 =1
 
integer, parameter, private my_t2 =35
 
real, dimension(my_t1:my_t2), private my_growth
 
real, parameter, private dmimin =.05e-3
 
real, parameter, private dmimax =1.e-3
 
integer, parameter, private mdimin =XMImin
 
integer, parameter, private mdimax =XMImax
 
real, parameter, private dmrmin =.05e-3
 
real, parameter, private dmrmax =.45e-3
 
integer, parameter, private mdrmin =XMRmin
 
integer, parameter, private mdrmax =XMRmax
 
integer, parameter, private indexsmin =100
 
real, parameter, private rerainmin =1.5*XMRmin
 
integer, parameter, private nrime =40
 
real, dimension(2:9, 0:nrime), private vel_rf
 
integer, parameter itlo =-60
 
integer, parameter ithi =40
 
integer, dimension(itlo:ithi, 4) nstats
 
real, dimension(itlo:ithi, 5) qmax
 
real, dimension(itlo:ithi, 22) qtot
 
real, parameter, private thom =T_ICE
 
real, parameter, private tnw =50.
 
real, parameter, private toler =1.0E-20
 

Function/Subroutine Documentation

subroutine module_microphysics::gsmcolumn ( real  ARAING,
real  ASNOWG,
  DTPG,
integer  I_index,
integer  J_index 
)

Definition at line 1132 of file module_bfmicrophysics.f.

References abfr, cbfr, ciacr, ciacw, condense(), deposit(), dmimax, ithi, itlo, mdimax, mdrmax, mdrmin, mic_step, my_growth, my_t1, my_t2, nrime, nstats, qmax, qtot, toler, and vel_rf.

Here is the call graph for this function:

subroutine module_microphysics::gsmconst ( real  DTPG,
integer  mype,
logical  first 
)

Definition at line 119 of file module_bfmicrophysics.f.

References abfr, c_n0r0, cbfr, ciacr, ciacw, dmrmax, dmrmin, ice_lookup(), mdimax, mdimin, mdrmax, mdrmin, mic_step, my_growth_rates(), nrime, rain_lookup(), and vel_rf.

119 !
120  implicit none
121 !-------------------------------------------------------------------------------
122 !--- SUBPROGRAM DOCUMENTATION BLOCK
123 ! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001
124 !-------------------------------------------------------------------------------
125 ! ABSTRACT:
126 ! * Reads various microphysical lookup tables used in COLUMN_MICRO
127 ! * Lookup tables were created "offline" and are read in during execution
128 ! * Creates lookup tables for saturation vapor pressure w/r/t water & ice
129 !-------------------------------------------------------------------------------
130 !
131 ! USAGE: CALL GSMCONST FROM SUBROUTINE GSMDRIVE AT MODEL START TIME
132 !
133 ! INPUT ARGUMENT LIST:
134 ! DTPH - physics time step (s)
135 !
136 ! OUTPUT ARGUMENT LIST:
137 ! NONE
138 !
139 ! OUTPUT FILES:
140 ! NONE
141 !
142 !
143 ! SUBROUTINES:
144 ! MY_GROWTH_RATES - lookup table for growth of nucleated ice
145 !
146 ! UNIQUE: NONE
147 !
148 ! LIBRARY: NONE
149 !
150 !
151 ! ATTRIBUTES:
152 ! LANGUAGE: FORTRAN 90
153 ! MACHINE : IBM SP
154 !
155  integer mype
156  real dtpg
157  logical first
158 !
159 !--- Parameters & data statement for local calculations
160 !
161  REAL, PARAMETER :: c1=1./3., dmr1=.1e-3, dmr2=.2e-3, dmr3=.32e-3, &
162  & n0r0=8.e6, n0s0=4.e6, rhol=1000., rhos=100., &
163  & xmr1=1.e6*dmr1, xmr2=1.e6*dmr2, xmr3=1.e6*dmr3
164  INTEGER, PARAMETER :: mdr1=xmr1, mdr2=xmr2, mdr3=xmr3
165 !
166  real dtph, bbfr
167  integer i
168 !
169 !--- Added on 5/16/01 for Moorthi
170 !
171  logical, parameter :: read_lookup=.false., write_lookup=.false.
172 !
173 !------------------------------------------------------------------------
174 ! ************* Parameters used in ETA model -- Not used in Global Model *****
175 !
176 !--- DPHD, DLMD are delta latitude and longitude at the model (NOT geodetic) equator
177 ! => "DX" is the hypotenuse of the model zonal & meridional grid increments.
178 !
179 ! DX=111.*(DPHD**2+DLMD**2)**.5 ! Resolution at MODEL equator (km)
180 ! DX=MIN(100., MAX(5., DX) )
181 !
182 !--- Assume the following functional relationship for key constants that
183 ! depend on grid resolution from DXmin (5 km) to DXmax (100 km) resolution:
184 !
185 ! DXmin=5.
186 ! DXmax=100.
187 ! DX=MIN(DXmax, MAX(DXmin, DX) )
188 !
189 !--- EXtune determines the degree to which the coefficients change with resolution.
190 ! The larger EXtune is, the more sensitive the parameter.
191 !
192 ! EXtune=1.
193 
194 !
195 !--- FXtune ==> F(DX) is the grid-resolution tuning parameter (from 0 to 1)
196 !
197 ! FXtune=((DXmax-DX)/(DXmax-DXmin))**EXtune
198 ! FXtune=MAX(0., MIN(1., FXtune))
199 !
200 !--- Calculate grid-averaged RH for the onset of condensation (RHgrd) based on
201 ! simple ***ASSUMED*** (user-specified) values at DXmax and at DXmin.
202 !
203 ! RH_DXmax=.90 !-- 90% RH at DXmax=100 km
204 ! RH_DXmin=.98 !-- 98% RH at DXmin=5 km
205 !
206 !--- Note that RHgrd is right now fixed throughout the domain!!
207 !
208 ! RHgrd=RH_DXmax+(RH_DXmin-RH_DXmax)*FXtune
209 ! ********************************************************************************
210 !
211 !
212  if (first) then
213 !
214 !--- Read in various lookup tables
215 !
216  if ( read_lookup ) then
217  OPEN (unit=1,file="eta_micro_lookup.dat",form="UNFORMATTED")
218  READ(1) ventr1
219  READ(1) ventr2
220  READ(1) accrr
221  READ(1) massr
222  READ(1) vrain
223  READ(1) rrate
224  READ(1) venti1
225  READ(1) venti2
226  READ(1) accri
227  READ(1) massi
228  READ(1) vsnowi
229  READ(1) vel_rf
230 ! read(1) my_growth ! Applicable only for DTPH=180 s for offline testing
231  CLOSE (1)
232  else
233  CALL ice_lookup ! Lookup tables for ice
234  CALL rain_lookup ! Lookup tables for rain
235  if (write_lookup) then
236  open(unit=1,file='micro_lookup.dat',form='unformatted')
237  write(1) ventr1
238  write(1) ventr2
239  write(1) accrr
240  write(1) massr
241  write(1) vrain
242  write(1) rrate
243  write(1) venti1
244  write(1) venti2
245  write(1) accri
246  write(1) massi
247  write(1) vsnowi
248  write(1) vel_rf
249 ! write(1) my_growth ! Applicable only for DTPH=180 s ????
250  CLOSE (1)
251  endif
252  endif
253 !!
254 !--- Constants associated with Biggs (1953) freezing of rain, as parameterized
255 ! following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
256 !
257  abfr=-0.66
258  bbfr=100.
259  cbfr=20.*pi*pi*bbfr*rhol*1.e-21
260 !
261 !--- QAUT0 is the threshold cloud content for autoconversion to rain
262 ! needed for droplets to reach a diameter of 20 microns (following
263 ! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM). It is
264 ! **STRONGLY** affected by the assumed droplet number concentrations
265 ! XNCW! For example, QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for
266 ! droplet number concentrations of 300, 200, and 100 cm**-3, respectively.
267 !
268 !--- Calculate grid-averaged XNCW based on simple ***ASSUMED*** (user-specified)
269 ! values at DXmax and at DXmin.
270 !
271 ! XNCW_DXmax=50.E6 !-- 50 /cm**3 at DXmax=100 km
272 ! XNCW_DXmin=200.E6 !-- 200 /cm**3 at DXmin=5 km
273 !
274 !--- Note that XNCW is right now fixed throughout the domain!!
275 !
276 ! XNCW=XNCW_DXmax+(XNCW_DXmin-XNCW_DXmax)*FXtune
277 !
278 ! QAUT0=PI*RHOL*XNCW*(20.E-6)**3/6.
279  qautx=pi*rhol*1.0e6*(20.e-6)**3/6.
280 !
281 !--- Based on rain lookup tables for mean diameters from 0.05 to 0.45 mm
282 ! * Four different functional relationships of mean drop diameter as
283 ! a function of rain rate (RR), derived based on simple fits to
284 ! mass-weighted fall speeds of rain as functions of mean diameter
285 ! from the lookup tables.
286 !
287  rr_drmin=n0r0*rrate(mdrmin) ! RR for mean drop diameter of .05 mm
288  rr_dr1=n0r0*rrate(mdr1) ! RR for mean drop diameter of .10 mm
289  rr_dr2=n0r0*rrate(mdr2) ! RR for mean drop diameter of .20 mm
290  rr_dr3=n0r0*rrate(mdr3) ! RR for mean drop diameter of .32 mm
291  rr_drmax=n0r0*rrate(mdrmax) ! RR for mean drop diameter of .45 mm
292 !
293  rqr_drmin=n0r0*massr(mdrmin) ! Rain content for mean drop diameter of .05 mm
294  rqr_dr1=n0r0*massr(mdr1) ! Rain content for mean drop diameter of .10 mm
295  rqr_dr2=n0r0*massr(mdr2) ! Rain content for mean drop diameter of .20 mm
296  rqr_dr3=n0r0*massr(mdr3) ! Rain content for mean drop diameter of .32 mm
297  rqr_drmax=n0r0*massr(mdrmax) ! Rain content for mean drop diameter of .45 mm
298  c_n0r0=pi*rhol*n0r0
299  cn0r0=1.e6/c_n0r0**.25
300  cn0r_dmrmin=1./(pi*rhol*dmrmin**4)
301  cn0r_dmrmax=1./(pi*rhol*dmrmax**4)
302 !
303  endif ! If (first) then loop ends here
304 !
305 ! Find out what microphysics time step should be
306 !
307  mic_step = max(1, int(dtpg/600.0+0.5))
308 ! mic_step = max(1, int(dtpg/300.0+0.5))
309  dtph = dtpg / mic_step
310  if (mype == 0) print *,' DTPG=',dtpg,' mic_step=',mic_step &
311  &, ' dtph=',dtph
312 !
313 !--- Calculates coefficients for growth rates of ice nucleated in water
314 ! saturated conditions, scaled by physics time step (lookup table)
315 !
316  CALL my_growth_rates (dtph)
317 !
318 !--- CIACW is used in calculating riming rates
319 ! The assumed effective collection efficiency of cloud water rimed onto
320 ! ice is =0.5 below:
321 !
322 !Moor CIACW=DTPH*0.25*PI*0.5*(1.E5)**C1 ! commented on 20050422
323 ! ice is =0.1 below:
324  ciacw=dtph*0.25*pi*0.1*(1.e5)**c1
325 ! CIACW = 0.0 ! Brad's suggestion 20040614
326 !
327 !--- CIACR is used in calculating freezing of rain colliding with large ice
328 ! The assumed collection efficiency is 1.0
329 !
330  ciacr=pi*dtph
331 !
332 !--- CRACW is used in calculating collection of cloud water by rain (an
333 ! assumed collection efficiency of 1.0)
334 !
335 !Moor CRACW=DTPH*0.25*PI*1.0 ! commented on 20050422
336 !
337 ! assumed collection efficiency of 0.1)
338  cracw=dtph*0.25*pi*0.1
339 ! CRACW = 0.0 ! Brad's suggestion 20040614
340 !
341  esw0=fpvsl(t0c) ! Saturation vapor pressure at 0C
342  rfmax=1.1**nrime ! Maximum rime factor allowed
343 !
344 !------------------------------------------------------------------------
345 !--------------- Constants passed through argument list -----------------
346 !------------------------------------------------------------------------
347 !
348 !--- Important parameters for self collection (autoconversion) of
349 ! cloud water to rain.
350 !
351 !--- CRAUT is proportional to the rate that cloud water is converted by
352 ! self collection to rain (autoconversion rate)
353 !
354  craut=1.-exp(-1.e-3*dtph)
355 !
356 ! IF (MYPE == 0)
357 ! & WRITE(6,"(A, A,F6.2,A, A,F5.4, A,F7.3,A, A,F6.2,A, A,F6.3,A)")
358 ! & 'KEY MICROPHYSICAL PARAMETERS FOR '
359 ! & ,'DX=',DX,' KM:'
360 ! & ,' FXtune=',FXtune
361 ! & ,' RHgrd=',100.*RHgrd,' %'
362 ! & ,' NCW=',1.E-6*XNCW,' /cm**3'
363 ! & ,' QAUT0=',1.E3*QAUT0,' g/kg'
364 !
365 !--- For calculating snow optical depths by considering bulk density of
366 ! snow based on emails from Q. Fu (6/27-28/01), where optical
367 ! depth (T) = 1.5*SWP/(Reff*DENS), SWP is snow water path, Reff
368 ! is effective radius, and DENS is the bulk density of snow.
369 !
370 ! SWP (kg/m**2)=(1.E-3 kg/g)*SWPrad, SWPrad in g/m**2 used in radiation
371 ! T = 1.5*1.E3*SWPrad/(Reff*DENS)
372 !
373 ! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
374 !
375 ! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI*(1.E-6*INDEXS)**3]
376 !
377  DO i=mdimin,mdimax
378 !MoorthiSDENS(I)=PI*1.5E-15*FLOAT(I*I*I)/MASSI(I)
379  sdens(i)=pi*1.0e-15*float(i*i*i)/massi(i)
380  ENDDO
381 !
382 !-----------------------------------------------------------------------
383 !

Here is the call graph for this function:

subroutine module_microphysics::ice_lookup ( )

Definition at line 434 of file module_bfmicrophysics.f.

References mdimax, mdimin, nrime, and vel_rf.

Referenced by gsmconst().

434 !
435  implicit none
436 !-----------------------------------------------------------------------------------
437 !
438 !---- Key diameter values in mm
439 !
440 !-----------------------------------------------------------------------------------
441 !
442 !---- Key concepts:
443 ! - Actual physical diameter of particles (D)
444 ! - Ratio of actual particle diameters to mean diameter (x=D/MD)
445 ! - Mean diameter of exponentially distributed particles, which is the
446 ! same as 1./LAMDA of the distribution (MD)
447 ! - All quantitative relationships relating ice particle characteristics as
448 ! functions of their diameter (e.g., ventilation coefficients, normalized
449 ! accretion rates, ice content, and mass-weighted fall speeds) are a result
450 ! of using composite relationships for ice crystals smaller than 1.5 mm
451 ! diameter merged with analogous relationships for larger sized aggregates.
452 ! Relationships are derived as functions of mean ice particle sizes assuming
453 ! exponential size spectra and assuming the properties of ice crystals at
454 ! sizes smaller than 1.5 mm and aggregates at larger sizes.
455 !
456 !-----------------------------------------------------------------------------------
457 !
458 !---- Actual ice particle diameters for which integrated distributions are derived
459 ! - DminI - minimum diameter for integration (.02 mm, 20 microns)
460 ! - DmaxI - maximum diameter for integration (2 cm)
461 ! - DdelI - interval for integration (1 micron)
462 !
463  real, parameter :: dmini=.02e-3, dmaxi=20.e-3, ddeli=1.e-6, &
464  & ximin=1.e6*dmini, ximax=1.e6*dmaxi
465  integer, parameter :: idimin=ximin, idimax=ximax
466 !
467 !---- Meaning of the following arrays:
468 ! - diam - ice particle diameter (m)
469 ! - mass - ice particle mass (kg)
470 ! - vel - ice particle fall speeds (m/s)
471 ! - vent1 - 1st term in ice particle ventilation factor
472 ! - vent2 - 2nd term in ice particle ventilation factor
473 !
474  real diam(idimin:idimax),mass(idimin:idimax),vel(idimin:idimax), &
475  & vent1(IDImin:IDImax),vent2(IDImin:IDImax)
476 !
477 !-----------------------------------------------------------------------------------
478 !
479 !---- Found from trial & error that the m(D) & V(D) mass & velocity relationships
480 ! between the ice crystals and aggregates overlapped & merged near a particle
481 ! diameter sizes of 1.5 mm. Thus, ice crystal relationships are used for
482 ! sizes smaller than 1.5 mm and aggregate relationships for larger sizes.
483 !
484  real, parameter :: d_crystal_max=1.5
485 !
486 !---- The quantity xmax represents the maximum value of "x" in which the
487 ! integrated values are calculated. For xmax=20., this means that
488 ! integrated ventilation, accretion, mass, and precipitation rates are
489 ! calculated for ice particle sizes less than 20.*mdiam, the mean particle diameter.
490 !
491  real, parameter :: xmax=20.
492 !
493 !-----------------------------------------------------------------------------------
494 !
495 !---- Meaning of the following arrays:
496 ! - mdiam - mean diameter (m)
497 ! - VENTI1 - integrated quantity associated w/ ventilation effects
498 ! (capacitance only) for calculating vapor deposition onto ice
499 ! - VENTI2 - integrated quantity associated w/ ventilation effects
500 ! (with fall speed) for calculating vapor deposition onto ice
501 ! - ACCRI - integrated quantity associated w/ cloud water collection by ice
502 ! - MASSI - integrated quantity associated w/ ice mass
503 ! - VSNOWI - mass-weighted fall speed of snow, used to calculate precip rates
504 !
505 !--- Mean ice-particle diameters varying from 50 microns to 1000 microns (1 mm),
506 ! assuming an exponential size distribution.
507 !
508  real mdiam
509 !
510 !-----------------------------------------------------------------------------------
511 !------------- Constants & parameters for ventilation factors of ice ---------------
512 !-----------------------------------------------------------------------------------
513 !
514 !---- These parameters are used for calculating the ventilation factors for ice
515 ! crystals between 0.2 and 1.5 mm diameter (Hall and Pruppacher, JAS, 1976).
516 ! From trial & error calculations, it was determined that the ventilation
517 ! factors of smaller ice crystals could be approximated by a simple linear
518 ! increase in the ventilation coefficient from 1.0 at 50 microns (.05 mm) to
519 ! 1.1 at 200 microns (0.2 mm), rather than using the more complex function of
520 ! 1.0 + .14*(Sc**.33*Re**.5)**2 recommended by Hall & Pruppacher.
521 !
522  real, parameter :: cvent1i=.86, cvent2i=.28
523 !
524 !---- These parameters are used for calculating the ventilation factors for larger
525 ! aggregates, where D>=1.5 mm (see Rutledge and Hobbs, JAS, 1983;
526 ! Thorpe and Mason, 1966).
527 !
528  real, parameter :: cvent1a=.65, cvent2a=.44
529 !
530  real m_agg,m_bullet,m_column,m_ice,m_plate
531 !
532 !---- Various constants
533 !
534  real, parameter :: c1=2./3., cexp=1./3.
535 !
536  logical :: iprint
537  logical, parameter :: print_diag=.false.
538 !
539 !-----------------------------------------------------------------------------------
540 !- Constants & parameters for calculating the increase in fall speed of rimed ice --
541 !-----------------------------------------------------------------------------------
542 !
543 !---- Constants & arrays for estimating increasing fall speeds of rimed ice.
544 ! Based largely on theory and results from Bohm (JAS, 1989, 2419-2427).
545 !
546 !-------------------- Standard atmosphere conditions at 1000 mb --------------------
547 !
548  real, parameter :: t_std=288., dens_std=1000.e2/(287.04*288.)
549 !
550 !---- These "bulk densities" are the actual ice densities in the ice portion of the
551 ! lattice. They are based on text associated w/ (12) on p. 2425 of Bohm (JAS,
552 ! 1989). Columns, plates, & bullets are assumed to have an average bulk density
553 ! of 850 kg/m**3. Aggregates were assumed to have a slightly larger bulk density
554 ! of 600 kg/m**3 compared with dendrites (i.e., the least dense, most "lacy" &
555 ! tenous ice crystal, which was assumed to be ~500 kg/m**3 in Bohm).
556 !
557  real, parameter :: dens_crystal=850., dens_agg=600.
558 !
559 !--- A value of Nrime=40 for a logarithmic ratio of 1.1 yields a maximum rime factor
560 ! of 1.1**40 = 45.26 that is resolved in these tables. This allows the largest
561 ! ice particles with a mean diameter of MDImax=1000 microns to achieve bulk
562 ! densities of 900 kg/m**3 for rimed ice.
563 !
564 ! integer, parameter :: Nrime=40
565  real m_rime, &
566  & rime_factor(0:Nrime), rime_vel(0:Nrime), &
567  & vel_rime(IDImin:IDImax,Nrime), ivel_rime(MDImin:MDImax,Nrime)
568 !
569  integer i, j, jj, k, icount
570  real c2, cbulk, cbulk_ice, px, dynvis_std, crime1 &
571  &, crime2, crime3, crime4, crime5, d, c_avg, c_agg &
572  &, c_bullet, c_column, c_plate, cl_agg, cl_bullet &
573  &, cl_column, cl_plate, v_agg, v_bullet, v_column &
574  &, v_plate, wd, ecc_column &
575  &, cvent1, cvent2, crime_best, rime_m1, rime_m2 &
576  &, x_rime, re_rime, smom3, pratei, expf &
577  &, bulk_dens, xmass, xmdiam, ecc_plate, dx
578 !
579 !-----------------------------------------------------------------------------------
580 !----------------------------- BEGIN EXECUTION -------------------------------------
581 !-----------------------------------------------------------------------------------
582 !
583 !
584  c2=1./sqrt(3.)
585 ! pi=acos(-1.)
586  cbulk=6./pi
587  cbulk_ice=900.*pi/6. ! Maximum bulk ice density allowed of 900 kg/m**3
588  px=.4**cexp ! Convert fall speeds from 400 mb (Starr & Cox) to 1000 mb
589 !
590 !--------------------- Dynamic viscosity (1000 mb, 288 K) --------------------------
591 !
592  dynvis_std=1.496e-6*t_std**1.5/(t_std+120.)
593  crime1=pi/24.
594  crime2=8.*9.81*dens_std/(pi*dynvis_std**2)
595  crime3=crime1*dens_crystal
596  crime4=crime1*dens_agg
597  crime5=dynvis_std/dens_std
598  do i=0,nrime
599  rime_factor(i)=1.1**i
600  enddo
601 !
602 !#######################################################################
603 ! Characteristics as functions of actual ice particle diameter
604 !#######################################################################
605 !
606 !---- M(D) & V(D) for 3 categories of ice crystals described by Starr
607 !---- & Cox (1985).
608 !
609 !---- Capacitance & characteristic lengths for Reynolds Number calculations
610 !---- are based on Young (1993; p. 144 & p. 150). c-axis & a-axis
611 !---- relationships are from Heymsfield (JAS, 1972; Table 1, p. 1351).
612 !
613  icount=60
614 !
615  if (print_diag) &
616  & write(7,"(2a)") '---- Increase in fall speeds of rimed ice', &
617  & ' particles as function of ice particle diameter ----'
618  do i=idimin,idimax
619  if (icount == 60 .and. print_diag) then
620  write(6,"(/2a/3a)") 'Particle masses (mg), fall speeds ', &
621  & '(m/s), and ventilation factors', &
622  & ' D(mm) CR_mass Mass_bull Mass_col Mass_plat ', &
623  & ' Mass_agg CR_vel V_bul CR_col CR_pla Aggreg', &
624  & ' Vent1 Vent2 '
625  write(7,"(3a)") ' <----------------------------------',&
626  & '--------------- Rime Factor --------------------------', &
627  & '--------------------------->'
628  write(7,"(a,23f5.2)") ' D(mm)',(rime_factor(k), k=1,5), &
629  & (rime_factor(k), k=6,40,2)
630  icount=0
631  endif
632  d=(float(i)+.5)*1.e-3 ! in mm
633  c_avg=0.
634  c_agg=0.
635  c_bullet=0.
636  c_column=0.
637  c_plate=0.
638  cl_agg=0.
639  cl_bullet=0.
640  cl_column=0.
641  cl_plate=0.
642  m_agg=0.
643  m_bullet=0.
644  m_column=0.
645  m_plate=0.
646  v_agg=0.
647  v_bullet=0.
648  v_column=0.
649  v_plate=0.
650  if (d < d_crystal_max) then
651 !
652 !---- This block of code calculates bulk characteristics based on average
653 ! characteristics of bullets, plates, & column ice crystals <1.5 mm size
654 !
655 !---- Mass-diameter relationships from Heymsfield (1972) & used
656 ! in Starr & Cox (1985), units in mg
657 !---- "d" is maximum dimension size of crystal in mm,
658 !
659 ! Mass of pure ice for spherical particles, used as an upper limit for the
660 ! mass of small columns (<~ 80 microns) & plates (<~ 35 microns)
661 !
662  m_ice=.48*d**3 ! Mass of pure ice for spherical particle
663 !
664  m_bullet=min(.044*d**3, m_ice)
665  m_column=min(.017*d**1.7, m_ice)
666  m_plate=min(.026*d**2.5, m_ice)
667 !
668  mass(i)=m_bullet+m_column+m_plate
669 !
670 !---- These relationships are from Starr & Cox (1985), applicable at 400 mb
671 !---- "d" is maximum dimension size of crystal in mm, dx in microns
672 !
673  dx=1000.*d ! Convert from mm to microns
674  if (dx <= 200.) then
675  v_column=8.114e-5*dx**1.585
676  v_bullet=5.666e-5*dx**1.663
677  v_plate=1.e-3*dx
678  else if (dx <= 400.) then
679  v_column=4.995e-3*dx**.807
680  v_bullet=3.197e-3*dx**.902
681  v_plate=1.48e-3*dx**.926
682  else if (dx <= 600.) then
683  v_column=2.223e-2*dx**.558
684  v_bullet=2.977e-2*dx**.529
685  v_plate=9.5e-4*dx
686  else if (dx <= 800.) then
687  v_column=4.352e-2*dx**.453
688  v_bullet=2.144e-2*dx**.581
689  v_plate=3.161e-3*dx**.812
690  else
691  v_column=3.833e-2*dx**.472
692  v_bullet=3.948e-2*dx**.489
693  v_plate=7.109e-3*dx**.691
694  endif
695 !
696 !---- Reduce fall speeds from 400 mb to 1000 mb
697 !
698  v_column=px*v_column
699  v_bullet=px*v_bullet
700  v_plate=px*v_plate
701 !
702 !---- DIFFERENT VERSION! CALCULATES MASS-WEIGHTED CRYSTAL FALL SPEEDS
703 !
704  vel(i)=(m_bullet*v_bullet+m_column*v_column+m_plate*v_plate)/ &
705  & mass(i)
706  mass(i)=mass(i)/3.
707 !
708 !---- Shape factor and characteristic length of various ice habits,
709 ! capacitance is equal to 4*PI*(Shape factor)
710 ! See Young (1993, pp. 143-152 for guidance)
711 !
712 !---- Bullets:
713 !
714 !---- Shape factor for bullets (Heymsfield, 1975)
715  c_bullet=.5*d
716 !---- Length-width functions for bullets from Heymsfield (JAS, 1972)
717  if (d > 0.3) then
718  wd=.25*d**.7856 ! Width (mm); a-axis
719  else
720  wd=.185*d**.552
721  endif
722 !---- Characteristic length for bullets (see first multiplicative term on right
723 ! side of eq. 7 multiplied by crystal width on p. 821 of Heymsfield, 1975)
724  cl_bullet=.5*pi*wd*(.25*wd+d)/(d+wd)
725 !
726 !---- Plates:
727 !
728 !---- Length-width function for plates from Heymsfield (JAS, 1972)
729  wd=.0449*d**.449 ! Width or thickness (mm); c-axis
730 !---- Eccentricity & shape factor for thick plates following Young (1993, p. 144)
731  ecc_plate=sqrt(1.-wd*wd/(d*d)) ! Eccentricity
732  c_plate=d*ecc_plate/asin(ecc_plate) ! Shape factor
733 !---- Characteristic length for plates following Young (1993, p. 150, eq. 6.6)
734  cl_plate=d+2.*wd ! Characteristic lengths for plates
735 !
736 !---- Columns:
737 !
738 !---- Length-width function for columns from Heymsfield (JAS, 1972)
739  if (d > 0.2) then
740  wd=.1973*d**.414 ! Width (mm); a-axis
741  else
742  wd=.5*d ! Width (mm); a-axis
743  endif
744 !---- Eccentricity & shape factor for columns following Young (1993, p. 144)
745  ecc_column=sqrt(1.-wd*wd/(d*d)) ! Eccentricity
746  c_column=ecc_column*d/log((1.+ecc_column)*d/wd) ! Shape factor
747 !---- Characteristic length for columns following Young (1993, p. 150, eq. 6.7)
748  cl_column=(wd+2.*d)/(c1+c2*d/wd) ! Characteristic lengths for columns
749 !
750 !---- Convert shape factor & characteristic lengths from mm to m for
751 ! ventilation calculations
752 !
753  c_bullet=.001*c_bullet
754  c_plate=.001*c_plate
755  c_column=.001*c_column
756  cl_bullet=.001*cl_bullet
757  cl_plate=.001*cl_plate
758  cl_column=.001*cl_column
759 !
760 !---- Make a smooth transition between a ventilation coefficient of 1.0 at 50 microns
761 ! to 1.1 at 200 microns
762 !
763  if (d > 0.2) then
764  cvent1=cvent1i
765  cvent2=cvent2i/3.
766  else
767  cvent1=1.0+.1*max(0., d-.05)/.15
768  cvent2=0.
769  endif
770 !
771 !---- Ventilation factors for ice crystals:
772 !
773  vent1(i)=cvent1*(c_bullet+c_plate+c_column)/3.
774  vent2(i)=cvent2*(c_bullet*sqrt(v_bullet*cl_bullet) &
775  & +c_plate*sqrt(v_plate*cl_plate) &
776  & +c_column*sqrt(v_column*cl_column) )
777  crime_best=crime3 ! For calculating Best No. of rimed ice crystals
778  else
779 !
780 !---- This block of code calculates bulk characteristics based on average
781 ! characteristics of unrimed aggregates >= 1.5 mm using Locatelli &
782 ! Hobbs (JGR, 1974, 2185-2197) data.
783 !
784 !----- This category is a composite of aggregates of unrimed radiating
785 !----- assemblages of dendrites or dendrites; aggregates of unrimed
786 !----- radiating assemblages of plates, side planes, bullets, & columns;
787 !----- aggregates of unrimed side planes (mass in mg, velocity in m/s)
788 !
789  m_agg=(.073*d**1.4+.037*d**1.9+.04*d**1.4)/3.
790  v_agg=(.8*d**.16+.69*d**.41+.82*d**.12)/3.
791  mass(i)=m_agg
792  vel(i)=v_agg
793 !
794 !---- Assume spherical aggregates
795 !
796 !---- Shape factor is the same as for bullets, = D/2
797  c_agg=.001*.5*d ! Units of m
798 !---- Characteristic length is surface area divided by perimeter
799 ! (.25*PI*D**2)/(PI*D**2) = D/4
800  cl_agg=.5*c_agg ! Units of m
801 !
802 !---- Ventilation factors for aggregates:
803 !
804  vent1(i)=cvent1a*c_agg
805  vent2(i)=cvent2a*c_agg*sqrt(v_agg*cl_agg)
806  crime_best=crime4 ! For calculating Best No. of rimed aggregates
807  endif
808 !
809 !---- Convert from shape factor to capacitance for ventilation factors
810 !
811  vent1(i)=4.*pi*vent1(i)
812  vent2(i)=4.*pi*vent2(i)
813  diam(i)=1.e-3*d ! Convert from mm to m
814  mass(i)=1.e-6*mass(i) ! Convert from mg to kg
815 !
816 !---- Calculate increase in fall speeds of individual rimed ice particles
817 !
818  do k=0,nrime
819 !---- Mass of rimed ice particle associated with rime_factor(k)
820  rime_m1=rime_factor(k)*mass(i)
821  rime_m2=cbulk_ice*diam(i)**3
822  m_rime=min(rime_m1, rime_m2)
823 !---- Best Number (X) of rimed ice particle combining eqs. (8) & (12) in Bohm
824  x_rime=crime2*m_rime*(crime_best/m_rime)**.25
825 !---- Reynolds Number for rimed ice particle using eq. (11) in Bohm
826  re_rime=8.5*(sqrt(1.+.1519*sqrt(x_rime))-1.)**2
827  rime_vel(k)=crime5*re_rime/diam(i)
828  enddo
829  do k=1,nrime
830  vel_rime(i,k)=rime_vel(k)/rime_vel(0)
831  enddo
832  if (print_diag) then
833  !
834  !---- Determine if statistics should be printed out.
835  !
836  iprint=.false.
837  if (d <= 1.) then
838  if (mod(i,10) == 0) iprint=.true.
839  else
840  if (mod(i,100) == 0) iprint=.true.
841  endif
842  if (iprint) then
843  write(6,"(f7.4,5e11.4,1x,5f7.4,1x,2e11.4)") &
844  & d,1.e6*mass(i),m_bullet,m_column,m_plate,m_agg, &
845  & vel(i),v_bullet,v_column,v_plate,v_agg, &
846  & vent1(i),vent2(i)
847  write(7,"(f7.4,23f5.2)") d,(vel_rime(i,k), k=1,5), &
848  & (vel_rime(i,k), k=6,40,2)
849  icount=icount+1
850  endif
851  endif
852  enddo
853 !
854 !#######################################################################
855 ! Characteristics as functions of mean particle diameter
856 !#######################################################################
857 !
858  venti1=0.
859  venti2=0.
860  accri=0.
861  massi=0.
862  vsnowi=0.
863  vel_rf=0.
864  ivel_rime=0.
865  icount=0
866  if (print_diag) then
867  icount=60
868  write(6,"(/2a)") '------------- Statistics as functions of ', &
869  & 'mean particle diameter -------------'
870  write(7,"(/2a)") '------ Increase in fall speeds of rimed ice', &
871  & ' particles as functions of mean particle diameter -----'
872  endif
873  do j=mdimin,mdimax
874  if (icount == 60 .AND. print_diag) then
875  write(6,"(/2a)") 'D(mm) Vent1 Vent2 ', &
876  & 'Accrete Mass Vel Dens '
877  write(7,"(3a)") ' <----------------------------------', &
878  & '--------------- Rime Factor --------------------------', &
879  & '--------------------------->'
880  write(7,"(a,23f5.2)") 'D(mm)',(rime_factor(k), k=1,5), &
881  & (rime_factor(k), k=6,40,2)
882  icount=0
883  endif
884  mdiam=deldmi*float(j) ! in m
885  smom3=0.
886  pratei=0.
887  rime_vel=0. ! Note that this array is being reused!
888  do i=idimin,idimax
889  dx=diam(i)/mdiam
890  if (dx <= xmax) then ! To prevent arithmetic underflows
891  expf=exp(-dx)*ddeli
892  venti1(j)=venti1(j)+vent1(i)*expf
893  venti2(j)=venti2(j)+vent2(i)*expf
894  accri(j)=accri(j)+diam(i)*diam(i)*vel(i)*expf
895  xmass=mass(i)*expf
896  do k=1,nrime
897  rime_vel(k)=rime_vel(k)+xmass*vel_rime(i,k)
898  enddo
899  massi(j)=massi(j)+xmass
900  pratei=pratei+xmass*vel(i)
901  smom3=smom3+diam(i)**3*expf
902  else
903  exit
904  endif
905  enddo
906  !
907  !--- Increased fall velocities functions of mean diameter (j),
908  ! normalized by ice content, and rime factor (k)
909  !
910  do k=1,nrime
911  ivel_rime(j,k)=rime_vel(k)/massi(j)
912  enddo
913  !
914  !--- Increased fall velocities functions of ice content at 0.1 mm
915  ! intervals (j_100) and rime factor (k); accumulations here
916  !
917  jj=j/100
918  if (jj >= 2 .AND. jj <= 9) then
919  do k=1,nrime
920  vel_rf(jj,k)=vel_rf(jj,k)+ivel_rime(j,k)
921  enddo
922  endif
923  bulk_dens=cbulk*massi(j)/smom3
924  venti1(j)=venti1(j)/mdiam
925  venti2(j)=venti2(j)/mdiam
926  accri(j)=accri(j)/mdiam
927  vsnowi(j)=pratei/massi(j)
928  massi(j)=massi(j)/mdiam
929  if (mod(j,10) == 0 .AND. print_diag) then
930  xmdiam=1.e3*mdiam
931  write(6,"(f6.3,4e11.4,f6.3,f8.3)") xmdiam,venti1(j),venti2(j),&
932  & accri(j),massi(j),vsnowi(j),bulk_dens
933  write(7,"(f6.3,23f5.2)") xmdiam,(ivel_rime(j,k), k=1,5), &
934  & (ivel_rime(j,k), k=6,40,2)
935  icount=icount+1
936  endif
937  enddo
938 !
939 !--- Average increase in fall velocities rimed ice as functions of mean
940 ! particle diameter (j, only need 0.1 mm intervals) and rime factor (k)
941 !
942  if (print_diag) then
943  write(7,"(/2a)") ' ------- Increase in fall speeds of rimed ', &
944  & 'ice particles at reduced, 0.1-mm intervals --------'
945  write(7,"(3a)") ' <----------------------------------', &
946  & '--------------- Rime Factor --------------------------', &
947  & '--------------------------->'
948  write(7,"(a,23f5.2)") 'D(mm)',(rime_factor(k), k=1,5), &
949  & (rime_factor(k), k=6,40,2)
950  endif
951  do j=2,9
952  vel_rf(j,0)=1.
953  do k=1,nrime
954  vel_rf(j,k)=.01*vel_rf(j,k)
955  enddo
956  if (print_diag) write(7,"(f4.1,2x,23f5.2)") 0.1*j, &
957  & (vel_rf(j,k), k=1,5),(vel_rf(j,k), k=6,40,2)
958  enddo
959 !
960 !-----------------------------------------------------------------------------------
961 !

Here is the caller graph for this function:

subroutine module_microphysics::my_growth_rates ( real  DTPH)

Definition at line 392 of file module_bfmicrophysics.f.

References my_growth.

Referenced by gsmconst().

392 !
393  implicit none
394 !
395 !--- Below are tabulated values for the predicted mass of ice crystals
396 ! after 600 s of growth in water saturated conditions, based on
397 ! calculations from Miller and Young (JAS, 1979). These values are
398 ! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
399 ! Young (1993). Values at temperatures colder than -27C were
400 ! assumed to be invariant with temperature.
401 !
402 !--- Used to normalize Miller & Young (1979) calculations of ice growth
403 ! over large time steps using their tabulated values at 600 s.
404 ! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
405 !
406  real dtph, dt_ice
407  REAL my_600(my_t1:my_t2)
408 !
409 !-- 20090714: These values are in g and need to be converted to kg below
410  DATA my_600 / &
411  & 5.5e-8, 1.4e-7, 2.8e-7, 6.e-7, 3.3e-6, & ! -1 to -5 deg C
412  & 2.e-6, 9.e-7, 8.8e-7, 8.2e-7, 9.4e-7, & ! -6 to -10 deg C
413  & 1.2e-6, 1.85e-6, 5.5e-6, 1.5e-5, 1.7e-5, & ! -11 to -15 deg C
414  & 1.5e-5, 1.e-5, 3.4e-6, 1.85e-6, 1.35e-6, & ! -16 to -20 deg C
415  & 1.05e-6, 1.e-6, 9.5e-7, 9.0e-7 , 9.5e-7, & ! -21 to -25 deg C
416  & 9.5e-7, 9.e-7, 9.e-7, 9.e-7, 9.e-7, & ! -26 to -30 deg C
417  & 9.e-7, 9.e-7, 9.e-7, 9.e-7, 9.e-7 / ! -31 to -35 deg C
418 !
419 !-----------------------------------------------------------------------
420 !
421  dt_ice=(dtph/600.)**1.5
422 ! MY_GROWTH=DT_ICE*MY_600 ! original version
423  my_growth=dt_ice*my_600*1.e-3 !-- 20090714: Convert from g to kg
424 !
425 !-----------------------------------------------------------------------
426 !

Here is the caller graph for this function:

subroutine module_microphysics::rain_lookup ( )

Definition at line 969 of file module_bfmicrophysics.f.

References mdrmax, and mdrmin.

Referenced by gsmconst().

969  implicit none
970 !
971 !--- Parameters & arrays for fall speeds of rain as a function of rain drop
972 ! diameter. These quantities are integrated over exponential size
973 ! distributions of rain drops at 1 micron intervals (DdelR) from minimum
974 ! drop sizes of .05 mm (50 microns, DminR) to maximum drop sizes of 10 mm
975 ! (DmaxR).
976 !
977  real, parameter :: dminr=.05e-3, dmaxr=10.e-3, ddelr=1.e-6, &
978  & xrmin=1.e6*dminr, xrmax=1.e6*dmaxr
979  integer, parameter :: idrmin=xrmin, idrmax=xrmax
980  real diam(idrmin:idrmax), vel(idrmin:idrmax)
981 !
982 !--- Parameters rain lookup tables, which establish the range of mean drop
983 ! diameters; from a minimum mean diameter of 0.05 mm (DMRmin) to a
984 ! maximum mean diameter of 0.45 mm (DMRmax). The tables store solutions
985 ! at 1 micron intervals (DelDMR) of mean drop diameter.
986 !
987  real mdiam, mass
988 !
989  logical, parameter :: print_diag=.false.
990 !
991  real d, cmass, pi2, expf
992  integer i, j, i1, i2
993 !
994 !-----------------------------------------------------------------------
995 !------- Fall speeds of rain as function of rain drop diameter ---------
996 !-----------------------------------------------------------------------
997 !
998  do i=idrmin,idrmax
999  diam(i)=float(i)*ddelr
1000  d=100.*diam(i) ! Diameter in cm
1001  if (d <= .42) then
1002  !
1003  !--- Rutledge & Hobbs (1983); vel (m/s), d (cm)
1004  !
1005  vel(i)=max(0., -.267+51.5*d-102.25*d*d+75.7*d**3)
1006  else if (d > 0.42 .and. d <= .58) then
1007  !
1008  !--- Linear interpolation of Gunn & Kinzer (1949) data
1009  !
1010  vel(i)=8.92+.25/(.58-.42)*(d-.42)
1011  else
1012  vel(i)=9.17
1013  endif
1014  enddo
1015  do i=1,100
1016  i1=(i-1)*100+idrmin
1017  i2=i1+90
1018  !
1019  !--- Print out rain fall speeds only for D<=5.8 mm (.58 cm)
1020  !
1021  if (diam(i1) > .58e-2) exit
1022  if (print_diag) then
1023  write(6,"(1x)")
1024  write(6,"('D(mm)-> ',10f7.3)") (1000.*diam(j), j=i1,i2,10)
1025  write(6,"('V(m/s)-> ',10f7.3)") (vel(j), j=i1,i2,10)
1026  endif
1027  enddo
1028 !
1029 !-----------------------------------------------------------------------
1030 !------------------- Lookup tables for rain processes ------------------
1031 !-----------------------------------------------------------------------
1032 !
1033 ! pi=acos(-1.)
1034  pi2=2.*pi
1035  cmass=1000.*pi/6.
1036  if (print_diag) then
1037  write(6,"(/'Diam - Mean diameter (mm)' &
1038  & /'VENTR1 - 1st ventilation coefficient (m**2)' &
1039  & /'VENTR2 - 2nd ventilation coefficient (m**3/s**.5)' &
1040  & /'ACCRR - accretion moment (m**4/s)' &
1041  & /'RHO*QR - mass content (kg/m**3) for N0r=8e6' &
1042  & /'RRATE - rain rate moment (m**5/s)' &
1043  & /'VR - mass-weighted rain fall speed (m/s)' &
1044  & /' Diam VENTR1 VENTR2 ACCRR ', &
1045  & 'RHO*QR RRATE VR')")
1046  endif
1047  do j=mdrmin,mdrmax
1048  mdiam=float(j)*deldmr
1049  ventr2(j)=0.
1050  accrr(j)=0.
1051  massr(j)=0.
1052  rrate(j)=0.
1053  do i=idrmin,idrmax
1054  expf=exp(-diam(i)/mdiam)*ddelr
1055  ventr2(j)=ventr2(j)+diam(i)**1.5*vel(i)**.5*expf
1056  accrr(j)=accrr(j)+diam(i)*diam(i)*vel(i)*expf
1057  massr(j)=massr(j)+diam(i)**3*expf
1058  rrate(j)=rrate(j)+diam(i)**3*vel(i)*expf
1059  enddo
1060  !
1061  !--- Derived based on ventilation, F(D)=0.78+.31*Schmidt**(1/3)*Reynold**.5,
1062  ! where Reynold=(V*D*rho/dyn_vis), V is velocity, D is particle diameter,
1063  ! rho is air density, & dyn_vis is dynamic viscosity. Only terms
1064  ! containing velocity & diameter are retained in these tables.
1065  !
1066  ventr1(j)=.78*pi2*mdiam**2
1067  ventr2(j)=.31*pi2*ventr2(j)
1068  !
1069  massr(j)=cmass*massr(j)
1070  rrate(j)=cmass*rrate(j)
1071  vrain(j)=rrate(j)/massr(j)
1072  if (print_diag) write(6,"(f6.3,5g12.5,f6.3)") 1000.*mdiam, &
1073  & ventr1(j),ventr2(j),accrr(j),8.e6*massr(j),rrate(j),vrain(j)
1074  enddo
1075 !
1076 !-----------------------------------------------------------------------
1077 !

Here is the caller graph for this function:

subroutine module_microphysics::rsipath ( integer  im,
integer  ix,
integer  ix2,
integer  levs,
real, dimension(ix,levs)  prsl,
real, dimension(ix,levs+1)  prsi,
real, dimension(ix,levs)  t,
real, dimension(ix,levs)  q,
  clw 
)

Definition at line 2634 of file module_bfmicrophysics.f.

References indexsmin, mdimax, mdimin, rerainmin, thom, and tnw.

2634  &, f_ice, f_rain, f_rime, flgmin &
2635  &, cwatp, cicep, rainp, snowp &
2636  &, recwat, rerain, resnow, lprnt, ipr)
2637 !
2638  implicit none
2639 !
2640 !--------------------CLOUD----------------------------------------------
2641  integer im, ix, ix2, levs, ipr
2642  real prsl(ix,levs), prsi(ix,levs+1), t(ix,levs), q(ix,levs) &
2643  &, clw(ix2,levs), f_ice(ix2,levs), f_rain(ix2,levs) &
2644  &, f_rime(ix2,levs) &
2645  &, cwatp(ix,levs), rainp(ix,levs), cicep(ix,levs) &
2646  &, snowp(ix,levs), recwat(ix,levs), resnow(ix,levs) &
2647  &, rerain(ix,levs)
2648  real flgmin
2649  real frice, frrain, qcice, qcwat, qrain, qsnow, qtot, sden &
2650  &, cpath, rho, dsnow, flarge, rimef, xsimass, nlice &
2651  &, tc, recw1, drain, xli, dum, NLImax, pfac, pp &
2652  &, snofac, tem
2653 !
2654  real, parameter :: cexp=1./3.
2655  integer i, l, indexs
2656  logical lprnt
2657 !
2658 
2659  recw1 = 620.3505 / tnw**cexp ! cloud droplet effective radius
2660 
2661  do l=1,levs
2662  do i=1,im
2663  !--- HYDROMETEOR'S OPTICAL PATH
2664  cwatp(i,l) = 0.
2665  cicep(i,l) = 0.
2666  rainp(i,l) = 0.
2667  snowp(i,l) = 0.
2668  !--- HYDROMETEOR'S EFFECTIVE RADIUS
2669  recwat(i,l) = recwmin
2670  rerain(i,l) = rerainmin
2671  resnow(i,l) = resnowmin
2672  ENDDO
2673  ENDDO
2674 
2675  do l=1,levs
2676  DO i=1,im
2677 
2678  ! Assume parameterized condensate is
2679  ! all water for T>=-10C,
2680  ! all ice for T<=-30C,
2681  ! and a linear mixture at -10C > T > -30C
2682  !
2683  ! * Determine hydrometeor composition of total condensate (QTOT)
2684  !
2685 ! pp = prsl(i,l) * 1000.0
2686  pp = prsl(i,l) / prsi(i,levs+1)
2687 ! pfac = max(0.25, sqrt(sqrt(min(1.0, pp*0.000025))))
2688 ! pfac = max(0.5, sqrt(sqrt(min(1.0, pp*0.000025))))
2689 ! pfac = max(0.5, sqrt(sqrt(min(1.0, pp*0.00002))))
2690 ! pfac = max(0.25, sqrt(sqrt(min(1.0, pp*0.00001))))
2691 ! pfac = max(0.25, sqrt(sqrt(min(1.0, pp))))
2692 ! pfac = max(0.1, sqrt(min(1.0, pp*0.00001)))
2693 ! pfac = max(0.5, sqrt(sqrt(min(1.0, pp*0.000033))))
2694 ! pfac = max(0.5, sqrt(sqrt(min(1.0, pp*0.00004))))
2695 !go pfac = max(0.5, (sqrt(min(1.0, pp*0.000025))))
2696  pfac = 1.0
2697  tc = t(i,l) - t0c
2698  qtot = clw(i,l)
2699  IF (qtot > epsq) THEN
2700  qcwat=0.
2701  qcice=0.
2702  qrain=0.
2703  qsnow=0.
2704  frice = max(0.0, min(1.0, f_ice(i,l)))
2705  frrain = max(0.0, min(1.0, f_rain(i,l)))
2706  IF(tc <= thom) then
2707  qcice = qtot
2708  ELSE
2709  qcice = frice * qtot
2710  qcwat = qtot - qcice
2711  qrain = frrain * qcwat
2712  qcwat = qcwat - qrain
2713  ENDIF
2714  !
2715  !--- Air density (RHO), model mass thickness (CPATH)
2716  !
2717  rho = prsl(i,l)/(rd*t(i,l)*(1.+eps1*q(i,l)))
2718  cpath = (prsi(i,l+1)-prsi(i,l))*(1000000.0/grav)
2719 
2720  !! CLOUD WATER
2721  !
2722  !--- Effective radius (RECWAT) & total water path (CWATP)
2723  ! Assume monodisperse distribution of droplets (no factor of 1.5)
2724  !
2725  IF(qcwat > 0.) THEN
2726  recwat(i,l) = max(recwmin, recw1*(rho*qcwat)**cexp)
2727  cwatp(i,l) = cpath*qcwat ! cloud water path
2728 ! tem = 5.0*(1+max(0.0,min(1.0,-0.05*tc)))
2729 ! RECWAT(I,L) = max(RECWAT(I,L), tem)
2730  ENDIF
2731 
2732  !! RAIN
2733  !
2734  !--- Effective radius (RERAIN) & total water path (RAINP)
2735  !--- Factor of 1.5 accounts for r**3/r**2 moments for exponentially
2736  ! distributed drops in effective radius calculations
2737  ! (from M.D. Chou's code provided to Y.-T. Hou)
2738  !
2739  IF(qrain > 0.) THEN
2740  drain = cn0r0*sqrt(sqrt((rho*qrain)))
2741  rerain(i,l) = 1.5*max(xmrmin, min(drain, xmrmax))
2742  rainp(i,l) = cpath*qrain ! rain water path
2743  ENDIF
2744 
2745  !! SNOW (large ice) & CLOUD ICE
2746  !
2747  !--- Effective radius (RESNOW) & total ice path (SNOWP)
2748  !--- Total ice path (CICEP) for cloud ice
2749  !--- Factor of 1.5 accounts for r**3/r**2 moments for exponentially
2750  ! distributed ice particles in effective radius calculations
2751  !
2752  !--- Separation of cloud ice & "snow" uses algorithm from
2753  ! subroutine GSMCOLUMN
2754  !
2755  IF(qcice > 0.) THEN
2756  !
2757  !--- Mean particle size following Houze et al. (JAS, 1979, p. 160),
2758  ! converted from Fig. 5 plot of LAMDAs. An analogous set of
2759  ! relationships also shown by Fig. 8 of Ryan (BAMS, 1996, p. 66),
2760  ! but with a variety of different relationships that parallel the
2761  ! Houze curves.
2762  !
2763 ! DUM=MAX(0.05, MIN(1., EXP(.0536*TC)) )
2764  dum=max(0.05, min(1., exp(.0564*tc)) )
2765  indexs=min(mdimax, max(mdimin, int(xmimax*dum) ) )
2766 ! indexs=max(INDEXSmin, indexs)
2767 ! NLImax=5.E3/sqrt(DUM) !- Ver3
2768  dum=max(flgmin*pfac, dum)
2769 ! DUM=MAX(FLGmin, DUM)
2770 ! NLImax=20.E3 !- Ver3
2771 ! NLImax=50.E3 !- Ver3 => comment this line out
2772  nlimax=10.e3/sqrt(dum) !- Ver3
2773 ! NLImax=5.E3/sqrt(DUM) !- Ver3
2774 ! NLImax=6.E3/sqrt(DUM) !- Ver3
2775 ! NLImax=7.5E3/sqrt(DUM) !- Ver3
2776 ! NLImax=20.E3/DUM !- Ver3
2777 ! NLImax=20.E3/max(0.2,DUM) !- Ver3
2778 ! NLImax=2.0E3/max(0.1,DUM) !- Ver3
2779 ! NLImax=2.5E3/max(0.1,DUM) !- Ver3
2780 ! NLImax=10.E3/max(0.2,DUM) !- Ver3
2781 ! NLImax=4.E3/max(0.2,DUM) !- Ver3
2782 !Moorthi DSNOW = XMImax*EXP(.0536*TC)
2783 !Moorthi INDEXS = MAX(INDEXSmin, MIN(MDImax, INT(DSNOW)))
2784  !
2785  !--- Assumed number fraction of large ice to total (large & small)
2786  ! ice particles, which is based on a general impression of the
2787  ! literature.
2788  !
2789  ! Small ice are assumed to have a mean diameter of 50 microns.
2790  !
2791  IF(tc >= 0.) THEN
2792  flarge=flg1p0
2793  ELSE
2794  flarge = dum
2795  ENDIF
2796 !------------------------Commented by Moorthi -----------------------------
2797 ! ELSEIF (TC >= -25.) THEN
2798 !
2799 !--- Note that absence of cloud water (QCWAT) is used as a quick
2800 ! substitute for calculating water subsaturation as in GSMCOLUMN
2801 !
2802 ! IF(QCWAT <= 0. .OR. TC < -8.
2803 ! & .OR. TC > -3.)THEN
2804 ! FLARGE=FLG0P2
2805 ! ELSE
2806 
2807 !--- Parameterize effects of rime splintering by increasing
2808 ! number of small ice particles
2809 !
2810 ! FLARGE=FLG0P1
2811 ! ENDIF
2812 ! ELSEIF (TC <= -50.) THEN
2813 ! FLARGE=.01
2814 ! ELSE
2815 ! FLARGE=.2*EXP(.1198*(TC+25))
2816 ! ENDIF
2817 !____________________________________________________________________________
2818 
2819  rimef=max(1., f_rime(i,l))
2820  xsimass=massi(mdimin)*(1.-flarge)/flarge
2821 ! if (lprnt) print *,' rimef=',rimef,' xsimass=',xsimass
2822 ! &,' indexs=',indexs,' massi=',massi(indexs),' flarge=',flarge
2823  nlice=rho*qcice/(xsimass+rimef*massi(indexs))
2824  !
2825  !--- From subroutine GSMCOLUMN:
2826  !--- Minimum number concentration for large ice of NLImin=10/m**3
2827  ! at T>=0C. Done in order to prevent unrealistically small
2828  ! melting rates and tiny amounts of snow from falling to
2829  ! unrealistically warm temperatures.
2830  !
2831  IF(tc >= 0.) THEN
2832  nlice=max(nlimin, nlice)
2833  ELSEIF (nlice > nlimax) THEN
2834  !
2835  !--- Ferrier 6/13/01: Prevent excess accumulation of ice
2836  !
2837  xli=(rho*qcice/nlimax-xsimass)/rimef
2838 
2839  IF(xli <= massi(450) ) THEN
2840  dsnow=9.5885e5*xli**.42066
2841  ELSE
2842  dsnow=3.9751e6*xli**.49870
2843  ENDIF
2844 
2845  indexs=min(mdimax, max(indexs, int(dsnow)))
2846  nlice=rho*qcice/(xsimass+rimef*massi(indexs))
2847  ENDIF
2848 
2849 ! if (tc > -20.0 .and. indexs >= indexsmin) then
2850 ! snofac = max(0.0, min(1.0, exp(1.0*(tc+20.0))))
2851 ! if (indexs >= indexsmin) then
2852 ! if (tc > -20.0 .or. indexs >= indexsmin) then
2853 ! if (tc > -40.0) then
2854 ! if (tc >= -40.0 .or. prsl(i,l) > 50.0) then
2855 !! if (tc >= -20.0) then
2856 ! if (tc >= -20.0 .or. prsl(i,l) > 50.0) then
2857 ! if ((tc >= -20.0 .or.
2858 ! & prsi(i,levs+1)-prsi(i,l) < 30.0)
2859  if (prsi(i,levs+1)-prsi(i,l) < 40.0 &
2860 ! if (prsi(i,levs+1)-prsi(i,l) < 70.0
2861  & .and. indexs >= indexsmin) then
2862 ! & prsi(i,levs)-prsl(i,l) < 20.0) then
2863 ! & prsi(i,levs)-prsl(i,l) < 30.0) then
2864 ! & prsi(i,levs)-prsl(i,l) < 40.0) then
2865 ! snofac = max(0.0, min(1.0, 0.05*(tc+40.0)))
2866 ! snofac = max(0.0, min(1.0, 0.1*(tc+25.0)))
2867 ! snofac = max(0.0, min(1.0, 0.0667*(tc+25.0)))
2868 ! if (indexs > indexsmin) then
2869  qsnow = min(qcice, nlice*rimef*massi(indexs)/rho)
2870 ! & * snofac
2871  endif
2872 ! qsnow = qcice
2873  qcice = max(0., qcice-qsnow)
2874 ! qsnow = 0.0
2875  cicep(i,l) = cpath*qcice ! cloud ice path
2876  resnow(i,l) = 1.5*float(indexs)
2877  sden = sdens(indexs)/rimef ! 1/snow density
2878  snowp(i,l) = cpath*qsnow*sden ! snow path / snow density
2879 ! SNOWP (I,L) = CPATH*QSNOW ! snow path / snow density
2880 ! if (lprnt .and. i == ipr) then
2881 ! print *,' L=',L,' snowp=',snowp(i,l),' cpath=',cpath
2882 ! &,' qsnow=',qsnow,' sden=',sden,' rimef=',rimef,' indexs=',indexs
2883 ! &,' sdens=',sdens(indexs),' resnow=',resnow(i,l)
2884 ! &,' qcice=',qcice,' cicep=',cicep(i,l)
2885 ! endif
2886 
2887 
2888  ENDIF ! END QCICE BLOCK
2889  ENDIF ! QTOT IF BLOCK
2890 
2891  ENDDO
2892  ENDDO
2893 !
subroutine module_microphysics::rsipath2 ( )
Parameters
[in]plyrreal, (IM,LEVS), model layer mean pressure in mb (100Pa)
[in]plvlreal, (IM,LEVS+1), model level pressure in mb (100Pa)
[in]tlyrreal, (IM,LEVS), model layer mean temperature in K
[in]qlyrreal, (IM,LEVS), layer specific humidity in gm/gm
[in]qcwatreal, (IM,LEVS), layer cloud liquid water condensate amount
[in]qcicereal, (IM,LEVS), layer cloud ice water condensate amount
[in]qrainreal, (IM,LEVS), layer rain drop water amount
[in]rrimereal, (IM,LEVS), mass ratio of total to unrimed ice ( >= 1 )
[in]IMinteger, horizontal dimention
[in]LEVSinteger, vertical layer dimensions
[in]iflipinteger, control flag for in/out vertical indexing
=0: index from toa to surface
=1: index from surface to toa
[in]flgminreal, minimum large ice fraction
[in]lprntlogical, logical check print control flag
[out]cwatpreal, (IM,LEVS), layer cloud liquid water path
[out]cicepreal, (IM,LEVS), layer cloud ice water path
[out]rainpreal, (IM,LEVS), layer rain water path
[out]snowpreal, (IM,LEVS), layer snow water path
[out]recwatreal, (IM,LEVS), layer cloud eff radius for liqid water (micron)
[out]rerainreal, (IM,LEVS), layer rain water effective radius (micron)
[out]resnowreal, (IM,LEVS), layer snow flake effective radius (micron)
[out]sndenreal, (IM,LEVS), 1/snow density

General Algorithm

\

  1. Set up pressure related arrays, convert unit from mb to cb (10Pa). The rest part uses cb in computation.
  2. Compute air density (rho), model mass thickness (cpath), and temperature in C (tc)
  3. For cloud water, compute effective raius (recwat) and total water path (cwatp): assume monodisperse distribution of droplets (no factor of 1.5)
  4. For rain water, compute effective radius (rerain) and total water path (rainp): factor of 1.5 accounts for \( r^3/r^2 \) moments for exponentially distributed drops in effective raius calculations (from m.d. chou's code provided to y.-t. hou)
  5. For cloud ice,
    - mean particle size following houze et al. (jas, 1979, p.160), converted from fig.5 plot of lamdas. An analogous set of relationships also shown by fig.8 of ryan (bams, 1996, p.66), but with a variaty of different relationships parallel the houze curves.

Definition at line 2927 of file module_bfmicrophysics.f.

References indexsmin, mdimax, mdimin, rerainmin, and tnw.

Referenced by module_radiation_clouds::progcld2().

2927 !...................................
2928 
2929 ! --- inputs:
2930  & ( plyr, plvl, tlyr, qlyr, qcwat, qcice, qrain, rrime, &
2931  & im, levs, iflip, flgmin, &
2932 ! --- outputs:
2933  & cwatp, cicep, rainp, snowp, recwat, rerain, resnow, snden &
2934  & )
2935 
2936 ! ================= subprogram documentation block ================ !
2937 ! !
2938 ! abstract: this program is a modified version of ferrier's original !
2939 ! "rsipath" subprogram. it computes layer's cloud liquid, ice, rain, !
2940 ! and snow water condensate path and the partical effective radius !
2941 ! for liquid droplet, rain drop, and snow flake. !
2942 ! !
2943 ! ==================== defination of variables ==================== !
2944 ! !
2945 ! input variables: !
2946 ! plyr (IM,LEVS) : model layer mean pressure in mb (100Pa) !
2947 ! plvl (IM,LEVS+1):model level pressure in mb (100Pa) !
2948 ! tlyr (IM,LEVS) : model layer mean temperature in k !
2949 ! qlyr (IM,LEVS) : layer specific humidity in gm/gm !
2950 ! qcwat (IM,LEVS) : layer cloud liquid water condensate amount !
2951 ! qcice (IM,LEVS) : layer cloud ice water condensate amount !
2952 ! qrain (IM,LEVS) : layer rain drop water amount !
2953 ! rrime (IM,LEVS) : mass ratio of total to unrimed ice ( >= 1 ) !
2954 ! IM : horizontal dimention !
2955 ! LEVS : vertical layer dimensions !
2956 ! iflip : control flag for in/out vertical indexing !
2957 ! =0: index from toa to surface !
2958 ! =1: index from surface to toa !
2959 ! flgmin : Minimum large ice fraction !
2960 ! lprnt : logical check print control flag !
2961 ! !
2962 ! output variables: !
2963 ! cwatp (IM,LEVS) : layer cloud liquid water path !
2964 ! cicep (IM,LEVS) : layer cloud ice water path !
2965 ! rainp (IM,LEVS) : layer rain water path !
2966 ! snowp (IM,LEVS) : layer snow water path !
2967 ! recwat(IM,LEVS) : layer cloud eff radius for liqid water (micron) !
2968 ! rerain(IM,LEVS) : layer rain water effective radius (micron) !
2969 ! resnow(IM,LEVS) : layer snow flake effective radius (micron) !
2970 ! snden (IM,LEVS) : 1/snow density !
2971 ! !
2972 ! !
2973 ! usage: call rsipath2 !
2974 ! !
2975 ! subroutines called: none !
2976 ! !
2977 ! program history log: !
2978 ! xx-xx-2001 b. ferrier - original program !
2979 ! xx-xx-2004 s. moorthi - modified for use in gfs model !
2980 ! 05-20-2004 y. hou - modified, added vertical index flag!
2981 ! to reduce data flipping, and rearrange code to !
2982 ! be comformable with radiation part programs. !
2983 ! !
2984 ! ==================== end of description ===================== !
2985 !
2986 
2987  implicit none
2988 
2989 ! --- constant parameter:
2990  real, parameter :: cexp= 1.0/3.0
2991 
2992 ! --- inputs:
2993  real, dimension(:,:), intent(in) :: &
2994  & plyr, plvl, tlyr, qlyr, qcwat, qcice, qrain, rrime
2995 
2996  integer, intent(in) :: im, levs, iflip
2997  real, dimension(:), intent(in) :: flgmin
2998 ! logical, intent(in) :: lprnt
2999 
3000 ! --- output:
3001  real, dimension(:,:), intent(out) :: &
3002  & cwatp, cicep, rainp, snowp, recwat, rerain, resnow, snden
3003 
3004 ! --- locals:
3005 ! real, dimension(IM,LEVS) :: delp, pp1, pp2
3006 
3007  real :: recw1, dsnow, qsnow, qqcice, flarge, xsimass, pfac, &
3008  & nlice, xli, nlimax, dum, tem, &
3009  & rho, cpath, rc, totcnd, tc
3010 
3011  integer :: i, k, indexs, ksfc, k1
3012 !
3013 !===> ... begin here
3014 !
3015  recw1 = 620.3505 / tnw**cexp ! cloud droplet effective radius
3016 
3017  do k = 1, levs
3018  do i = 1, im
3019  !--- hydrometeor's optical path
3020  cwatp(i,k) = 0.0
3021  cicep(i,k) = 0.0
3022  rainp(i,k) = 0.0
3023  snowp(i,k) = 0.0
3024  snden(i,k) = 0.0
3025  !--- hydrometeor's effective radius
3026  recwat(i,k) = recwmin
3027  rerain(i,k) = rerainmin
3028  resnow(i,k) = resnowmin
3029  enddo
3030  enddo
3031 
3034 ! --- set up pressure related arrays, convert unit from mb to cb (10Pa)
3035 ! cause the rest part uses cb in computation
3036 
3037  if (iflip == 0) then ! data from toa to sfc
3038  ksfc = levs + 1
3039  k1 = 0
3040  else ! data from sfc to top
3041  ksfc = 1
3042  k1 = 1
3043  endif ! end_if_iflip
3044 !
3045  do k = 1, levs
3046  do i = 1, im
3047  totcnd = qcwat(i,k) + qcice(i,k) + qrain(i,k)
3048  qsnow = 0.0
3049  if(totcnd > epsq) then
3051 ! --- air density (rho), model mass thickness (cpath), temperature in c (tc)
3052 
3053  rho = 0.1 * plyr(i,k) &
3054  & / (rd* tlyr(i,k) * (1.0 + eps1*qlyr(i,k)))
3055  cpath = abs(plvl(i,k+1) - plvl(i,k)) * (100000.0 / grav)
3056  tc = tlyr(i,k) - t0c
3057 
3060 ! cloud water
3061 !
3062 ! --- effective radius (recwat) & total water path (cwatp):
3063 ! assume monodisperse distribution of droplets (no factor of 1.5)
3064 
3065  if (qcwat(i,k) > 0.0) then
3066  recwat(i,k) = max(recwmin,recw1*(rho*qcwat(i,k))**cexp)
3067  cwatp(i,k) = cpath * qcwat(i,k) ! cloud water path
3068 ! tem = 5.0*(1.0 + max(0.0, min(1.0,-0.05*tc)))
3069 ! recwat(i,k) = max(recwat(i,k), tem)
3070  endif
3074 ! rain
3075 !
3076 ! --- effective radius (rerain) & total water path (rainp):
3077 ! factor of 1.5 accounts for r**3/r**2 moments for exponentially
3078 ! distributed drops in effective radius calculations
3079 ! (from m.d. chou's code provided to y.-t. hou)
3080 
3081  if (qrain(i,k) > 0.0) then
3082  tem = cn0r0 * sqrt(sqrt(rho*qrain(i,k)))
3083  rerain(i,k) = 1.5 * max(xmrmin, min(xmrmax, tem))
3084  rainp(i,k) = cpath * qrain(i,k) ! rain water path
3085  endif
3086 
3087 ! snow (large ice) & cloud ice
3088 !
3089 ! --- effective radius (resnow) & total ice path (snowp) for snow, and
3090 ! total ice path (cicep) for cloud ice:
3091 ! factor of 1.5 accounts for r**3/r**2 moments for exponentially
3092 ! distributed ice particles in effective radius calculations
3093 ! separation of cloud ice & "snow" uses algorithm from subroutine gsmcolumn
3094 
3095 ! pfac = max(0.5, sqrt(sqrt(min(1.0, pp1(i,k)*0.00004))))
3096 !go pfac = max(0.5, (sqrt(min(1.0, pp1(i,k)*0.000025))))
3097  pfac = 1.0
3098 
3099  if (qcice(i,k) > 0.0) then
3100 
3106 ! --- mean particle size following houze et al. (jas, 1979, p. 160),
3107 ! converted from fig. 5 plot of lamdas. an analogous set of
3108 ! relationships also shown by fig. 8 of ryan (bams, 1996, p. 66),
3109 ! but with a variety of different relationships that parallel
3110 ! the houze curves.
3111 
3112 ! dum = max(0.05, min(1.0, exp(0.0536*tc) ))
3113  dum = max(0.05, min(1.0, exp(0.0564*tc) ))
3114  indexs = min(mdimax, max(mdimin, int(xmimax*dum) ))
3115  dum=max(flgmin(i)*pfac, dum)
3116 
3117 !!\n - assumed number of fraction of large ice to total (large+small) ice
3118 !! particles, which is based on a general impression of the literature. small ice
3119 !! are assumed to have a mean diameter of 50 microns
3120 ! --- assumed number fraction of large ice to total (large & small) ice
3121 ! particles, which is based on a general impression of the literature.
3122 ! small ice are assumed to have a mean diameter of 50 microns.
3123 
3124  if (tc >= 0.0) then
3125  flarge = flg1p0
3126  else
3127  flarge = dum
3128 ! flarge = max(FLGmin*pfac, dum)
3129  endif
3130 !------------------------commented by moorthi -----------------------------
3131 ! elseif (tc >= -25.0) then
3132 !
3133 ! --- note that absence of cloud water (qcwat) is used as a quick
3134 ! substitute for calculating water subsaturation as in gsmcolumn
3135 !
3136 ! if (qcwat(i,k) <= 0.0 .or. tc < -8.0 &
3137 ! & .or. tc > -3.0) then
3138 ! flarge = FLG0P2
3139 ! else
3140 !
3141 ! --- parameterize effects of rime splintering by increasing
3142 ! number of small ice particles
3143 !
3144 ! flarge = FLG0P1
3145 ! endif
3146 ! elseif (tc <= -50.0) then
3147 ! flarge = 0.01
3148 ! else
3149 ! flarge = 0.2 * exp(0.1198*(tc+25.0))
3150 ! endif
3151 !____________________________________________________________________________
3152 
3153  xsimass = massi(mdimin) * (1.0 - flarge) / flarge
3154 ! nlimax = 20.0e3 !- ver3
3155 ! NLImax=50.E3 !- Ver3 => comment this line out
3156  nlimax=10.e3/sqrt(dum) !- Ver3
3157 ! NLImax=5.E3/sqrt(DUM) !- Ver3
3158 ! NLImax=6.E3/sqrt(DUM) !- Ver3
3159 ! NLImax=7.5E3/sqrt(DUM) !- Ver3
3160 
3161 ! indexs = min(MDImax, max(MDImin, int(XMImax*dum) ))
3162 !moorthi dsnow = XMImax * exp(0.0536*tc)
3163 !moorthi indexs = max(INDEXSmin, min(MDImax, int(dsnow)))
3164 
3165 ! if (lprnt) print *,' rrime=',rrime,' xsimass=',xsimass, &
3166 ! & ' indexs=',indexs,' massi=',massi(indexs),' flarge=',flarge
3167 
3168  tem = rho * qcice(i,k)
3169  nlice = tem / (xsimass +rrime(i,k)*massi(indexs))
3170 
3171 ! --- from subroutine gsmcolumn:
3172 ! minimum number concentration for large ice of NLImin=10/m**3
3173 ! at t>=0c. done in order to prevent unrealistically small
3174 ! melting rates and tiny amounts of snow from falling to
3175 ! unrealistically warm temperatures.
3176 
3177  if (tc >= 0.0) then
3178 
3179  nlice = max(nlimin, nlice)
3180 
3181  elseif (nlice > nlimax) then
3182 
3183 ! --- ferrier 6/13/01: prevent excess accumulation of ice
3184 
3185  xli = (tem/nlimax - xsimass) / rrime(i,k)
3186 
3187  if (xli <= massi(450) ) then
3188  dsnow = 9.5885e5 * xli**0.42066
3189  else
3190  dsnow = 3.9751e6 * xli** 0.49870
3191  endif
3192 
3193  indexs = min(mdimax, max(indexs, int(dsnow)))
3194  nlice = tem / (xsimass + rrime(i,k)*massi(indexs))
3195 
3196  endif ! end if_tc block
3197 
3198 ! if (abs(plvl(i,ksfc)-plvl(i,k+k1)) < 300.0 &
3199 ! if (abs(plvl(i,ksfc)-plvl(i,k+k1)) < 400.0 &
3200 ! if (plvl(i,k+k1) > 600.0 &
3201 ! & .and. indexs >= INDEXSmin) then
3202 ! if (tc > -20.0 .and. indexs >= indexsmin) then
3203  if (plvl(i,ksfc) > 850.0 .and. &
3204 ! & plvl(i,k+k1) > 600.0 .and. indexs >= indexsmin) then
3205  & plvl(i,k+k1) > 700.0 .and. indexs >= indexsmin) then ! 20060516
3206 !! if (plvl(i,ksfc) > 800.0 .and. &
3207 !! & plvl(i,k+k1) > 700.0 .and. indexs >= indexsmin) then
3208 ! if (plvl(i,ksfc) > 700.0 .and. &
3209 ! & plvl(i,k+k1) > 600.0 .and. indexs >= indexsmin) then
3210  qsnow = min( qcice(i,k), &
3211  & nlice*rrime(i,k)*massi(indexs)/rho )
3212  endif
3213 
3214  qqcice = max(0.0, qcice(i,k)-qsnow)
3215  cicep(i,k) = cpath * qqcice ! cloud ice path
3216  resnow(i,k) = 1.5 * float(indexs)
3217  snden(i,k) = sdens(indexs) / rrime(i,k) ! 1/snow density
3218  snowp(i,k) = cpath*qsnow ! snow path
3219 ! snowp (i,k) = cpath*qsnow*snden(i,k) ! snow path / snow density
3220 
3221 ! if (lprnt .and. i == ipr) then
3222 ! if (i == 2) then
3223 ! print *,' L=',k,' snowp=',snowp(i,k),' cpath=',cpath, &
3224 ! & ' qsnow=',qsnow,' sden=',snden(i,k),' rrime=',rrime(i,k),&
3225 ! & ' indexs=',indexs,' sdens=',sdens(indexs),' resnow=', &
3226 ! & resnow(i,k),' qcice=',qqcice,' cicep=',cicep(i,k)
3227 ! endif
3228 
3229  endif ! end if_qcice block
3230  endif ! end if_totcnd block
3231 
3232  enddo
3233  enddo
3234 !
3235 !...................................

Here is the caller graph for this function:

Variable Documentation

real, private module_microphysics::abfr
private

Definition at line 13 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), and gsmconst().

13  real,private :: abfr, cbfr, ciacw, ciacr, c_n0r0, &
real, private module_microphysics::c_n0r0
private

Definition at line 13 of file module_bfmicrophysics.f.

Referenced by gsmconst().

real, private module_microphysics::cbfr
private

Definition at line 13 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), and gsmconst().

real, private module_microphysics::ciacr
private

Definition at line 13 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), and gsmconst().

real, private module_microphysics::ciacw
private

Definition at line 13 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), and gsmconst().

real, parameter, private module_microphysics::dmimax =1.e-3
private

Definition at line 35 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

real, parameter, private module_microphysics::dmimin =.05e-3
private

Definition at line 35 of file module_bfmicrophysics.f.

35  REAL, PRIVATE,PARAMETER :: dmimin=.05e-3, dmimax=1.e-3, &
real, parameter, private module_microphysics::dmrmax =.45e-3
private

Definition at line 48 of file module_bfmicrophysics.f.

Referenced by gsmconst().

real, parameter, private module_microphysics::dmrmin =.05e-3
private

Definition at line 48 of file module_bfmicrophysics.f.

Referenced by gsmconst().

48  REAL, PRIVATE,PARAMETER :: dmrmin=.05e-3, dmrmax=.45e-3, &
integer, parameter, private module_microphysics::indexsmin =100
private

Definition at line 60 of file module_bfmicrophysics.f.

Referenced by rsipath(), and rsipath2().

60  INTEGER, PRIVATE, PARAMETER :: indexsmin=100
integer, parameter module_microphysics::ithi =40

Definition at line 87 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

integer, parameter module_microphysics::itlo =-60

Definition at line 87 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

87  INTEGER, PARAMETER :: itlo=-60, ithi=40
integer, parameter, private module_microphysics::mdimax =XMImax
private

Definition at line 38 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), gsmconst(), ice_lookup(), rsipath(), and rsipath2().

integer, parameter, private module_microphysics::mdimin =XMImin
private

Definition at line 38 of file module_bfmicrophysics.f.

Referenced by gsmconst(), ice_lookup(), rsipath(), and rsipath2().

38  INTEGER, PRIVATE,PARAMETER :: mdimin=xmimin, mdimax=xmimax
integer, parameter, private module_microphysics::mdrmax =XMRmax
private

Definition at line 52 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), gsmconst(), and rain_lookup().

integer, parameter, private module_microphysics::mdrmin =XMRmin
private

Definition at line 52 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), gsmconst(), and rain_lookup().

52  INTEGER, PRIVATE,PARAMETER :: mdrmin=xmrmin, mdrmax=xmrmax
integer, private module_microphysics::mic_step
private

Definition at line 18 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), and gsmconst().

18  integer, private :: mic_step
real, dimension(my_t1:my_t2), private module_microphysics::my_growth
private

Definition at line 28 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), and my_growth_rates().

28  REAL,PRIVATE,DIMENSION(MY_T1:MY_T2) :: my_growth
integer, parameter, private module_microphysics::my_t1 =1
private

Definition at line 27 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

27  INTEGER, PRIVATE,PARAMETER :: my_t1=1, my_t2=35
integer, parameter, private module_microphysics::my_t2 =35
private

Definition at line 27 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

integer, parameter, private module_microphysics::nrime =40
private

Definition at line 82 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), gsmconst(), and ice_lookup().

82  INTEGER, PRIVATE,PARAMETER :: nrime=40
integer, dimension(itlo:ithi,4) module_microphysics::nstats

Definition at line 88 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

88  INTEGER nstats(itlo:ithi,4)
real, dimension(itlo:ithi,5) module_microphysics::qmax

Definition at line 89 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

89  REAL qmax(itlo:ithi,5), qtot(itlo:ithi,22)
real, dimension(itlo:ithi,22) module_microphysics::qtot

Definition at line 89 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

real, parameter, private module_microphysics::rerainmin =1.5*XMRmin
private

Definition at line 61 of file module_bfmicrophysics.f.

Referenced by rsipath(), and rsipath2().

61  REAL, PRIVATE, PARAMETER :: rerainmin=1.5*xmrmin &
real, parameter, private module_microphysics::thom =T_ICE
private

Definition at line 99 of file module_bfmicrophysics.f.

Referenced by rsipath().

99  REAL, PRIVATE, PARAMETER :: thom=t_ice, tnw=50., toler=1.0e-20 &
real, parameter, private module_microphysics::tnw =50.
private

Definition at line 99 of file module_bfmicrophysics.f.

Referenced by rsipath(), and rsipath2().

real, parameter, private module_microphysics::toler =1.0E-20
private

Definition at line 99 of file module_bfmicrophysics.f.

Referenced by gsmcolumn().

real, dimension(2:9,0:nrime), private module_microphysics::vel_rf
private

Definition at line 83 of file module_bfmicrophysics.f.

Referenced by gsmcolumn(), gsmconst(), and ice_lookup().

83  REAL, DIMENSION(2:9,0:Nrime),PRIVATE :: vel_rf