!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------- !------------ FMS version number and tagname for this file ----------- ! $Id: dsd.f90,v 1.1.2.1.2.1 2009/08/10 10:48:14 rsh Exp $ ! $Name: mom4p1_pubrel_dec2009_nnz $ subroutine dsd(Q,Re,D,N,nsizes,dtype,rho_a,tc, & dmin,dmax,apm,bpm,rho_c,p1,p2,p3,fc,scaled) use array_lib use math_lib implicit none ! Purpose: ! Create a discrete drop size distribution ! Part of QuickBeam v1.03 by John Haynes ! http://reef.atmos.colostate.edu/haynes/radarsim ! ! Inputs: ! [Q] hydrometeor mixing ratio (g/kg) ! [Re] Optional Effective Radius (microns). 0 = use default. ! [D] discrete drop sizes (um) ! [nsizes] number of elements of [D] ! [dtype] distribution type ! [rho_a] ambient air density (kg m^-3) ! [tc] temperature (C) ! [dmin] minimum size cutoff (um) ! [dmax] maximum size cutoff (um) ! [rho_c] alternate constant density (kg m^-3) ! [p1],[p2],[p3] distribution parameters ! ! Input/Output: ! [fc] scaling factor for the distribution ! [scaled] has this hydrometeor type been scaled? ! [apm] a parameter for mass (kg m^[-bpm]) ! [bmp] b params for mass ! ! Outputs: ! [N] discrete concentrations (cm^-3 um^-1) ! or, for monodisperse, a constant (1/cm^3) ! ! Requires: ! function infind ! ! Created: ! 11/28/05 John Haynes (haynes@atmos.colostate.edu) ! Modified: ! 01/31/06 Port from IDL to Fortran 90 ! 07/07/06 Rewritten for variable DSD's ! 10/02/06 Rewritten using scaling factors (Roger Marchand and JMH) ! ----- INPUTS ----- integer*4, intent(in) :: nsizes integer, intent(in) :: dtype real*8, intent(in) :: Q,D(nsizes),rho_a,tc,dmin,dmax, & rho_c,p1,p2,p3 ! ----- INPUT/OUTPUT ----- real*8, intent(inout) :: fc(nsizes),apm,bpm,Re logical, intent(inout) :: scaled ! ----- OUTPUTS ----- real*8, intent(out) :: N(nsizes) ! ----- INTERNAL ----- real*8 :: & N0,D0,vu,np,dm,ld, & ! gamma, exponential variables dmin_mm,dmax_mm,ahp,bhp, & ! power law variables rg,log_sigma_g, & ! lognormal variables rho_e ! particle density (kg m^-3) real*8 :: tmp1, tmp2 real*8 :: pi,rc integer k,lidx,uidx pi = acos(-1.0) ! // if density is constant, store equivalent values for apm and bpm if ((rho_c > 0) .and. (apm < 0)) then apm = (pi/6)*rho_c bpm = 3. endif select case(dtype) ! ---------------------------------------------------------! ! // modified gamma ! ! ---------------------------------------------------------! ! :: N0 = total number concentration (m^-3) ! :: np = fixed number concentration (kg^-1) ! :: D0 = characteristic diameter (um) ! :: dm = mean diameter (um) ! :: vu = distribution width parameter case(1) if (abs(p1+1) < 1E-8) then ! // D0, vu are given vu = p3 if(Re.le.0) then dm = p2 D0 = gamma(vu)/gamma(vu+1)*dm else D0 = 2.0*Re*gamma(vu+2)/gamma(vu+3) endif if (scaled .eqv. .false.) then fc = ( & ((D*1E-6)**(vu-1)*exp(-1*D/D0)) / & (apm*((D0*1E-6)**(vu+bpm))*gamma(vu+bpm)) & ) * 1E-12 scaled = .true. endif N = fc*rho_a*(Q*1E-3) elseif (abs(p2+1) < 1E-8) then ! // N0, vu are given np = p1 vu = p3 tmp1 = (Q*1E-3)**(1./bpm) if (scaled .eqv. .false.) then fc = (D*1E-6 / (gamma(vu)/(apm*np*gamma(vu+bpm)))** & (1./bpm))**vu scaled = .true. endif N = ( & (rho_a*np*fc*(D*1E-6)**(-1.))/(gamma(vu)*tmp1**vu) * & exp(-1.*fc**(1./vu)/tmp1) & ) * 1E-12 else ! // vu isn't given print *, 'Error: Must specify a value for vu' stop endif ! ---------------------------------------------------------! ! // exponential ! ! ---------------------------------------------------------! ! :: N0 = intercept parameter (m^-4) ! :: ld = slope parameter (um) case(2) if (abs(p1+1) > 1E-8) then ! // N0 has been specified, determine ld N0 = p1 if(Re>0) then ! if Re is set and No is set than the distribution is fully defined. ! so we assume Re and No have already been chosen consistant with ! the water content, Q. ! print *,'using Re pass ...' ld = 1.5/Re ! units 1/um N = ( & N0*exp(-1*ld*D) & ) * 1E-12 else tmp1 = 1./(1.+bpm) if (scaled .eqv. .false.) then fc = ((apm*gamma(1.+bpm)*N0)**tmp1)*(D*1E-6) scaled = .true. endif N = ( & N0*exp(-1.*fc*(1./(rho_a*Q*1E-3))**tmp1) & ) * 1E-12 endif elseif (abs(p2+1) > 1E-8) then ! // ld has been specified, determine N0 ld = p2 if (scaled .eqv. .false.) then fc = (ld*1E6)**(1.+bpm)/(apm*gamma(1+bpm))* & exp(-1.*(ld*1E6)*(D*1E-6))*1E-12 scaled = .true. endif N = fc*rho_a*(Q*1E-3) else ! // ld will be determined from temperature, then N0 follows ld = 1220*10.**(-0.0245*tc)*1E-6 N0 = ((ld*1E6)**(1+bpm)*Q*1E-3*rho_a)/(apm*gamma(1+bpm)) N = ( & N0*exp(-1*ld*D) & ) * 1E-12 endif ! ---------------------------------------------------------! ! // power law ! ! ---------------------------------------------------------! ! :: ahp = Ar parameter (m^-4 mm^-bhp) ! :: bhp = br parameter ! :: dmin_mm = lower bound (mm) ! :: dmax_mm = upper bound (mm) case(3) ! :: br parameter if (abs(p1+2) < 1E-8) then ! :: if p1=-2, bhp is parameterized according to Ryan (2000), ! :: applicatable to cirrus clouds if (tc < -30) then bhp = -1.75+0.09*((tc+273)-243.16) elseif ((tc >= -30) .and. (tc < -9)) then bhp = -3.25-0.06*((tc+273)-265.66) else bhp = -2.15 endif elseif (abs(p1+3) < 1E-8) then ! :: if p1=-3, bhp is parameterized according to Ryan (2000), ! :: applicable to frontal clouds if (tc < -35) then bhp = -1.75+0.09*((tc+273)-243.16) elseif ((tc >= -35) .and. (tc < -17.5)) then bhp = -2.65+0.09*((tc+273)-255.66) elseif ((tc >= -17.5) .and. (tc < -9)) then bhp = -3.25-0.06*((tc+273)-265.66) else bhp = -2.15 endif else ! :: otherwise the specified value is used bhp = p1 endif ! :: Ar parameter dmin_mm = dmin*1E-3 dmax_mm = dmax*1E-3 ! :: commented lines are original method with constant density ! rc = 500. ! (kg/m^3) ! tmp1 = 6*rho_a*(bhp+4) ! tmp2 = pi*rc*(dmax_mm**(bhp+4))*(1-(dmin_mm/dmax_mm)**(bhp+4)) ! ahp = (Q*1E-3)*1E12*tmp1/tmp2 ! :: new method is more consistent with the rest of the distributions ! :: and allows density to vary with particle size tmp1 = rho_a*(Q*1E-3)*(bhp+bpm+1) tmp2 = apm*(dmax_mm**bhp*dmax**(bpm+1)-dmin_mm**bhp*dmin**(bpm+1)) ahp = tmp1/tmp2 * 1E24 ! ahp = tmp1/tmp2 lidx = infind(D,dmin) uidx = infind(D,dmax) do k=lidx,uidx N(k) = ( & ahp*(D(k)*1E-3)**bhp & ) * 1E-12 enddo ! print *,'test=',ahp,bhp,ahp/(rho_a*Q),D(100),N(100),bpm,dmin_mm,dmax_mm ! ---------------------------------------------------------! ! // monodisperse ! ! ---------------------------------------------------------! ! :: D0 = particle diameter (um) case(4) if (scaled .eqv. .false.) then D0 = p1 rho_e = (6/pi)*apm*(D0*1E-6)**(bpm-3) fc(1) = (6./(pi*D0**3*rho_e))*1E12 scaled = .true. endif N(1) = fc(1)*rho_a*(Q*1E-3) ! ---------------------------------------------------------! ! // lognormal ! ! ---------------------------------------------------------! ! :: N0 = total number concentration (m^-3) ! :: np = fixed number concentration (kg^-1) ! :: rg = mean radius (um) ! :: log_sigma_g = ln(geometric standard deviation) case(5) if (abs(p1+1) < 1E-8) then ! // rg, log_sigma_g are given log_sigma_g = p3 tmp2 = (bpm*log_sigma_g)**2. if(Re.le.0) then rg = p2 else rg =Re*exp(-2.5*(log_sigma_g**2)) endif if (scaled .eqv. .false.) then fc = 0.5 * ( & (1./((2.*rg*1E-6)**(bpm)*apm*(2.*pi)**(0.5) * & log_sigma_g*D*0.5*1E-6)) * & exp(-0.5*((log(0.5*D/rg)/log_sigma_g)**2.+tmp2)) & ) * 1E-12 scaled = .true. endif N = fc*rho_a*(Q*1E-3) elseif (abs(p2+1) < 1E-8) then ! // N0, log_sigma_g are given Np = p1 log_sigma_g = p3 N0 = np*rho_a tmp1 = (rho_a*(Q*1E-3))/(2.**bpm*apm*N0) tmp2 = exp(0.5*bpm**2.*(log_sigma_g))**2. rg = ((tmp1/tmp2)**(1/bpm))*1E6 N = 0.5*( & N0 / ((2.*pi)**(0.5)*log_sigma_g*D*0.5*1E-6) * & exp((-0.5*(log(0.5*D/rg)/log_sigma_g)**2.)) & ) * 1E-12 else ! // vu isn't given print *, 'Error: Must specify a value for sigma_g' stop endif end select end subroutine dsd