PROGRAM make_p3_lookuptable1 !______________________________________________________________________________________ ! ! This program creates the lookup tables for ice microphysical processes used by ! the P3 microphysics scheme. ! ! Note: compile with double-precision (pgf90 -r8 create_p3_lookupTable_1.f90) ! ! P3 lookup table version: v4.1 ! Last modified : 2019-03-13 !______________________________________________________________________________________ implicit none !-- integer :: i_rhor ! index for rho_rime loop [1 .. n_rhor] integer :: i_Fr ! index for rime-mass-fraction loop [1 .. n_Fr] integer :: i_Qnorm ! index for normalized qitot loop [1 .. n_Qnorm] integer :: i_Drscale ! index for scaled mean rain size loop [1 .. n_Drscale] integer, parameter :: n_rhor = 5 ! number of indices for i_rhor loop (1st; outer) integer, parameter :: n_Fr = 4 ! number of indices for i_Fr loop (2nd) integer, parameter :: n_Qnorm = 50 ! number of indices for i_Qnorm loop (3rd) integer, parameter :: n_Drscale = 30 ! number of indices for scaled mean rain size (4th; inner) !== integer :: i,ii,iii,jj,kk,kkk,dumii real :: N,q,qdum,dum1,dum2,cs1,ds1,lam,n0,lamf,qerror,del0,c0,c1,c2,dd,sum1,sum2, & sum3,sum4,xx,a0,b0,a1,b1,dum,bas1,aas1,aas2,bas2,gammq,gamma,d1,d2,delu,lamold, & cap,lamr,dia,amg,dv,n0dum,sum5,sum6,sum7,sum8,dg,cg,bag,aag,dcritg,dcrits, & dcritr,Fr,csr,dsr,duml,dum3,dum4,rhodep,cgpold,m1,m2,m3,dt,mu_r,initlamr,lamv, & rdumii,lammin,lammax,pi,g,p,t,rho,mu,mu_i,ds,cs,bas,aas,dcrit !-- outputs from lookup table (i.e. "f1prxx" read by access_lookup_table in s/r p3_main) real, dimension(n_Qnorm,n_Fr) :: uns,ums,refl,dmm,rhomm,nagg,nrwat,qsave,nsave,vdep, & eff,lsave,a_100,n_100,vdep1,nsmall,nlarge,nrwats real, dimension(n_Qnorm,n_Drscale,n_Fr) :: qrrain,nrrain,nsrain,qsrain,ngrain !== real, dimension(n_Drscale) :: lamrs real, dimension(10000) :: fall1 real, dimension(1000) :: fall2,fallr,num,numi real, dimension(n_rhor) :: cgp,crp real, dimension(150) :: mu_r_table real, parameter :: Dm_max = 2000.e-6 ! max. mean ice size for lambda limiter real, parameter :: Dm_min = 2.e-6 ! min. mean ice size for lambda limiter ! modify for constant mu_r, set value here real, parameter :: mu_r_constant = 0. real, parameter :: thrd = 1./3. real, parameter :: sxth = 1./6. ! end of variable declaration !------------------------------------------------------------------------------- ! set constants and parameters ! assume 600 hPa, 253 K for p and T for fallspeed calcs (for reference air density) pi = 3.14159 !=acos(-1.) !pi = 3.14159265359 !=acos(-1.) <-- to be replaced with this one (for double precision) g = 9.861 ! gravity p = 60000. ! air pressure (pa) t = 253.15 ! temp (K) rho = p/(287.15*t) ! air density (kg m-3) mu = 1.496E-6*t**1.5/(t+120.)/rho ! viscosity of air dv = 8.794E-5*t**1.81/p ! diffusivity of water vapor in air dt = 10. ! time step for collection (s) ! parameters for surface roughness of ice particle used for fallspeed ! see mitchell and heymsfield 2005 del0 = 5.83 c0 = 0.6 c1 = 4./(del0**2*c0**0.5) c2 = del0**2/4. !--- specified mass-dimension relationship (cgs units) for unrimed crystals: ! ms = cs*D^ds ! ! for graupel: ! mg = cg*D^dg no longer used, since bulk volume is predicted !=== !---- Choice of m-D parameters for large unrimed ice: ! Heymsfield et al. 2006 ! ds=1.75 ! cs=0.0040157+6.06e-5*(-20.) ! sector-like branches (P1b) ! ds=2.02 ! cs=0.00142 ! bullet-rosette ! ds=2.45 ! cs=0.00739 ! side planes ! ds=2.3 ! cs=0.00419 ! radiating assemblages of plates (mitchell et al. 1990) ! ds=2.1 ! cs=0.00239 ! aggreagtes of side planes, bullets, etc. (Mitchell 1996) ! ds=2.1 ! cs=0.0028 !-- ** note: if using one of the above (.i.e. not brown and francis, which is already in mks units), ! then uncomment line below to convert from cgs units to mks ! cs=cs*100.**ds/1000. !== ! Brown and Francis (1995) ds = 1.9 ! cs = 0.01855 ! original (pre v2.3), based on assumption of Dmax cs = 0.0121 ! scaled value based on assumtion of Dmean from Hogan et al. 2012, JAMC !==== ! specify m-D parameter for fully rimed ice ! note: cg is not constant, due to variable density dg = 3. !--- projected area-diam relationship (mks units) for unrimed crystals: ! note: projected area = aas*D^bas ! sector-like branches (P1b) ! bas = 1.97 ! aas = 0.55*100.**bas/(100.**2) ! bullet-rosettes ! bas = 1.57 ! aas = 0.0869*100.**bas/(100.**2) ! aggreagtes of side planes, bullets, etc. bas = 1.88 aas = 0.2285*100.**bas/(100.**2) !=== !--- projected area-diam relationship (mks units) for fully rimed ice: ! note: projected area = aag*D^bag ! assumed non-spherical ! bag = 2.0 ! aag = 0.625*100.**bag/(100.**2) ! assumed spherical: bag = 2. aag = pi*0.25 !=== dcrit = (pi/(6.*cs)*900.)**(1./(ds-3.)) ! check to make sure projected area at dcrit not greater than than of solid sphere ! stop and give warning message if that is the case if (pi/4.*dcrit**2.lt.aas*dcrit**bas) then print*,'STOP, area > area of solid ice sphere, unrimed' stop endif open(unit=1,file='./p3_lookup_table_1.dat-v4.1',status='unknown') !......................................................... ! generate lookup table for mu (for rain) ! ! space of a scaled q/N -- initlamr do i = 1,150 ! loop over lookup table values initlamr = (real(i)*2.)*1.e-6 + 250.e-6 initlamr = 1./initlamr ! iterate to get mu_r ! mu_r-lambda relationship is from Cao et al. (2008), eq. (7) ! start with first guess, mu_r = 0 mu_r = 0. do ii = 1,50 lamr = initlamr*(gamma(mu_r+4.)/(6.*gamma(mu_r+1.)))**thrd ! new estimate for mu_r based on lambda ! set max lambda in formula for mu to 20 mm-1, so Cao et al. ! formula is not extrapolated beyond Cao et al. data range dum = min(20.,lamr*1.e-3) mu_r = max(0.,-0.0201*dum**2+0.902*dum-1.718) ! if lambda is converged within 0.1%, then exit loop if (ii.ge.2) then if (abs((lamold-lamr)/lamr).lt.0.001) goto 111 endif lamold = lamr enddo !ii-loop 111 continue ! assign lookup table values ! modify for constant mu_r = 1 ! mu_r_table(i) = mu_r mu_r_table(i) = mu_r_constant enddo !i-loop !......................................................... ! main loop over rime density ! alpha parameter of m-D for rimed ice crp(1) = 50.*pi*sxth crp(2) = 250.*pi*sxth crp(3) = 450.*pi*sxth crp(4) = 650.*pi*sxth crp(5) = 900.*pi*sxth do i_rhor = 1,n_rhor print*,(pi/(4.*aas))**(1./(bas-2.)) print*,(aas/aag)**(1./(bag-bas)) !------------------------------------------------------------------------ ! find threshold with rimed mass added ! loop over rimed mass fraction (4 points) ! Fr below are values of rime mass fraction for the lookup table ! specific values in model are interpolated between these four points !- note: add this code eventually (outside of i_Fr loop) ! Fr(1) = 0. ! Fr(2) = 0.333 ! Fr(3) = 0.667 ! Fr(4) = 1. != do i_Fr = 1,n_Fr ! loop for rime mass fraction, Fr !-- these lines to be replaced by Fr(i_Fr) initialization outside of loops if (i_Fr.eq.1) Fr = 0. if (i_Fr.eq.2) Fr = 0.333 if (i_Fr.eq.3) Fr = 0.667 if (i_Fr.eq.4) Fr = 1. !== ! calculate mass-dimension relationship for partially-rimed crystals ! msr = csr*D^dsr ! formula from P3 Part 1 (JAS) ! dcritr is critical size separating fully-rimed from partially-rime ice ! first guess, set cgp=crp cgp(i_rhor) = crp(i_rhor) ! case of no riming (Fr = 0), then we need to set dcrits and dcritr to arbitrary large values if (i_Fr.eq.1) then dcrits = 1.e+6 dcritr = dcrits csr = cs dsr = ds ! case of partial riming (Fr between 0 and 1) elseif (i_Fr.eq.2.or.i_Fr.eq.3) then do dcrits = (cs/cgp(i_rhor))**(1./(dg-ds)) dcritr = ((1.+Fr/(1.-Fr))*cs/cgp(i_rhor))**(1./(dg-ds)) csr = cs*(1.+Fr/(1.-Fr)) dsr = ds ! get mean density of vapor deposition/aggregation grown ice rhodep = 1./(dcritr-dcrits)*6.*cs/(pi*(ds-2.))*(dcritr**(ds-2.)-dcrits**(ds-2.)) ! get density of fully-rimed ice as rime mass fraction weighted rime density plus ! density of vapor deposition/aggregation grown ice cgpold = cgp(i_rhor) cgp(i_rhor) = crp(i_rhor)*Fr+rhodep*(1.-Fr)*pi*sxth if (abs((cgp(i_rhor)-cgpold)/cgp(i_rhor)).lt.0.01) goto 115 enddo 115 continue ! case of complete riming (Fr=1.0) else ! set threshold size between partially-rimed and fully-rimed ice as arbitrary large dcrits = (cs/cgp(i_rhor))**(1./(dg-ds)) dcritr = 1.e+6 ! here is the "arbitrary large" csr = cgp(i_rhor) dsr = dg endif print*,'dcrit,dcrits',i_rhor,dcrit,dcrits !--------------------------------------------------------------------------------------- ! set up particle fallspeed arrays ! fallspeed is a function of mass dimension and projected area dimension relationships ! following mitchell and heymsfield (2005), jas ! set up array of particle fallspeed to make computationally efficient ! for high-resolution (in diameter space), ice fallspeed is stored in 'fall1' array (m/s) ! for lower-resolution (in diameter space), ice fallspeed is stored in 'fall2' array (m/s) ! rain fallspeed is stored in 'fallr' (m/s) dd = 2.e-6 ! m ! loop over particle size do jj = 1,10000 ! particle size (m) ! d1 = real(jj)*dd - 0.5*dd d1 = real(jj)*2.*1.e-6 - 1.e-6 ! kept like this to preserve bit-matching !----- get mass-size and projected area-size relationships for given size (d1) ! call get_mass_size if (d1.le.dcrit) then cs1 = pi*sxth*900. ds1 = 3. bas1 = 2. aas1 = pi/4. else if (d1.gt.dcrit.and.d1.le.dcrits) then cs1 = cs ds1 = ds bas1 = bas aas1 = aas else if (d1.gt.dcrits.and.d1.le.dcritr) then cs1 = cgp(i_rhor) ds1 = dg bas1 = bag aas1 = aag else if (d1.gt.dcritr) then cs1 = csr ds1 = dsr if (i_Fr.eq.1) then aas1 = aas bas1 = bas else ! for projected area, keep bas1 constant, but modify aas1 according to rimed fraction bas1 = bas dum1 = aas*d1**bas dum2 = aag*d1**bag m1 = cs1*d1**ds1 m2 = cs*d1**ds m3 = cgp(i_rhor)*d1**dg ! linearly interpolate based on particle mass dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) ! dum3 = (1.-Fr)*dum1+Fr*dum2 aas1 = dum3/(d1**bas) endif endif !===== !----------------------------------------------------------- ! neglect turbulent correction for aggregates for now ! correction for turbulence ! if (d1.lt.500.e-6) then a0 = 0. b0 = 0. ! else ! a0=1.7e-3 ! b0=0.8 ! end if ! fall speed for particle ! best number xx = 2.*cs1*g*rho*d1**(ds1+2.-bas1)/(aas1*(mu*rho)**2) ! drag terms b1 = c1*xx**0.5/(2.*((1.+c1*xx**0.5)**0.5-1.)*(1.+c1*xx**0.5)**0.5)-a0*b0*xx** & b0/(c2*((1.+c1*xx**0.5)**0.5-1.)**2) a1 = (c2*((1.+c1*xx**0.5)**0.5-1.)**2-a0*xx**b0)/xx**b1 ! velocity in terms of drag terms fall1(jj) = a1*mu**(1.-2.*b1)*(2.*cs1*g/(rho*aas1))**b1*d1**(b1*(ds1-bas1+2.)-1.) enddo !jj-loop !......................................................... ! fallspeed array for ice-ice and ice-rain collision calculations do jj = 1,1000 ! particle size d1 = real(jj)*20.*1.e-6 - 10.e-6 if (d1.le.dcrit) then cs1 = pi*sxth*900. ds1 = 3. bas1 = 2. aas1 = pi/4. else if (d1.gt.dcrit.and.d1.le.dcrits) then cs1 = cs ds1 = ds bas1 = bas aas1 = aas else if (d1.gt.dcrits.and.d1.le.dcritr) then cs1 = cgp(i_rhor) ds1 = dg bas1 = bag aas1 = aag else if (d1.gt.dcritr) then cs1 = csr ds1 = dsr if (i_Fr.eq.1) then aas1 = aas bas1 = bas else ! for area, ! keep bas1 constant, but modify aas1 according ! to rimed fraction bas1 = bas dum1 = aas*d1**bas dum2 = aag*d1**bag ! dum3 = (1.-Fr)*dum1+Fr*dum2 m1 = cs1*d1**ds1 m2 = cs*d1**ds m3 = cgp(i_rhor)*d1**dg ! linearly interpolate based on particle mass dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) aas1 = dum3/(d1**bas) endif endif ! correction for turbulence ! if (d1.lt.500.e-6) then a0 = 0. b0 = 0. ! else ! a0=1.7e-3 ! b0=0.8 ! end if ! fall speed for ice ! Best number xx = 2.*cs1*g*rho*d1**(ds1+2.-bas1)/(aas1*(mu*rho)**2) ! drag terms b1 = c1*xx**0.5/(2.*((1.+c1*xx**0.5)**0.5-1.)*(1.+c1*xx**0.5)**0.5)-a0*b0*xx** & b0/(c2*((1.+c1*xx**0.5)**0.5-1.)**2) a1 = (c2*((1.+c1*xx**0.5)**0.5-1.)**2-a0*xx**b0)/xx**b1 ! velocity in terms of drag terms fall2(jj) = a1*mu**(1.-2.*b1)*(2.*cs1*g/(rho*aas1))**b1*d1**(b1*(ds1-bas1+2.)-1.) !--------------------------------------------------------------- ! fall speed for rain particle ! dia = d1*1.e+6 ! HM fix, bug noted by Melissa dia = d1 ! diameter m amg = pi*sxth*997.*dia**3 ! mass kg ! HM fix, bug noted by Melissa amg = amg*1000. ! convert kg to g if(dia.le.134.43e-6) then dum2 = 4.5795e5*amg**(2.*thrd) goto 101 endif if(dia.lt.1511.64e-6) then dum2 = 4.962e3*amg**thrd goto 101 endif if(dia.lt.3477.84e-6) then dum2 = 1.732e3*amg**sxth goto 101 endif dum2 = 917. 101 continue fallr(jj) = dum2*1.e-2 ! convert (cm s-1) to (m s-1) enddo !jj-loop !--------------------------------------------------------------------------------- ! loop around normalized Q (Qnorm) do i_Qnorm = 1,n_Qnorm ! lookup table values of normalized Qnorm ! (range of mean mass diameter from ~ 1 micron to 1 cm) q = 261.7**((i_Qnorm+10)*0.1)*1.e-18 ! uncomment below to test and print proposed values of qovn ! print*,i_Qnorm,(6./(pi*500.)*q)**0.3333 ! end do ! stop ! test values ! N=5.e3 ! q=0.01e-3 print*,'&&&&&&&&&&&i_rhor',i_rhor print*,'***************',i_Qnorm print*,'Fr',Fr print*,'q,N',q,N ! initialize qerror to arbitrarily large value qerror = 1.e+20 !..................................................................................... ! find parameters for gamma distribution ! size distribution for ice is assumed to be ! N(D) = n0 * D^mu_i * exp(-lam*D) ! for the given q and N, we need to find n0, mu_i, and lam ! approach for finding lambda: ! cycle through a range of lambda, find closest lambda to produce correct q ! start with lam, range of lam from 100 to 1 x 10^7 is large enough to ! cover full range over mean size from approximately 1 micron to 1 cm do ii = 1,9000 lam = 1.0013**ii*100. ! get 'mu' parameter (Heymsfield 2003) ! division by 100 is to convert m-1 to cm-1 mu_i = 0.076*(lam/100.)**0.8-2. ! make sure mu_i >= 0, otherwise size dist is infinity at D = 0 mu_i = max(mu_i,0.) ! set upper limit at 6 mu_i = min(mu_i,6.) !-- for lambda limiter: ! set min lam corresponding to 2000 micron for mean size ! dum = Dm_max+Fr*(3000.e-6) dum = Dm_max lam = max(lam,(mu_i+1.)/dum) ! set max lam corresponding to 2 micron mean size lam = min(lam,(mu_i+1.)/Dm_min) ! this range corresponds to range of lam of 500 to 5000000 !== ! get normalized n0 n0 = lam**(mu_i+1.)/(gamma(mu_i+1.)) ! calculate integral for each of the 4 parts of the size distribution ! check difference with respect to Qnorm ! dum1 is integral from 0 to dcrit (solid ice) ! dum2 is integral from dcrit to dcrits (unrimed large ice) ! dum3 is integral from dcrits to dcritr (fully rimed ice) ! dum4 is integral from dcritr to inf (partially rimed) ! set up m-D relationship for solid ice with D < Dcrit cs1 = pi*sxth*900. ds1 = 3. dum1 = lam**(-ds1-mu_i-1.)*gamma(mu_i+ds1+1.)*(1.-gammq(mu_i+ds1+1.,dcrit*lam)) dum2 = lam**(-ds-mu_i-1.)*gamma(mu_i+ds+1.)*(gammq(mu_i+ds+1.,dcrit*lam)) dum = lam**(-ds-mu_i-1.)*gamma(mu_i+ds+1.)*(gammq(mu_i+ds+1.,dcrits*lam)) dum2 = dum2-dum dum2 = max(dum2,0.) dum3 = lam**(-dg-mu_i-1.)*gamma(mu_i+dg+1.)*(gammq(mu_i+dg+1.,dcrits*lam)) dum = lam**(-dg-mu_i-1.)*gamma(mu_i+dg+1.)*(gammq(mu_i+dg+1.,dcritr*lam)) dum3 = dum3-dum dum3 = max(dum3,0.) dum4 = lam**(-dsr-mu_i-1.)*gamma(mu_i+dsr+1.)*(gammq(mu_i+dsr+1.,dcritr*lam)) ! sum of the integrals from the 4 regions of the size distribution qdum = n0*(cs1*dum1+cs*dum2+cgp(i_rhor)*dum3+csr*dum4) ! numerical integration for test to make sure incomplete gamma function is working ! sum1 = 0. ! dd = 1.e-6 ! do iii=1,50000 ! dum=real(iii)*dd ! if (dum.lt.dcrit) then ! sum1 = sum1+n0*dum**mu_i*cs1*dum**ds1* ! 1 exp(-lam*dum)*dd ! else ! sum1 = sum1+n0*dum**mu_i*cs*dum**ds* ! 1 exp(-lam*dum)*dd ! end if ! end do ! print*,'sum1=',sum1 ! stop if (ii.eq.1) then qerror = abs(q-qdum) lamf = lam endif ! find lam with smallest difference between Qnorm and estimate of Qnorm, assign to lamf if (abs(q-qdum).lt.qerror) then lamf = lam qerror = abs(q-qdum) endif enddo !ii-loop ! check and print relative error in q to make sure it is not too large ! note: large error is possible if size bounds are exceeded!!!!!!!!!! print*,'qerror (%)',qerror/q*100. ! find n0 based on final lam value ! set final lamf to 'lam' variable ! this is the value of lam with the smallest qerror lam = lamf ! recalculate mu_i based on final lam ! get 'mu' parameter (Heymsfield 2003) ! division by 100 is to convert m-1 to cm-1 mu_i = 0.076*(lam/100.)**0.8-2. ! make sure mu_i >= 0, otherwise size dist is infinity at D = 0 mu_i = max(mu_i,0.) ! set upper limit at 6 mu_i = min(mu_i,6.) ! n0 = N*lam**(mu_i+1.)/(gamma(mu_i+1.)) ! find n0 from lam and Qnorm ! this is done instead of finding n0 from lam and N, since N ! may need to be adjusted to constrain mean size within reasonable bounds dum1 = lam**(-ds1-mu_i-1.)*gamma(mu_i+ds1+1.)*(1.-gammq(mu_i+ds1+1.,dcrit*lam)) dum2 = lam**(-ds-mu_i-1.)*gamma(mu_i+ds+1.)*(gammq(mu_i+ds+1.,dcrit*lam)) dum = lam**(-ds-mu_i-1.)*gamma(mu_i+ds+1.)*(gammq(mu_i+ds+1.,dcrits*lam)) dum2 = dum2-dum dum3 = lam**(-dg-mu_i-1.)*gamma(mu_i+dg+1.)*(gammq(mu_i+dg+1.,dcrits*lam)) dum = lam**(-dg-mu_i-1.)*gamma(mu_i+dg+1.)*(gammq(mu_i+dg+1.,dcritr*lam)) dum3 = dum3-dum dum4 = lam**(-dsr-mu_i-1.)*gamma(mu_i+dsr+1.)*(gammq(mu_i+dsr+1.,dcritr*lam)) n0 = q/(cs1*dum1+cs*dum2+cgp(i_rhor)*dum3+csr*dum4) print*,'lam,N0:',lam,n0 print*,'mu_i:',mu_i ! test final lam, N0 values ! sum1 = 0. ! dd = 1.e-6 ! do iii=1,50000 ! dum=real(iii)*dd ! if (dum.lt.dcrit) then ! sum1 = sum1+n0*dum**mu_i*cs1*dum**ds1*exp(-lam*dum)*dd ! elseif (dum.ge.dcrit.and.dum.lt.dcrits) then ! sum1 = sum1+n0*dum**mu_i*cs*dum**ds*exp(-lam*dum)*dd ! elseif (dum.ge.dcrits.and.dum.lt.dcritr) then ! sum1 = sum1+n0*dum**mu_i*cg*dum**dg*exp(-lam*dum)*dd ! elseif (dum.ge.dcritr) then ! sum1 = sum1+n0*dum**mu_i*csr*dum**dsr*exp(-lam*dum)*dd ! endif ! enddo ! print*,'sum1=',sum1 ! stop ! At this point, we have solved for all of the size distribution parameters !-------------------------------------------------------------------- ! find max/min N to constrain mean size (i.e. lambda limiter), this is stored ! and passed to lookup table, so that N can be adjusted during ! the simulation to constrain mean size ! limit min size to 2 micron (diameter) duml = (mu_i+1.)/Dm_min dum1 = duml**(-ds1-mu_i-1.)*gamma(mu_i+ds1+1.)*(1.-gammq(mu_i+ds1+1.,dcrit*duml)) dum2 = duml**(-ds-mu_i-1.)*gamma(mu_i+ds+1.)*(gammq(mu_i+ds+1.,dcrit*duml)) dum = duml**(-ds-mu_i-1.)*gamma(mu_i+ds+1.)*(gammq(mu_i+ds+1.,dcrits*duml)) dum2 = dum2-dum dum3 = duml**(-dg-mu_i-1.)*gamma(mu_i+dg+1.)*(gammq(mu_i+dg+1.,dcrits*duml)) dum = duml**(-dg-mu_i-1.)*gamma(mu_i+dg+1.)*(gammq(mu_i+dg+1.,dcritr*duml)) dum3 = dum3-dum dum4 = duml**(-dsr-mu_i-1.)*gamma(mu_i+dsr+1.)*(gammq(mu_i+dsr+1.,dcritr*duml)) n0dum = q/(cs1*dum1+cs*dum2+cgp(i_rhor)*dum3+csr*dum4) ! nlarge is used by P3_MAIN to apply the lambda limiter (upper limit) nlarge(i_Qnorm,i_Fr) = n0dum/(duml**(mu_i+1.)/(gamma(mu_i+1.))) ! limit max size (number-weighted) to 2 mm ! dum = Dm_max+Fr*(3000.e-6) dum = Dm_max duml = (mu_i+1.)/(dum) dum1 = duml**(-ds1-mu_i-1.)*gamma(mu_i+ds1+1.)*(1.-gammq(mu_i+ds1+1.,dcrit*duml)) dum2 = duml**(-ds-mu_i-1.)*gamma(mu_i+ds+1.)*(gammq(mu_i+ds+1.,dcrit*duml)) dum = duml**(-ds-mu_i-1.)*gamma(mu_i+ds+1.)*(gammq(mu_i+ds+1.,dcrits*duml)) dum2 = dum2-dum dum3 = duml**(-dg-mu_i-1.)*gamma(mu_i+dg+1.)*(gammq(mu_i+dg+1.,dcrits*duml)) dum = duml**(-dg-mu_i-1.)*gamma(mu_i+dg+1.)*(gammq(mu_i+dg+1.,dcritr*duml)) dum3 = dum3-dum dum4 = duml**(-dsr-mu_i-1.)*gamma(mu_i+dsr+1.)*(gammq(mu_i+dsr+1.,dcritr*duml)) n0dum = q/(cs1*dum1+cs*dum2+cgp(i_rhor)*dum3+csr*dum4) ! nsmall is used by P3_MAIN to apply the lambda limiter (lower limit) nsmall(i_Qnorm,i_Fr) = n0dum/(duml**(mu_i+1.)/(gamma(mu_i+1.))) print*,'nmax,nmin',nsmall(i_Qnorm,i_Fr),nlarge(i_Qnorm,i_Fr) !..................................................................................... ! begin moment and microphysical process calculations for the lookup table !..................................................................................... ! mass- and number-weighted mean fallspeed (m/s) ! add reflectivity !..................................................................................... ! assume conditions for t and p as assumed above (giving rhos), then in microphysics scheme ! multiply by density correction factor (rhos/rho)^0.54, from Heymsfield et al. 2006 ! fallspeed formulation from Mitchell and Heymsfield 2005 dd = 2.e-6 ! dD for numerical integration, units of m ! initialize for numerical integration sum1 = 0. sum2 = 0. sum3 = 0. sum4 = 0. ! reflectivity sum5 = 0. ! mass mean size sum6 = 0. ! mass-weighted mean density sum7 = 0. ! numerically integrate over size distribution do ii = 1,10000 ! dum = particle size ! dum = real(ii)*dd - 0.5*dd dum = real(ii)*dd - 1.e-6 ! units of m !assign mass-size parameters (depending on size at ii) if (dum.le.dcrit) then ds1 = 3. cs1 = pi*sxth*900. else if (dum.gt.dcrit.and.dum.le.dcrits) then ds1 = ds cs1 = cs elseif (dum.gt.dcrits.and.dum.le.dcritr) then ds1 = dg cs1 = cgp(i_rhor) elseif (dum.gt.dcritr) then ds1 = dsr cs1 = csr endif ! numerator of number-weighted velocity - sum1 sum1 = sum1+fall1(ii)*dum**mu_i*exp(-lam*dum)*dd ! numerator of mass-weighted velocity - sum2 sum2 = sum2+fall1(ii)*cs1*dum**(ds1+mu_i)*exp(-lam*dum)*dd ! total number and mass for weighting above fallspeeds ! don't need to include n0 and cs since these parameters are ! in both numerator and denominator !denominator of number-weighted V sum3 = sum3+dum**mu_i*exp(-lam*dum)*dd !denominator of mass-weighted V sum4 = sum4+cs1*dum**(ds1+mu_i)*exp(-lam*dum)*dd ! reflectivity (integral of mass moment squared) sum5 = sum5+n0*(cs1/917.)**2*(6./pi)**2*dum**(2.*ds1+mu_i)*exp(-lam*dum)*dd ! numerator of mass-weighted mean size sum6 = sum6+cs1*dum**(ds1+mu_i+1.)*exp(-lam*dum)*dd ! numerator of mass-weighted density ! particle density is defined as mass divided by volume of sphere with same D sum7 = sum7+(cs1*dum**ds1)**2/(pi*sxth*dum**3)*dum**mu_i*exp(-lam*dum)*dd enddo !ii-loop ! save mean fallspeeds for lookup table uns(i_Qnorm,i_Fr) = sum1/sum3 ums(i_Qnorm,i_Fr) = sum2/sum4 refl(i_Qnorm,i_Fr) = sum5 dmm(i_Qnorm,i_Fr) = sum6/sum4 rhomm(i_Qnorm,i_Fr) = sum7/sum4 !..................................................................................... ! self-aggregation !..................................................................................... sum1 = 0. dd = 20.e-6 do jj=1000,1,-1 ! set up binned distribution of ice ! d1 = real(jj)*dd - 0.5*dd d1 = real(jj)*dd - 10.e-6 !kept like this to preserve bit-matching num(jj) = n0*d1**mu_i*exp(-lam*d1)*dd enddo !jj-loop ! loop over exponential size distribution ! note: collection of ice within the same bin is neglected do jj = 1000,1,-1 do kk = 1,jj-1 ! particle size d1 = real(jj)*20.*1.e-6 - 10.e-6 d2 = real(kk)*20.*1.e-6 - 10.e-6 if (d1.le.dcrit) then bas1 = 2. aas1 = pi*0.25 elseif (d1.gt.dcrit.and.d1.le.dcrits) then bas1 = bas aas1 = aas else if (d1.gt.dcrits.and.d1.le.dcritr) then bas1 = bag aas1 = aag else if (d1.gt.dcritr) then cs1 = csr ds1 = dsr if (i_Fr.eq.1) then aas1 = aas bas1 = bas else ! for area, keep bas1 constant, but modify aas1 according to rimed fraction bas1 = bas dum1 = aas*d1**bas dum2 = aag*d1**bag m1 = cs1*d1**ds1 m2 = cs*d1**ds m3 = cgp(i_rhor)*d1**dg ! linearly interpolate based on particle mass dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) aas1 = dum3/(d1**bas) endif endif ! parameters for particle 2 if (d2.le.dcrit) then bas2 = 2. aas2 = pi/4. elseif (d2.gt.dcrit.and.d2.le.dcrits) then bas2 = bas aas2 = aas elseif (d2.gt.dcrits.and.d2.le.dcritr) then bas2 = bag aas2 = aag elseif (d2.gt.dcritr) then cs1 = csr ds1 = dsr if (i_Fr.eq.1) then aas1 = aas bas1 = bas else ! for area, ! keep bas1 constant, but modify aas1 according ! to rimed fraction bas2 = bas dum1 = aas*d2**bas dum2 = aag*d2**bag m1 = cs1*d1**ds1 m2 = cs*d1**ds m3 = cgp(i_rhor)*d1**dg ! linearly interpolate based on particle mass dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) aas2 = dum3/(d2**bas) endif endif ! absolute value, differential fallspeed delu = abs(fall2(jj)-fall2(kk)) ! note: in micro code we have to multiply by air density ! correction factor for fallspeed, and collection efficiency ! sum for integral ! sum1 = # of collision pairs ! the assumption is that each collision pair reduces crystal ! number mixing ratio by 1 kg^-1 s^-1 per kg/m^3 of air (this is ! why we need to multiply by air density, to get units of ! 1/kg^-1 s^-1) sum1 = sum1+(aas1*d1**bas1+aas2*d2**bas2)*delu*num(jj)*num(kk) ! sum1 = sum1+(aas1*d1**bas1+aas2*d2**bas2)*delu*num(jj)*num(kk) ! remove collected particles from distribution over time period dt, update num ! note -- dt is time scale for removal, not model time step ! num(kk) = num(kk)-(aas1*d1**bas1+aas2*d2**bas2)*delu*num(jj)*num(kk)*dt ! num(kk) = max(num(kk),0.) ! write(6,'(2i5,8e15.5)')jj,kk,sum1,num(jj),num(kk),delu,aas1,d1,aas2,d2 ! num(kk)=num(kk)-(aas1*d1**bas1+aas2*d2**bas2)*delu*num(jj)*num(kk)*0.1*0.5 ! num(kk)=max(num(kk),0.) ! sum1 = sum1+0.5*(aas1*d1**bas1+aas2*d2**bas2)*delu*n0*n0*(d1+d2)**mu_i*exp(-lam*(d1+d2))*dd**2 enddo !kk-loop enddo !jj-loop ! save for output nagg(i_Qnorm,i_Fr) = sum1 print*,'nagg',nagg(i_Qnorm,i_Fr) !..................................................................................... ! collection of cloud droplets !..................................................................................... ! note: in microhpysics code needs to be multiplied by collection efficiency Eci ! note: also needs to be multiplied by air density correction factor for fallspeed, ! air density, and cloud water mixing ratio or number concentration ! initialize sum for integral sum1 = 0. sum2 = 0. ! dd (dD for numerical integration) (m) dd = 2.e-6 ! loop over exponential size distribution (from 1 micron to 2 cm) do jj = 1,10000 ! particle size or dimension (m) for numerical integration d1 = real(jj)*2.*1.e-6 - 1.e-6 ! get mass-dimension and projected area-dimension relationships ! for different ice types across the size distribution based on critical dimensions ! separating these ice types (see Fig. 2, morrison and grabowski 2008) ! mass = cs1*D^ds1 ! projected area = bas1*D^bas1 if (d1.le.dcrit) then cs1 = pi*sxth*900. ds1 = 3. bas1 = 2. aas1 = pi*0.25 elseif (d1.gt.dcrit.and.d1.le.dcrits) then cs1 = cs ds1 = ds bas1 = bas aas1 = aas elseif (d1.gt.dcrits.and.d1.le.dcritr) then cs1 = cgp(i_rhor) ds1 = dg bas1 = bag aas1 = aag elseif (d1.gt.dcritr) then cs1 = csr ds1 = dsr if (i_Fr.eq.1) then aas1 = aas bas1 = bas else ! for area, ! keep bas1 constant, but modify aas1 according to rimed fraction bas1 = bas dum1 = aas*d1**bas dum2 = aag*d1**bag m1 = cs1*d1**ds1 m2 = cs*d1**ds m3 = cgp(i_rhor)*d1**dg ! linearly interpolate based on particle mass dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) aas1 = dum3/(d1**bas) endif endif ! sum for integral ! include assumed ice particle size threshold for riming of 100 micron ! note: sum1 (nrwat) is the scaled collection rate, units of m^3 kg^-1 s^-1 ! sum2 (nrwats) is mass of snow times scaled collection rate, ! units of m^3 s^-1 if (d1.ge.100.e-6) then sum1 = sum1+aas1*d1**bas1*fall1(jj)*n0*d1**mu_i*exp(-lam*d1)*dd ! sum2 = sum2+aas1*d1**bas1*fall1(jj)*n0*d1**mu_i*exp(-lam*d1)*dd*cs1*d1**ds1 endif enddo !jj-loop ! save for output ! note: read in as 'f1pr4' in P3_MAIN nrwat(i_Qnorm,i_Fr) = sum1 ! nrwats(i_Qnorm,i_Fr) = sum2 print*,'nrwat',nrwat(i_Qnorm,i_Fr) !..................................................................................... ! collection of rain !..................................................................................... ! note: in microphysics code, we need to multiply rate by n0r, collection efficiency, ! air density, and air density correction factor ! this approach implicitly assumes that the PSDs are constant during the microphysics ! time step, this could produce errors if the time step is large. In particular, ! more mass or number could be removed than is available. This will be taken care ! of by conservation checks in the microphysics code. ! loop around lambda for rain do i_Drscale = 1,n_Drscale dum = 1.24**i_Drscale*10.e-6 ! assumed lamv for tests ! dum = 7.16e-5 ! Note: the lookup table for rain is based on lamv, i.e., ! inverse volume mean diameter lamv = 1./dum lamrs(i_Drscale) = lamv !------------------------------------------------------------ ! get mu_r from lamr dum = 1./lamv ! if (dum.lt.282.e-6) then ! mu_r = 8.282 ! elseif (dum.ge.282.e-6 .and. dum.lt.502.e-6) then !! interpolate ! rdumii = (dum-250.e-6)*1.e6*0.5 ! rdumii = max(rdumii,1.) ! rdumii = min(rdumii,150.) ! dumii = int(rdumii) ! dumii = min(149,dumii) ! mu_r = mu_r_table(dumii)+(mu_r_table(dumii+1)-mu_r_table(dumii))*(rdumii-real(dumii)) ! elseif (dum.ge.502.e-6) then ! mu_r = 0. ! endif ! modify for constant mu_r mu_r = mu_r_constant ! recalculate slope based on mu_r ! LAMR = (pi*sxth*rhow*nr(i_Qnorm,k)*gamma(mu_r+4.)/(qr(i_Qnorm,k)*gamma(mu_r+1.)))**thrd ! this is done by re-scaling lamv to account for DSD shape (mu_r) lamr = lamv*(gamma(mu_r+4.)/(6.*gamma(mu_r+1.)))**thrd ! set maximum value for rain lambda ! lammax = (mu_r+1.)/10.e-6 lammax = (mu_r+1.)*1.e+5 ! set to small value since breakup is explicitly included (mean size 5 mm) ! lammin = (mu_r+1.)/5000.e-6 lammin = (mu_r+1.)*200. lamr = min(lamr,lammax) lamr = max(lamr,lammin) ! initialize sum sum1 = 0. sum2 = 0. sum6 = 0. ! total rain ! sum8=0. dd = 20.e-6 do jj=1,1000 ! particle size d1 = real(jj)*20.*1.e-6 - 10.e-6 ! num is the scaled binned rain size distribution ! we need to multiply by n0r to get unscaled distribution num(jj) = d1**mu_r*exp(-lamr*d1)*dd ! get (unscaled) binned ice size distribution ! bug fix hm, 2/10/12 numi(jj) = n0*d1**mu_i*exp(-lam*d1)*dd enddo !jj-loop ! loop over rain and ice size distributions do jj = 1,1000 do kk = 1,1000 ! particle size d1 = real(jj)*20.*1.e-6 - 10.e-6 ! ice d2 = real(kk)*20.*1.e-6 - 10.e-6 ! rain if (d1.le.dcrit) then cs1 = pi*sxth*900. ds1 = 3. bas1 = 2. aas1 = pi*0.25 elseif (d1.gt.dcrit.and.d1.le.dcrits) then cs1 = cs ds1 = ds bas1 = bas aas1 = aas elseif (d1.gt.dcrits.and.d1.le.dcritr) then cs1 = cgp(i_rhor) ds1 = dg bas1 = bag aas1 = aag else if (d1.gt.dcritr) then cs1 = csr ds1 = dsr if (i_Fr.eq.1) then aas1 = aas bas1 = bas else ! for area, keep bas1 constant, but modify aas1 accordingto rimed fraction bas1 = bas dum1 = aas*d1**bas dum2 = aag*d1**bag m1 = cs1*d1**ds1 m2 = cs*d1**ds m3 = cgp(i_rhor)*d1**dg ! linearly interpolate based on particle mass dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) aas1 = dum3/(d1**bas) endif endif ! absolute value, differential fallspeed delu = abs(fall2(jj)-fallr(kk)) ! collection of rain mass and number ! allow collection of rain both when rain fallspeed > ice fallspeed ! and ice fallspeed > rain fallspeed ! this is applied below freezing to calculate total amount of rain ! mass and number that collides with ice and freezes ! if (fall2(jj).ge.fallr(kk)) then ! sum for integral ! change in rain N (units of m^4 s^-1 kg^-1), thus need to multiply ! by air density (units kg m^-3) and n0r (units kg^-1 m^-1) in the ! microphysics code ! sum1 = sum1+(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* & ! exp(-lam*d1)* &dd*num(kk) sum1 = sum1+(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* & exp(-lam*d1)*dd*num(kk) ! sum1 = sum1+min((aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* & ! exp(-lam*d1)*dd*num(kk),num(kk)) ! change in rain q (units of m^4 s^-1), again need to multiply ! by air density and n0r in the microphysics code ! sum2 = sum2+(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* & ! exp(-lam*d1)*dd*num(kk)*pi*sxth*997.*d2**3 sum2 = sum2+(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* & exp(-lam*d1)*dd*num(kk)*pi*sxth*997.*d2**3 ! remove collected rain drops from distribution ! num(kk) = num(kk)-(aas1*d1**bas1+pi*0.25*d2**2)*delu*n0*d1**mu_i* & ! exp(-lam*d1)*dd*num(kk)*dt ! num(kk) = max(num(kk),0.) !...................................................... ! now calculate collection of ice mass by rain ! ice collecting rain ! again, allow collection both when ice fallspeed > rain fallspeed ! and when rain fallspeed > ice fallspeed ! this is applied to conditions above freezing to calculate ! acceleration of melting due to collisions with liquid (rain) ! if (fall2(jj).ge.fallr(kk)) then ! collection of ice number ! sum5=sum5+(aas1*d1**bas1+pi/4.*d2**2)*delu* ! 1 exp(-lamr*d2)*dd*numi(jj) ! collection of ice mass (units of m^4 s^-1), need to multiply by air density and n0r in microphysics code sum6 = sum6+(aas1*d1**bas1+pi*0.25*d2**2)*delu*d2**mu_r* & exp(-lamr*d2)*dd*numi(jj)*cs1*d1**ds1 ! remove collected snow from distribution ! numi(jj) = numi(jj)-(aas1*d1**bas1+pi*0.25*d2**2)*delu*d2**mu_r* & ! exp(-lamr*d2)*dd*numi(jj)*dt ! numi(jj) = max(numi(jj),0.) enddo !kk-loop enddo !jj-loop ! save for output nrrain(i_Qnorm,i_Drscale,i_Fr) = sum1 qrrain(i_Qnorm,i_Drscale,i_Fr) = sum2 qsrain(i_Qnorm,i_Drscale,i_Fr) = sum6 enddo !i_Drscale-loop (loop around lambda for rain) ! !..................................................................................... ! vapor deposition/melting/wet growth !..................................................................................... ! vapor deposition including ventilation effects ! note: in microphysics code we need to multiply by air density ! and (mu/dv)^0.3333*(rhofac/mu)^0.5, where rhofac is air density ! correction factor sum1 = 0. sum2 = 0. dd = 2.e-6 ! loop over exponential size distribution do jj = 1,10000 ! particle size d1 = real(jj)*2.*1.e-6 - 1.e-6 ! get capacitance for different ice regimes if (d1.le.dcrit) then cap = 1. ! for small spherical crystal use sphere elseif (d1.gt.dcrit.and.d1.le.dcrits) then cap = 0.48 ! field et al. 2006 elseif (d1.gt.dcrits.and.d1.le.dcritr) then cap = 1. ! for graupel assume sphere elseif (d1.gt.dcritr) then cs1 = csr ds1 = dsr if (i_Fr.eq.1) then cap = 0.48 else dum1 = 0.48 dum2 = 1. m1 = cs1*d1**ds1 m2 = cs*d1**ds m3 = cgp(i_rhor)*d1**dg ! linearly interpolate to get capacitance based on particle mass dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) cap = dum3 endif endif ! hm 1/19/13 for ventilation, only include fallspeed ! and size effects, the rest of the term ! Sc^1/3 x Re^1/2 is multiplied in-line in the model code ! to allow effects of atmospheric conditions on ventilation ! dum = (mu/dv)**0.333333*(fall1(jj)*d1/mu)**0.5 dum = (fall1(jj)*d1)**0.5 ! ventilation from hall and pruppacher 1976 ! hm 1/19/13 only include ventilation for super-100 micron particles ! if (dum.lt.1.) then ! units are m^3 kg^-1 s^-1, thus multiplication by air density in ! microphysics code gives s^-1 if (d1.lt.100.e-6) then sum1 = sum1+cap*n0*d1**(mu_i+1.)*exp(-lam*d1)*dd else ! sum1 = sum1+cap*n0*(0.65+0.44*dum)*d1**(mu_i+1.)*exp(-lam*d1)*dd sum1 = sum1+cap*n0*0.65*d1**(mu_i+1.)*exp(-lam*d1)*dd sum2 = sum2+cap*n0*0.44*dum*d1**(mu_i+1.)*exp(-lam*d1)*dd endif enddo !jj-loop vdep(i_Qnorm,i_Fr) = sum1 vdep1(i_Qnorm,i_Fr) = sum2 print*,'vdep',vdep(i_Qnorm,i_Fr) !..................................................................................... ! ice effective radius ! use definition of Francis et al. (1994), e.g., Eq. 3.11 in Fu (1996) J. Climate !..................................................................................... sum1 = 0. sum2 = 0. dd = 2.e-6 ! loop over exponential size distribution do jj = 1,10000 ! particle size d1 = real(jj)*2.*1.e-6 - 1.e-6 if (d1.le.dcrit) then cs1 = pi*sxth*900. ds1 = 3. bas1 = 2. aas1 = pi*0.25 elseif (d1.gt.dcrit.and.d1.le.dcrits) then cs1 = cs ds1 = ds bas1 = bas aas1 = aas elseif (d1.gt.dcrits.and.d1.le.dcritr) then cs1 = cgp(i_rhor) ds1 = dg bas1 = bag aas1 = aag elseif (d1.gt.dcritr) then cs1 = csr ds1 = dsr if (i_Fr.eq.1) then bas1 = bas aas1 = aas else ! for area, keep bas1 constant, but modify aas1 according to rimed fraction bas1 = bas dum1 = aas*d1**bas dum2 = aag*d1**bag m1 = cs1*d1**ds1 m2 = cs*d1**ds m3 = cgp(i_rhor)*d1**dg ! linearly interpolate based on particle mass dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) aas1 = dum3/(d1**bas) endif endif ! n0 not included below becuase it is in both numerator and denominator ! cs1=pi*sxth*917. ! ds1=3. ! aas1=pi/4.*2. ! bas1=2. sum1 = sum1+cs1*d1**ds1*d1**mu_i*exp(-lam*d1)*dd sum2 = sum2+aas1*d1**bas1*d1**mu_i*exp(-lam*d1)*dd ! if (d1.ge.100.e-6) then ! sum3 = sum3+n0*d1**mu_i*exp(-lam*d1)*dd ! sum4 = sum4+n0*aas1*d1**bas1*d1**mu_i*exp(-lam*d1)*dd ! endif enddo !jj-loop ! calculate eff radius ! eff(i_Qnorm,i_Fr) = sum1/(1.7321*916.7*sum2) ! hm 4/9/09, calculate effective size following Fu (1996) ! eff(i_Qnorm,i_Fr) = sum1/(1.1547*916.7*sum2) ! hm, calculate for eff rad for twp ice eff(i_Qnorm,i_Fr) = 3.*sum1/(4.*sum2*916.7) ! a_100(i_Qnorm,i_Fr)=sum4 ! n_100(i_Qnorm,i_Fr)=sum3 print*,'eff rad',eff(i_Qnorm,i_Fr) !..................................................................................... 522 continue !-- this column is not used (for v2.2 and after) ! (kept temporarly in order to preserve order of columns in lookup_table_1) nsave(i_Qnorm,i_Fr) = 1. ! HM, set to 1 for consistency w/ old lookup table file (to verify bit matching) !== qsave(i_Qnorm,i_Fr) = q lsave(i_Qnorm,i_Fr) = lam enddo !i_Qnorm-loop enddo !i_Fr-loop ! output variables to ascii lookup table 222 format(2i5,15e15.5) 223 format(2i5,6e15.5) !-- ice table do i_Fr=1,n_Fr do i_Qnorm = 1,n_Qnorm ! set values less than 1.e-99 set to 0, otherwise the 'E' is left off in write statements for floting point ! numbers using some compilers if (qsave(i_Qnorm,i_Fr).lt.1.e-99) qsave(i_Qnorm,i_Fr)=0. if (nsave(i_Qnorm,i_Fr).lt.1.e-99) nsave(i_Qnorm,i_Fr)=0. if (uns(i_Qnorm,i_Fr).lt.1.e-99) uns(i_Qnorm,i_Fr)=0. if (ums(i_Qnorm,i_Fr).lt.1.e-99) ums(i_Qnorm,i_Fr)=0. if (nagg(i_Qnorm,i_Fr).lt.1.e-99) nagg(i_Qnorm,i_Fr)=0. if (nrwat(i_Qnorm,i_Fr).lt.1.e-99) nrwat(i_Qnorm,i_Fr)=0. if (vdep(i_Qnorm,i_Fr).lt.1.e-99) vdep(i_Qnorm,i_Fr)=0. if (eff(i_Qnorm,i_Fr).lt.1.e-99) eff(i_Qnorm,i_Fr)=0. if (nlarge(i_Qnorm,i_Fr).lt.1.e-99) nlarge(i_Qnorm,i_Fr)=0. if (nsmall(i_Qnorm,i_Fr).lt.1.e-99) nsmall(i_Qnorm,i_Fr)=0. if (lsave(i_Qnorm,i_Fr).lt.1.e-99) lsave(i_Qnorm,i_Fr)=0. if (refl(i_Qnorm,i_Fr).lt.1.e-99) refl(i_Qnorm,i_Fr)=0. if (vdep1(i_Qnorm,i_Fr).lt.1.e-99) vdep1(i_Qnorm,i_Fr)=0. if (dmm(i_Qnorm,i_Fr).lt.1.e-99) dmm(i_Qnorm,i_Fr)=0. if (rhomm(i_Qnorm,i_Fr).lt.1.e-99) rhomm(i_Qnorm,i_Fr)=0. write(1,222)i_rhor,i_Fr,qsave(i_Qnorm,i_Fr),nsave(i_Qnorm,i_Fr),uns(i_Qnorm,i_Fr), & ums(i_Qnorm,i_Fr),nagg(i_Qnorm,i_Fr),nrwat(i_Qnorm,i_Fr),vdep(i_Qnorm,i_Fr), & eff(i_Qnorm,i_Fr),nlarge(i_Qnorm,i_Fr),nsmall(i_Qnorm,i_Fr),lsave(i_Qnorm,i_Fr), & refl(i_Qnorm,i_Fr),vdep1(i_Qnorm,i_Fr),dmm(i_Qnorm,i_Fr),rhomm(i_Qnorm,i_Fr) enddo !i-loop !-- ice-rain collection table: do i_Qnorm = 1,n_Qnorm do i_Drscale=1,n_Drscale ! set values less than 1.e-99 set to 0, otherwise the 'E' is left off in write statements for floting point ! numbers using some compilers if (qsave(i_Qnorm,i_Fr).lt.1.e-99) qsave(i_Qnorm,i_Fr)=0. if (nsave(i_Qnorm,i_Fr).lt.1.e-99) nsave(i_Qnorm,i_Fr)=0. if (lamrs(i_Drscale).lt.1.e-99) lamrs(i_Drscale)=0. if (nrrain(i_Qnorm,i_Drscale,i_Fr).lt.1.e-99) nrrain(i_Qnorm,i_Drscale,i_Fr)=0. if (qrrain(i_Qnorm,i_Drscale,i_Fr).lt.1.e-99) qrrain(i_Qnorm,i_Drscale,i_Fr)=0. if (qsrain(i_Qnorm,i_Drscale,i_Fr).lt.1.e-99) qsrain(i_Qnorm,i_Drscale,i_Fr)=0. write(1,223) i_rhor,i_Fr,qsave(i_Qnorm,i_Fr),nsave(i_Qnorm,i_Fr),lamrs(i_Drscale), & nrrain(i_Qnorm,i_Drscale,i_Fr),qrrain(i_Qnorm,i_Drscale,i_Fr), & qsrain(i_Qnorm,i_Drscale,i_Fr) enddo !i_Drscale-loop enddo !i-loop enddo !i_Fr-loop (riming fraction) enddo !i_rhor-loop (main loop over variable rime density) END PROGRAM make_p3_lookuptable1 !______________________________________________________________________________________ ! Incomplete gamma function ! from Numerical Recipes in Fortran 77: The Art of ! Scientific Computing function gammq(a,x) real a,gammq,x ! USES gcf,gser ! Returns the incomplete gamma function Q(a,x) = 1-P(a,x) real gammcf,gammser,gln if (x.lt.0..or.a.le.0) pause 'bad argument in gammq' if (x.lt.a+1.) then call gser(gamser,a,x,gln) gammq=1.-gamser else call gcf(gammcf,a,x,gln) gammq=gammcf end if return end !------------------------------------- subroutine gser(gamser,a,x,gln) integer itmax real a,gamser,gln,x,eps parameter(itmax=100,eps=3.e-7) integer n real ap,del,sum,gamma gln = log(gamma(a)) if (x.le.0.) then if (x.lt.0.) pause 'x < 0 in gser' gamser = 0. return end if ap=a sum=1./a del=sum do n=1,itmax ap=ap+1. del=del*x/ap sum=sum+del if (abs(del).lt.abs(sum)*eps) goto 1 end do pause 'a too large, itmax too small in gser' 1 gamser=sum*exp(-x+a*log(x)-gln) return end !------------------------------------- subroutine gcf(gammcf,a,x,gln) integer itmax real a,gammcf,gln,x,eps,fpmin parameter(itmax=100,eps=3.e-7,fpmin=1.e-30) integer i real an,b,c,d,del,h,gamma gln=log(gamma(a)) b=x+1.-a c=1./fpmin d=1./b h=d do i=1,itmax an=-i*(i-a) b=b+2. d=an*d+b if(abs(d).lt.fpmin) d=fpmin c=b+an/c if(abs(c).lt.fpmin) c=fpmin d=1./d del=d*c h = h*del if(abs(del-1.).lt.eps)goto 1 end do pause 'a too large, itmax too small in gcf' 1 gammcf=exp(-x+a*log(x)-gln)*h return end !------------------------------------- REAL FUNCTION gamma(X) !D DOUBLE PRECISION FUNCTION Dgamma(X) !---------------------------------------------------------------------- ! ! THIS ROUTINE CALCULATES THE gamma FUNCTION FOR A REAL ARGUMENT X. ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1. ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE gamma ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2. ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE ! MACHINE-DEPENDENT CONSTANTS. ! ! !---------------------------------------------------------------------- ! ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS ! ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS ! XBIG - THE LARGEST ARGUMENT FOR WHICH gamma(X) IS REPRESENTABLE ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION ! gamma(XBIG) = BETA**MAXEXP ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER; ! APPROXIMATELY BETA**MAXEXP ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1.0+EPS .GT. 1.0 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1/XMININ IS MACHINE REPRESENTABLE ! ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: ! ! BETA MAXEXP XBIG ! ! CRAY-1 (S.P.) 2 8191 966.961 ! CYBER 180/855 ! UNDER NOS (S.P.) 2 1070 177.803 ! IEEE (IBM/XT, ! SUN, ETC.) (S.P.) 2 128 35.040 ! IEEE (IBM/XT, ! SUN, ETC.) (D.P.) 2 1024 171.624 ! IBM 3033 (D.P.) 16 63 57.574 ! VAX D-FORMAT (D.P.) 2 127 34.844 ! VAX G-FORMAT (D.P.) 2 1023 171.489 ! ! XINF EPS XMININ ! ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 ! CYBER 180/855 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 ! IEEE (IBM/XT, ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 ! IEEE (IBM/XT, ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308 ! !---------------------------------------------------------------------- ! ! ERROR RETURNS ! ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED ! TO BE FREE OF UNDERFLOW AND OVERFLOW. ! ! ! INTRINSIC FUNCTIONS REQUIRED ARE: ! ! INT, DBLE, EXP, LOG, REAL, SIN ! ! ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS, ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON ! (ED.), SPRINGER VERLAG, BERLIN, 1976. ! ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND ! SONS, NEW YORK, 1968. ! ! LATEST MODIFICATION: OCTOBER 12, 1989 ! ! AUTHORS: W. J. CODY AND L. STOLTZ ! APPLIED MATHEMATICS DIVISION ! ARGONNE NATIONAL LABORATORY ! ARGONNE, IL 60439 ! !---------------------------------------------------------------------- INTEGER I,N LOGICAL PARITY REAL & !D DOUBLE PRECISION C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, & TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DIMENSION C(7),P(8),Q(8) !---------------------------------------------------------------------- ! MATHEMATICAL CONSTANTS !---------------------------------------------------------------------- DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/, & SQRTPI/0.9189385332046727417803297E0/, & PI/3.1415926535897932384626434E0/ !D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/, !D 1 SQRTPI/0.9189385332046727417803297D0/, !D 2 PI/3.1415926535897932384626434D0/ !---------------------------------------------------------------------- ! MACHINE DEPENDENT PARAMETERS !---------------------------------------------------------------------- DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/, & XINF/3.4E38/ !D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/, !D 1 XINF/1.79D308/ !---------------------------------------------------------------------- ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX ! APPROXIMATION OVER (1,2). !---------------------------------------------------------------------- DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, & -3.79804256470945635097577E+2,6.29331155312818442661052E+2, & 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, & -3.61444134186911729807069E+4,6.64561438202405440627855E+4/ DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, & -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, & 2.25381184209801510330112E+4,4.75584627752788110767815E+3, & -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/ !D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1, !D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2, !D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, !D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ !D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2, !D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3, !D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3, !D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ !---------------------------------------------------------------------- ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). !---------------------------------------------------------------------- DATA C/-1.910444077728E-03,8.4171387781295E-04, & -5.952379913043012E-04,7.93650793500350248E-04, & -2.777777777777681622553E-03,8.333333333333333331554247E-02, & 5.7083835261E-03/ !D DATA C/-1.910444077728D-03,8.4171387781295D-04, !D 1 -5.952379913043012D-04,7.93650793500350248D-04, !D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, !D 3 5.7083835261D-03/ !---------------------------------------------------------------------- ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT !---------------------------------------------------------------------- CONV(I) = REAL(I) !D CONV(I) = DBLE(I) PARITY=.FALSE. FACT=ONE N=0 Y=X IF(Y.LE.ZERO)THEN !---------------------------------------------------------------------- ! ARGUMENT IS NEGATIVE !---------------------------------------------------------------------- Y=-X Y1=AINT(Y) RES=Y-Y1 IF(RES.NE.ZERO)THEN IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE. FACT=-PI/SIN(PI*RES) Y=Y+ONE ELSE RES=XINF GOTO 900 ENDIF ENDIF !---------------------------------------------------------------------- ! ARGUMENT IS POSITIVE !---------------------------------------------------------------------- IF(Y.LT.EPS)THEN !---------------------------------------------------------------------- ! ARGUMENT .LT. EPS !---------------------------------------------------------------------- IF(Y.GE.XMININ)THEN RES=ONE/Y ELSE RES=XINF GOTO 900 ENDIF ELSEIF(Y.LT.TWELVE)THEN Y1=Y IF(Y.LT.ONE)THEN !---------------------------------------------------------------------- ! 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- Z=Y Y=Y+ONE ELSE !---------------------------------------------------------------------- ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY !---------------------------------------------------------------------- N=INT(Y)-1 Y=Y-CONV(N) Z=Y-ONE ENDIF !---------------------------------------------------------------------- ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 !---------------------------------------------------------------------- XNUM=ZERO XDEN=ONE DO 260 I=1,8 XNUM=(XNUM+P(I))*Z XDEN=XDEN*Z+Q(I) 260 CONTINUE RES=XNUM/XDEN+ONE IF(Y1.LT.Y)THEN !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- RES=RES/Y1 ELSEIF(Y1.GT.Y)THEN !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0 !---------------------------------------------------------------------- DO 290 I=1,N RES=RES*Y Y=Y+ONE 290 CONTINUE ENDIF ELSE !---------------------------------------------------------------------- ! EVALUATE FOR ARGUMENT .GE. 12.0, !---------------------------------------------------------------------- IF(Y.LE.XBIG)THEN YSQ=Y*Y SUM=C(7) DO 350 I=1,6 SUM=SUM/YSQ+C(I) 350 CONTINUE SUM=SUM/Y-Y+SQRTPI SUM=SUM+(Y-HALF)*LOG(Y) RES=EXP(SUM) ELSE RES=XINF GOTO 900 ENDIF ENDIF !---------------------------------------------------------------------- ! FINAL ADJUSTMENTS AND RETURN !---------------------------------------------------------------------- IF(PARITY)RES=-RES IF(FACT.NE.ONE)RES=FACT/RES 900 gamma=RES !D900 Dgamma = RES RETURN ! ---------- LAST LINE OF gamma ---------- END !-------------------------------------------------------------------------- ! subroutine get_mass_size ! ! !----- get mass-size and projected area-size relationships for given size (d1) ! if (d1.le.dcrit) then ! cs1 = pi*sxth*900. ! ds1 = 3. ! bas1 = 2. ! aas1 = pi/4. ! else if (d1.gt.dcrit.and.d1.le.dcrits) then ! cs1 = cs ! ds1 = ds ! bas1 = bas ! aas1 = aas ! else if (d1.gt.dcrits.and.d1.le.dcritr) then ! cs1 = cgp(i_rhor) ! ds1 = dg ! bas1 = bag ! aas1 = aag ! else if (d1.gt.dcritr) then ! cs1 = csr ! ds1 = dsr ! if (i_Fr.eq.1) then ! aas1 = aas ! bas1 = bas ! else ! ! ! for projected area, keep bas1 constant, but modify aas1 according to rimed fraction ! bas1 = bas ! dum1 = aas*d1**bas ! dum2 = aag*d1**bag ! m1 = cs1*d1**ds1 ! m2 = cs*d1**ds ! m3 = cgp(i_rhor)*d1**dg ! ! linearly interpolate based on particle mass ! dum3 = dum1+(m1-m2)*(dum2-dum1)/(m3-m2) ! ! dum3 = (1.-Fr)*dum1+Fr*dum2 ! aas1 = dum3/(d1**bas) ! endif ! endif ! !===== ! ! end subroutine get_mass_size