module module_mp_jensen_ishmael !------------------------------------------------------------------------------------------------------! ! This is the initial release of the Ice-Spheroids Habit Model with Aspect-ratio Evolution (ISHMAEL). ! ! This code calculates the moisture tendencies of water vapor, cloud water, rainwater, and three ice ! ! species which are planar-nucleated particles, columnar-nucleated particles, and aggregates. Planar- ! ! and columnar-nucleated ice particles are parameterized as spheroids (two growth dimensions are ! ! tracked), and ice particle shape and density evolution are predicted for vapor growth, sublimation, ! ! riming, and melting. The model can produce various vapor-grown habits (which currently depend on the ! ! environmental temperature and supersaturation), and the model predicts the change in ice particle ! ! shape and density during riming. Ice particle fall speeds are diagnosed based on ice particle ! ! properties. Aggregates are currently parameterized as spheroids with a fixed shape and density. Ice ! ! is initiated as planar or columnar based on temperature. After nucleation, the shapes of particles ! ! are not constrained, i.e., planar-nucleated particles can attain columnar shapes if the growth ! ! conditions allow it. ! ! ! ! Predicting two ice particle dimensions requires two prognostic variables. Thus mass, number, and two ! ! size-related mixing ratio are predicted for each ice species. That makes 12 prognostic variables ! ! for ice. Just throw more processors at it. Rain mass and number and cloud mass are also predicted ! ! See Jensen et al. 2017 (JAS), Jensen et al. 2018 (MWR), and Jensen et al. 2018 (JAS) for details. ! ! ! ! Author: Anders A. Jensen, NCAR-RAL, ajensen@ucar.edu, 303-497-8484 ! ! Last modified: 07 June 2019 (ISHMAEL version 0.1) ! !------------------------------------------------------------------------------------------------------! use module_wrf_error implicit none !.. You are welcome !.. Constants real, private, parameter :: & PI = 3.14159265, & !.. Pi = 3.14159265... until forever G_HOME = 9.8, & !.. Gravity (m s^-1) RD = 287.15, & !.. Dry air gas constant (J kg^-1 K^-1) RV = 461.5, & !.. Moist air gas constant (J kg^-1 K^-1) CP = 1004.0, & !.. Heat capacity Air (J kg^-1 K^-1) CPW = 4218.0, & !.. Heat capacity Water (J kg^-1 K^-1) RCP = 0.285856574, & !.. RD / CP RHOW = 1000.0, & !.. Water density (kg m^-3) T0 = 273.15, & !.. STP temperature (K) P0 = 101325.0, & !.. STP pressure (Pa) RHOI = 920.0, & !.. Bulk ice density (kg m^-3) R0 = 1.27494 !.. Air density at 1000 mb, T=T0 (kg m^-3) !.. The lookup table arrays are itab for ice-cloud water collection !.. and itabr for ice-rain collection. The lookup tables are read in !.. the subroutine jensen_ishmael_init which occurs once. real, private, allocatable, dimension(:,:,:,:,:) :: itab, itabr !.. Inherent growth ratio data is stored in the igrdata array !.. The gamma_tab is a tabulation of the gamma function needed !.. because ice particle shape evolves. real, private, allocatable, dimension(:) :: igrdata, gamma_tab !.. Aggregation lookup tables arrays real, private, allocatable, dimension(:,:,:) :: coltab, coltabn !.. Added for double precision integer, parameter :: LUT_KIND_R4 = selected_real_kind(6) !.. 4 byte real real(kind=LUT_KIND_R4), private, allocatable, dimension(:,:,:,:,:) :: itab_o, itabr_o real(kind=LUT_KIND_R4), private, allocatable, dimension(:) :: gamma_tab_o contains !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. This subroutine currently reads in the ice-cloud collection, !.. ice-rain collecion, gamma calculations, and inherent growth ratio !.. data lookup tables and calculates the ice-ice aggregation lookup table. subroutine jensen_ishmael_init implicit none integer :: i1, i2, i3, i4, icat integer, parameter :: ncat=7,npair=35,ndn=60 integer :: ipair(ncat,ncat) real, dimension(ncat,9) :: dstprms real, dimension(ncat) :: gnu,cfmas,pwmas,cfvt,pwvt,tablo,tabhi !.. Collection pairs ipair=& reshape((/0,0,3,7,11,15,19 ,1,0,5,9,13,17,21 & ,2,4,22,24,0,0,0 ,6,8,23,28,0,0,0 & ,10,12,25,29,35,0,0 ,14,16,26,30,32,0,0 & ,18,20,27,31,33,34,0/), & (/ncat,ncat/)) !.. Ice-ice collection lookup table properties (JYH) dstprms=& !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! cloud rain planar columnar aggreg graup hail ! 1 2 3 4 5 6 7 ---> ncat !----------------------------------------------------------------------- reshape((/.5, .5, .318, .318, .5, .5, .5, & !.. shape (capacitance shape) 524., 524., .333, .333, .496, 157., 471., & !.. cfmas (m=cfmas*D**pwmas) 3., 3., 2.4, 2.4, 2.4, 3., 3., & !.. pwmas 3173., 149., 4.836, 4.836, 3.084, 93.3, 161., & !.. cfvt (v = cfvt*D**pwvt) 2., .5, 0.25, .25, .2, .5, .5, & !.. pwvt 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, & !.. tablo (size limits for table) 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, & !.. tabhi .1e-5, .1e-4, .1e-5, .1e-5, .1e-4, .1e-4, .1e-4, & !.. dnmin (characteristic diam .1e-2, .1e-1, .1e-1, .1e-1, .1e-1, .1e-1, .1e-1/), & !.. dnmax limits) (/ncat,9/)) do icat=1,ncat cfmas(icat)=dstprms(icat,2) pwmas(icat)=dstprms(icat,3) cfvt (icat)=dstprms(icat,4) pwvt (icat)=dstprms(icat,5) tablo(icat)=dstprms(icat,6) tabhi(icat)=dstprms(icat,7) gnu(icat) = 4.0 enddo !.. Read in ice-cloud water collection lookup table if(.not.allocated(itab_o)) then allocate(itab_o(51,51,51,11,2)) open(20, FILE='ishmael-qi-qc.bin', FORM='unformatted') do i1 = 1,51 do i2 = 1,51 do i3 = 1,51 do i4 = 1,11 read(20) itab_o(i1,i2,i3,i4,1), & itab_o(i1,i2,i3,i4,2) enddo enddo enddo enddo close(20) print*, 'Jensen_ISHMAEL lookup table 1 of 3' endif !.. Read in ice-rain collection lookup table if(.not.allocated(itabr_o)) then allocate(itabr_o(51,51,51,11,6)) open(30, FILE='ishmael-qi-qr.bin', FORM='unformatted') do i1 = 1,51 do i2 = 1,51 do i3 = 1,51 do i4 = 1,11 read(30) itabr_o(i1,i2,i3,i4,1), & itabr_o(i1,i2,i3,i4,2), & itabr_o(i1,i2,i3,i4,3), & itabr_o(i1,i2,i3,i4,4), & itabr_o(i1,i2,i3,i4,5), & itabr_o(i1,i2,i3,i4,6) enddo enddo enddo enddo close(30) print*, 'Jensen_ISHMAEL lookup table 2 of 3' endif !.. Read in gamma function lookup table for speedy code if(.not.allocated(gamma_tab_o)) then allocate(gamma_tab_o(505001)) open(40, FILE='ishmael-gamma-tab.bin', FORM='unformatted') do i1 = 1,505001 read(40) gamma_tab_o(i1) enddo close(40) print*, 'Jensen_ISHMAEL lookup table 3 of 3' endif !.. Read in inhernent growth ratio data (-1C to -60C) !.. This currently assumes columns below -20C (plates are used later in the code) !.. From Chen and Lamb 1994 if(.not.allocated(igrdata)) then allocate(igrdata(60)) igrdata = & (/ 0.910547, 0.81807, 0.6874, 0.60127, 1.59767, 2.32423, 2.08818, & 1.61921, 1.15865, 0.863071, 0.617586, 0.453917, 0.351975, 0.28794, & 0.269298, 0.28794, 0.333623, 0.418883, 0.56992, 0.796458, 1.14325, & 1.64103, 1.90138, 1.82653, 1.61921, 1.47436, 1.32463, 1.25556, & 1.22239, 1.206, 1.11522, 1.10751, 1.10738, 1.11484, 1.12234, & 1.12221, 1.14529, 1.16884, 1.20104, 1.22573, 1.25094, 1.27666, & 1.31183, 1.3388, 1.35704, 1.37553, 1.38479, 1.39411, 1.40349, & 1.41294, 1.42245, 1.43202, 1.44166, 1.45137, 1.46114, 1.47097, & 1.48087, 1.50105, 1.50087, 1.51098 /) endif !.. make collection table (JYH) !.. called only on initialization of the model !.. coltab = mass mixing ratio collection !.. coltabn = number mixing ratio collection !.. ndn = table elements !.. ncat = # of hydrometeor categories if(.not.allocated(coltab)) then allocate(coltab(ndn,ndn,npair)) endif if(.not.allocated(coltabn)) then allocate(coltabn(ndn,ndn,npair)) endif !.. Call to build the aggregation lookup table (JYH) call mkcoltb(ndn,ncat,coltab,coltabn,ipair,gnu,tablo,tabhi,cfmas,pwmas,cfvt,pwvt) print*, 'Jensen_ISHMAEL aggregation lookup table built' !.. Check that tables are allocated (Thanks to JYH and Shaofeng) if(.not.allocated(itab)) allocate(itab(51,51,51,11,2)) if(.not.allocated(itabr)) allocate(itabr(51,51,51,11,6)) if(.not.allocated(gamma_tab)) allocate(gamma_tab(505001)) itab = itab_o; itabr = itabr_o; gamma_tab = gamma_tab_o end subroutine jensen_ishmael_init !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. The main subroutine. subroutine MP_JENSEN_ISHMAEL(ITIMESTEP, DT_IN, P, DZ, & TH, QV, QC, QR, NR, QI1, NI1, AI1, CI1, & QI2, NI2, AI2, CI2, & QI3, NI3, AI3, CI3, & IDS,IDE, JDS,JDE, KDS,KDE, & IMS,IME, JMS,JME, KMS,KME, & ITS,ITE, JTS,JTE, KTS,KTE, & RAINNC, RAINNCV, SNOWNC, SNOWNCV, & diag_effc3d, diag_effi3d, diag_dbz3d, & diag_vmi3d_1, diag_di3d_1, diag_rhopo3d_1, diag_phii3d_1, & diag_vmi3d_2, diag_di3d_2, diag_rhopo3d_2, diag_phii3d_2, & diag_vmi3d_3, diag_di3d_3, diag_rhopo3d_3, diag_phii3d_3, & diag_itype_1,diag_itype_2,diag_itype_3 & ) !.. ITIMESTEP - time step counter (integer) !.. DT_IN - model time step (s) !.. P - air pressure (Pa) !.. DZ - difference in height over interface (m) !.. TH - potential temperature (K) !.. QV - water vapor mixing ratio (kg kg^-1) !.. QC - cloud water mixing ratio (kg kg^-1) !.. QR - rain water mixing ratio (kg kg^-1) !.. NR - rain number concentration (kg^-1) !.. QI1 - planar-nucleated mixing ratio (kg kg^-1) !.. NI1 - planar-nucleated number mixing ratio (kg^-1) !.. AI1 - planar-nucleated a^2c mixing ratio (m^3 kg^-1) !.. CI1 - planar-nucleated c^2a mixing ratio (m^3 kg^-1) !.. QI2 - columnar-nucleated mixing ratio (kg kg^-1) !.. NI2 - columnar-nucleated number mixing ratio (kg^-1) !.. AI2 - columnar-nucleated a^2c mixing ratio (m^3 kg^-1) !.. CI2 - columnar-nucleated c^2a mixing ratio (m^3 kg^-1) !.. QI3 - aggregate mixing ratio (kg kg^-1) !.. NI3 - aggregate number mixing ratio (kg^-1) !.. AI3 - aggregate a^2c mixing ratio (m^3 kg^-1) !.. CI3 - aggregate c^2a mixing ratio (m^3 kg^-1) !.. IDS,IDE - I domain start/end locations !.. JDS,JDE - J domain start/end locations !.. KDS,KDE - K domain start/end locations !.. IMS,IME - I memory start/end locations !.. JMS,JME - J memory start/end locations !.. KMS,KME - K memory start/end locations !.. ITS,ITE - I tile start/end locations !.. JTS,JTE - J tile start/end locations !.. KTS,KTE - K tile start/end locations !.. RAINNC - total accumulated precipitation (mm) !.. RAINNCV - total accumulated precipitation in a time step (mm) !.. SNOWNC - total ice accumulated precipitation (mm) !.. SNOWNCV - total ice accumulated precipitation in a time step (mm) !.. DIAG_EFFC3D - Effective cloud water radius for coupling to radiation (m) !.. DIAG_EFFI3D - Effective ice radius for coupling to radiation (m) !.. DIAG_DBZ3D - 10-cm reflectivity (dBZ) !.. DIAG_VMI3D_1 - planar-nucleated mass-weighted fall speeds (m s^-1) !.. DIAG_VMI3D_2 - columnar-nucleated mass-weighted fall speeds (m s^-1) !.. DIAG_VMI3D_3 - aggregate mass-weighted fall speeds (m s^-1) !.. DIAG_DI3D_1 - planar-nucleated mass-weighted maximum diameter (m) !.. DIAG_DI3D_2 - columnar-nucleated mass-weighted maximum diameter (m) !.. DIAG_DI3D_3 - aggregate mass-weighted maximum diameter (m) !.. DIAG_RHOPO3D_1 - planar-nucleated mass-weighted effective density (kg m^-3) !.. DIAG_RHOPO3D_2 - columnar-nucleated mass-weighted effective density (kg m^-3) !.. DIAG_RHOPO3D_3 - aggregate mass-weighted effective density (kg m^-3) !.. DIAG_PHII3D_1 - planar-nucleated number-weighted aspect ratio !.. DIAG_PHII3D_2 - columnar-nucleated number-weighted aspect ratio !.. DIAG_PHII3D_3 - aggregate number-weighted aspect ratio !.. DIAG_ITYPE_1 - currently not used !.. DIAG_ITYPE_2 - currently not used !.. DIAG_ITYPE_3 - currently not used implicit none !.. Grid variables integer, intent(in) :: & ids, ide, jds, jde, kds, kde , & ims, ime, jms, jme, kms, kme , & its, ite, jts, jte, kts, kte !.. Time and timestep integer, intent(in) :: ITIMESTEP real, intent(in) :: dt_in !.. 3D variables real, dimension(ims:ime, kms:kme, jms:jme), intent(in) :: p, dz real, dimension(ims:ime, kms:kme, jms:jme), intent(inout) :: & th, qv, qc, qr, nr, qi1, ni1, ai1, ci1, qi2, ni2, ai2, ci2, & qi3, ni3, ai3, ci3, & diag_effc3d, diag_effi3d, diag_dbz3d, & diag_vmi3d_1, diag_di3d_1, diag_rhopo3d_1, diag_phii3d_1, & diag_vmi3d_2, diag_di3d_2, diag_rhopo3d_2, diag_phii3d_2, & diag_vmi3d_3, diag_di3d_3, diag_rhopo3d_3, diag_phii3d_3, & diag_itype_1, diag_itype_2, diag_itype_3 !.. 2D variables real, dimension(ims:ime, jms:jme), intent(inout) :: & rainnc, rainncv, snownc, snowncv integer :: I,K,J !.. Define the ice species !.. Do not change these without changing/knowing how nucleation works integer, parameter :: CAT = 3 !.. Number of ice categories (species) integer, parameter :: ICE1 = 1 !.. Ice initiated as planar integer, parameter :: ICE2 = 2 !.. Ice initiated as columnar integer, parameter :: ICE3 = 3 !.. Aggregates !.. Use tile locations for our variables !.. Make everthing 1D for our loop !.. These variables go in and out of newmicro real, dimension(kts:kte) :: & th1d, qv1d, qc1d, qr1d, nr1d, p1d, dz1d, & effc1d, effi1d, dbz1d real, dimension(cat, kts:kte) :: & qi1d, ni1d, ai1d, ci1d, & vmi1d, di1d, rhopo1d, phii1d, icetype1d !.. Precipitation real :: qrpre1d real, dimension(cat) :: qipre1d !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Accumulated precipitation in one timestep (zero out before microphysics) qrpre1d =0. qipre1d(ICE1)=0. qipre1d(ICE2)=0. qipre1d(ICE3)=0. !.. Loop over the grid and send each vertical column into the microphysics routine do j=jts,jte !.. j loop (north-south) do i=its,ite !.. i loop (horizontal) !.. Get the 1D variables from WRF and zero out process rates do k=kts,kte !.. k loop (vertical) !.. Thermo,rain,cloud,vapor,etc TH1D(k) = TH(i,k,j) QV1D(k) = QV(i,k,j) QC1D(k) = QC(i,k,j) QR1D(k) = QR(i,k,j) NR1D(k) = NR(i,k,j) P1D(k) = P(i,k,j) DZ1D(k) = DZ(i,k,j) EFFC1D(k) = 0. EFFI1D(k) = 0. DBZ1D(k) = -35. !.. Planar-nucleated QI1D(ICE1,k) = QI1(i,k,j) NI1D(ICE1,k) = NI1(i,k,j) AI1D(ICE1,k) = AI1(i,k,j) CI1D(ICE1,k) = CI1(i,k,j) VMI1D(ICE1,k) = 0. DI1D(ICE1,k) = 0. RHOPO1D(ICE1,k) = 0. PHII1D(ICE1,k) = 0. ICETYPE1D(ICE1,k) = 0. !.. Columnar-nucleated QI1D(ICE2,k) = QI2(i,k,j) NI1D(ICE2,k) = NI2(i,k,j) AI1D(ICE2,k) = AI2(i,k,j) CI1D(ICE2,k) = CI2(i,k,j) VMI1D(ICE2,k) = 0. DI1D(ICE2,k) = 0. RHOPO1D(ICE2,k) = 0. PHII1D(ICE2,k) = 0. ICETYPE1D(ICE2,k) = 0. !.. Aggregates QI1D(ICE3,k) = QI3(i,k,j) NI1D(ICE3,k) = NI3(i,k,j) AI1D(ICE3,k) = AI3(i,k,j) CI1D(ICE3,k) = CI3(i,k,j) VMI1D(ICE3,k) = 0. DI1D(ICE3,k) = 0. RHOPO1D(ICE3,k) = 0. PHII1D(ICE3,k) = 0. ICETYPE1D(ICE3,k) = 0. enddo !.. (vertical k-loop) !.. Call me Ishmael. call me_ishmael(itimestep, kts, kte, i, j, dt_in, p1d, dz1d, & th1d, qv1d, qc1d, qr1d, nr1d, qi1d, ni1d, ai1d, ci1d, & CAT, ICE1, ICE2, ICE3, effc1d, effi1d, dbz1d, icetype1d, & vmi1d, di1d, rhopo1d, phii1d, qrpre1d, qipre1d) !.. Send the 1D columns from newmicro back to WRF do k=kts,kte !.. k loop (vertical) TH(i,k,j) = TH1D(k) QV(i,k,j) = QV1D(k) QC(i,k,j) = QC1D(k) QR(i,k,j) = QR1D(k) NR(i,k,j) = NR1D(k) DIAG_EFFC3D(i,k,j) = max((min(EFFC1D(k),50.e-6)),2.51e-6) DIAG_EFFI3D(i,k,j) = max((min(EFFI1D(k),999.e-6)),5.e-6) DIAG_DBZ3D(i,k,j) = max(dbz1d(k),-35.) !.. Planar-nucleated QI1(i,k,j) = QI1D(ICE1,k) NI1(i,k,j) = NI1D(ICE1,k) AI1(i,k,j) = AI1D(ICE1,k) CI1(i,k,j) = CI1D(ICE1,k) DIAG_VMI3D_1(i,k,j) = vmi1d(ICE1,k) DIAG_DI3D_1(i,k,j) = di1d(ICE1,k) DIAG_RHOPO3D_1(i,k,j)= rhopo1d(ICE1,k) DIAG_PHII3D_1(i,k,j) = phii1d(ICE1,k) DIAG_ITYPE_1(i,k,j) = icetype1d(ICE1,k) !.. Columnar-nucleated QI2(i,k,j) = QI1D(ICE2,k) NI2(i,k,j) = NI1D(ICE2,k) AI2(i,k,j) = AI1D(ICE2,k) CI2(i,k,j) = CI1D(ICE2,k) DIAG_VMI3D_2(i,k,j) = vmi1d(ICE2,k) DIAG_DI3D_2(i,k,j) = di1d(ICE2,k) DIAG_RHOPO3D_2(i,k,j)= rhopo1d(ICE2,k) DIAG_PHII3D_2(i,k,j) = phii1d(ICE2,k) DIAG_ITYPE_2(i,k,j) = icetype1d(ICE2,k) !.. Aggregates QI3(i,k,j) = QI1D(ICE3,k) NI3(i,k,j) = NI1D(ICE3,k) AI3(i,k,j) = AI1D(ICE3,k) CI3(i,k,j) = CI1D(ICE3,k) DIAG_VMI3D_3(i,k,j) = vmi1d(ICE3,k) DIAG_DI3D_3(i,k,j) = di1d(ICE3,k) DIAG_RHOPO3D_3(i,k,j)= rhopo1d(ICE3,k) DIAG_PHII3D_3(i,k,j) = phii1d(ICE3,k) DIAG_ITYPE_3(i,k,j) = icetype1d(ICE3,k) enddo !.. Accumulated precipitation and precipitation rates !.. Frozen precipitation SNOWNC(i,j) = SNOWNC(i,j) + qipre1d(ICE1) + qipre1d(ICE2) + qipre1d(ICE3) SNOWNCV(i,j) = qipre1d(ICE1) + qipre1d(ICE2) + qipre1d(ICE3) !.. Total precipitation RAINNC(i,j) = RAINNC(i,j) + qrpre1d + qipre1d(ICE1) + qipre1d(ICE2) + qipre1d(ICE3) RAINNCV(i,j) = qrpre1d + qipre1d(ICE1) + qipre1d(ICE2) + qipre1d(ICE3) enddo ! i loop (horizontal) enddo ! j loop (north-south) end subroutine MP_JENSEN_ISHMAEL !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. The microphysics subroutine me_ishmael(it, kts, kte, i, j, dt, pres_e, dzmic, & theta, qv, qc, qr, nr, qi, ni, ai, ci, CAT, ICE1, ICE2, ICE3, & effc1d, effi1d, dbz1d, icetype1d, vmi1d, di1d, rhopo1d, phii1d, & qrpre1d, qipre1d) implicit none !.. User logicals (all should be TRUE except for testing) logical, parameter :: RIMING = .TRUE. ! Riming logical, parameter :: RAIN_ICE = .TRUE. ! Rain-ice collection (for testing) logical, parameter :: FREEZE_QC = .TRUE. ! Freeze cloud drops (Homogeneous) logical, parameter :: SPLINTERS = .TRUE. ! Rime Splinters logical, parameter :: WET_GROW = .TRUE. ! Wet growth !.. Input/Output variables from WRF integer, intent(in) :: it, i, j !.. Integer timestep, model i and j integer, intent(in) :: kts, kte !.. Tile k-start and k-end integer, intent(in) :: CAT, ICE1, ICE2, ICE3 !.. Number of categories, planar and columnar index real, intent(in) :: dt !.. Timestep real, dimension(kts:kte), intent(in) :: pres_e, dzmic !.. Pressure and dz real, dimension(kts:kte), intent(inout) :: theta !.. Theta (K) real, dimension(kts:kte), intent(inout) :: qv !.. Vapor mixing ratio (kg kg^-1) real, dimension(kts:kte), intent(inout) :: qc !.. Cloud water mixing ratio (kg kg^-1) real, dimension(kts:kte), intent(inout) :: qr !.. Rain mixing ratio (kg kg^-1) real, dimension(kts:kte), intent(inout) :: nr !.. Rain number mixing ratio (# kg^-1) real, dimension(cat, kts:kte), intent(inout) :: qi !.. Ice mass mixing ratio (kg kg^-1) real, dimension(cat, kts:kte), intent(inout) :: ni !.. Ice number mixing ratio (# kg^-1) real, dimension(cat, kts:kte), intent(inout) :: ai !.. Ice a^2c mixing ratio (m^3 kg^-1) real, dimension(cat, kts:kte), intent(inout) :: ci !.. Ice c^2a mixing ratio (m^3 kg^-1) !.. Ice diagnostics real, dimension(kts:kte), intent(inout) :: effc1d, effi1d, dbz1d real, dimension(cat, kts:kte), intent(inout) :: rhopo1d real, dimension(cat, kts:kte), intent(inout) :: phii1d real, dimension(cat, kts:kte), intent(inout) :: vmi1d real, dimension(cat, kts:kte), intent(inout) :: di1d real, dimension(cat, kts:kte), intent(inout) :: icetype1d real :: dbzr,dbzsum real :: qrpre1d !.. Rain precipitation accumulation (mm) real, dimension(cat) :: qipre1d !.. Ice liquid-equivalent precipitation accumulation (mm) !.. The Parameters real, parameter :: NU = 4.0 !.. Ice distribution parameter (keep at 4) !.. If changed, build new lookup tables real, parameter :: QSMALL = 1.e-12 !.. Smallest ice mass (kg kg^-1) real, parameter :: QNSMALL= 1.25e-7 !.. Smallest ice number (# kg^-1) real, parameter :: QASMALL= 1.e-24 !.. Smallest ice volume (m^-3 kg^-1) real, parameter :: ao = 0.1e-6 !.. alphstr=ao^(1-deltastr) real, parameter :: AR = 149.1 !.. 'a' parameter in rain fallspeed-size relation (m^(1-br)/s) real, parameter :: BR = 0.5 !.. 'b' parameter in rain fallspeed-size relation real, parameter :: F1R = 0.78 !.. Ventilation coefficient for rain real, parameter :: F2R = 0.308 !.. Ventilation coefficient for rain real, parameter :: LAMMAXR = 1./20.E-6 !.. Max allowed value for rain slope parameter (m^-1) real, parameter :: LAMMINR = 1./2800.E-6 !.. Min allowed value for rain slope parameter (m^-1) !.. Local/Dummy variables logical :: vgflag !.. Vapor growth flag to determine if vapor growth !.. is producing subsaturated conditions logical, dimension(cat) :: has_ice !.. true if ice exists for an ice species logical, dimension(cat) :: dry_growth !.. Wet growth/dry growth riming flag integer :: k, cc, ccvar !.. k is the k-loop index, cc is the ice species index integer :: gi !.. Index for the Gamma(NU) lookup integer :: gi2, gi3 !.. Index for the Gamma(NU) lookup integer :: numice !.. Number of ice species currently with ice integer :: current_index !.. Used to determine number of ice species wi th mass > QSMALL real, allocatable, dimension(:) :: icearray !.. Array with size of number !.. of ice species with mass > QSMALL real :: wetg !.. wet growth check real :: gamma_arg !.. Gamma argument: Gamma(gamma_arg) real, dimension(cat) :: ani, cni, rni !.. Characteristic a-,c-,r-axis sizes real, dimension(cat) :: aniold, cniold, dsold real :: alphstr !.. alpha_* = ao^(1-delta_*) real :: alphv !.. alphv = 4/3*PI*alpha_* real :: betam !.. betam = 2 + delta_* real :: anf, cnf, rnf, iwcf !.. a-,c-,r-axis and iwc after vapor growth real :: anfr, cnfr, rnfr, iwcfr !.. a-,c-,r-axis and iwc after riming real :: phibr, phifr !.. aspect ratio before and after riming real :: capgam !.. distribution-weighted capacitance (m) real :: i_dt !.. inverse timestep real :: i_cp !.. inverse heat capacity of air real :: sink, source, ratio !.. Terms used to prevent over-depletion real :: rrr !.. rain characteristic radius !.. Fall speed parameters real :: phiivt, bl, al, aa, ba, qe, xn, bx, xm, f_c1, f_c2, bm, am, nre !.. Local variables for riming real :: vir, vif, vfr, rhobarrime, rimedr, rimedrr, rimec1, rhobarwg !.. Local variables for ice nucleation, multiplication real :: anuc, nibnuc(CAT), amass, cmass, fmult, nucfrac real :: a_demott, b_demott, c_demott, d_demott real :: curnum, ratel, ratekg, inrate !.. Local variables for vapor growth real :: vi, iwci, nim3dum, prdsum, qcrimesum, qcrimesumr !.. Local variables for melting real :: mfact, nimelt !.. Local variables for evaporation real :: epsr, dumqv, dumqvs, arn !.. Some dummy variables real :: dum, dum1, dumt, dumgam, dumgam2 real :: tmpsum, tmpsumr, tmpmelt !.. Local variables for Radiation real :: lrsig, sig, lwc, ncm3dum, r_n, phirad, radw1, radw2, radw3, totaliwc !.. Local variables for Aggregation real :: dn1, dn2, dn3, phiagg1, phiagg2, phiagg3 !.. Fall speeds (some not used) and effectice ice radius variables real, dimension(kts:kte) :: vtrm, vtrn, vtrmc, vtrnc, nc real, dimension(cat, kts:kte) :: vtrni1, effi, vtrmi1, vtrzi1 !.. Lookup table values real :: rrri, rrni, rqci, rdsi, rrho, proc(cat,2), procr(cat,6) integer :: irri, irni, iqci, idsi, irho, iti real :: qi_qc_nrm(cat), qi_qc_nrd(cat) real :: qi_qr_nrm(cat), qi_qr_nrd(cat), qi_qr_nrn(cat) !.. Sedimentation variables integer :: nn, nstep real, dimension(kts:kte) :: fluxqr, fluxnr real, dimension(cat, kts:kte) :: fluxqi, fluxni, fluxai, fluxci real, dimension(cat) :: falltndqi, falltndni, falltndai, falltndci real :: falltndqr, falltndnr, maxfall logical :: sedi logical :: domicro !.. The rest of the variables real, dimension(kts:kte) :: rhoair, i_rhoair !.. Air density (kg m^-3) and inverse air density real :: qcrimefrac(cat) !.. Qc_Rime/(Qc_Rime + Qr_Rime) real :: gammnu !.. Gamma(NU) real :: i_gammnu !.. 1/Gamma(NU) real :: fourthirdspi !.. 4/3 * PI real :: mnuccd !.. Deposition nucleation rate (kg kg^-1 s^-1) real :: nnuccd !.. Deposition nucleation rate (# kg^-1) real :: pcc !.. Cloud mass growth rate (kg kg^-1 s^-1) real :: qagg(cat) !.. Aggregation mass tendency (kg kg^-1) real :: nagg(cat) !.. Aggregation number tendency (# kg^-1) real :: prd(cat) !.. Vapor growth rate (kg kg^-1 s^-1) real :: nrd(cat) !.. Vapor growth number rate sublimation (# kg^-1 s^-1) real :: ard(cat) !.. A-axis vapor growth rate mixing ratio (m kg^-1 s^-1) real :: crd(cat) !.. C-axis vapor growth rate mixing ratio (m kg^-1 s^-1) real :: prdr(cat) !.. Riming growth rate (kg kg^-1 s^-1) real :: ardr(cat) !.. A-axis riming growth rate mixing ratio (m kg^-1 s^-1) real :: crdr(cat) !.. C-axis riming growth rate mixing ratio (m kg^-1 s^-1) real :: qmlt(cat) !.. Mass melting rate (kg kg^-1 s^-1) real :: nmlt(cat) !.. Number melting rate (# kg^-1 s^-1) real :: amlt(cat) !.. A-axis melting rate (m kg^-1 s^-1) real :: cmlt(cat) !.. C-axis melting rate (m kg^-1 s^-1) real :: rimesum(cat) !.. Ice-cloud rime mass rate (kg m^-3 s^-1) real :: rimesumr(cat) !.. Ice-rain rime mass rate (kg m^-3 s^-1) real :: rimetotal(cat) !.. Total rime mass rate (kg m^-3 s^-1) real :: dQRfzri(cat) !.. Change in rain mass from ice-rain (kg kg^-1 s^-1) real :: dQIfzri(cat) !.. Change in ice mass from ice-rain (kg kg^-1 s^-1) real :: dNfzri(cat) !.. Change in ice number from ice-rain (# kg^-1 s^-1) real :: dQImltri(cat) !.. Change in ice mass from ice-rain melt (kg kg^-1 s^-1) real :: dNmltri(cat) !.. Change in ice number from ice-rain melt (# kg^-1 s^-1) real :: numrateri(cat) !.. Number rate from ice-rain (value from lookkup table) real :: rainrateri(cat) !.. Rain mass rate from ice-rain (value from lookkup table) real :: icerateri(cat) !.. Ice mass rate from ice-rain (value from lookkup table) real :: totalnuc(cat) !.. Total mass nucleated to ice (kg kg^-1 s^-1) real :: totalnucn(cat) !.. Total number nucleated to ice (# kg^-1 s^-1) real :: nim !.. Number of homogeneously frozen cloud drops (# kg^-1 s^-1) real :: mim !.. Mass of homogeneously frozen cloud drops (kg kg^-1 s^-1) real :: nimr !.. Number of homogeneously frozen rain drops (# kg^-1 s^-1) real :: mimr !.. Mass of homogeneously frozen rain drops (kg kg^-1 s^-1) real :: deltastr(cat) !.. Delta_* (Average inherent growth ratio) real :: rhobar(cat) !.. Ice effective density (kg m^-3) real :: rhodepout(cat) !.. Ice effective density (kg m^-3) real :: dsdepout(cat) !.. Ice effective density (kg m^-3) real :: rhorimeout(cat) !.. Ice effective density (kg m^-3) real :: gdentotal !.. Total rime density from ice-cloud and ice-rain (kg m^-3) real :: gdenavg(cat) !.. Rime density from ice-cloud (kg m^-3) real :: gdenavgr(cat) !.. Rime density from ice-rain (kg m^-3) real :: fv(cat) !.. Ventilation coefficient for diffusion real :: fh(cat) !.. Ventilation coefficient for heat real :: temp, i_temp !.. Air temperature (K) and inverse real :: qvs !.. Liquid saturation vapor mixing ratio (kg kg^-1) real :: qvi !.. Ice saturation vapor mixing ratio (kg kg^-1) real :: sup !.. Liquid supersaturation ratio real :: sui !.. Ice supersaturation ratio real :: svpi !.. Ice saturation vapor pressure real :: svpl !.. Liquid saturation vapor pressure real :: qs0 !.. Saturation vapor mixing ratio at 0C (kg kg^-1) real :: xxls !.. Latent heat of sublimation (J kg^-1) real :: xxlv !.. Latent heat of vaporization (J kg^-1) real :: xxlf !.. Latent heat of fusion (J kg^-1) real :: ab !.. Latent heat correction parameter for growth real :: abi !.. Latent heat correction parameter for growth (ice) real :: mu !.. Dynamic visocity of air (kg m^-1 s^-1) real :: dv !.. Diffusivity of water vapor in air (m^2 s^-1) real :: kt !.. Thermal conductivity (W m^-1 K^-1) real :: nsch !.. Schmidt number real :: npr !.. Prandlt number real :: lamr !.. Rain slope parameter n(D)=n0rr*exp(-lamr*D)dD real :: n0rr !.. Rain distribution parameter real :: prc !.. Autoconversion cloud water to rain (kg kg^-1 s^-1) real :: nprc !.. Autoconversion cloud number loss (# kg^-1 s^-1) real :: nprc1 !.. Autoconversion rain number gain (# kg^-1 s^-1) real :: pre !.. Rain evpaoration rate (kg kg^-1 s^-1) real :: npre !.. Rain evaporation number (# kg^-1 s^-1) real :: pra !.. Rain accretion of cloud (kg kg^-1 s^-1) real :: npra !.. Rain accretion of cloud number (# kg^-1 s^-1) real :: nragg !.. Rain self-collection (# kg^-1 s^-1) real :: igr !.. Inherent growth ratio real :: qmult(cat) !.. Ice multiplication rate (kg kg^-1 s^-1) real :: nmult(cat) !.. Ice number multiplication rate (# kg^-1 s^-1) real :: mbiggr !.. Biggs freezing (kg kg^-1 s^-1) real :: nbiggr !.. Biggs number freezing (# kg^-1 s^-1) !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. NU=4 is constant so gammnu=6 gammnu = gamma(NU) i_gammnu = 1./gammnu !.. Inverse gamma i_dt = 1./dt !.. Inverse timestep i_cp = 1./CP !.. Inverse heat capacity fourthirdspi = 4./3.*PI !.. Four divided by three times pi !.. No need to loop through microphysics and sedimentation unless !.. condensate exists sedi = .false. !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! do k=kts,kte !.. Grid K loop !.. Check the column for condensate or supersaturation wrt liquid. !.. If it is there, then do the microphysics domicro= .false. temp = theta(k)/(100000./pres_e(k))**(RCP) svpl = polysvp(temp,0) qvs = 0.622*svpl/(pres_e(k)-svpl) sup = qv(k)/qvs-1. svpi = polysvp(temp,1) qvi = 0.622*svpi/(pres_e(k)-svpi) if(temp.gt.T0) qvi=qvs sui = qv(k)/qvi-1. if(qc(k).gt.QSMALL.or.qr(k).gt.QSMALL.or.qi(ICE1,k).gt.QSMALL.or.& qi(ICE2,k).gt.QSMALL.or.qi(ICE3,k).gt.QSMALL.or.sup.ge.0.) then domicro= .true. endif !.. If condensate or supersaturation wrt liquid, do microphysics if(domicro) then if(allocated(icearray)) deallocate(icearray) numice = 0 !.. Number of species with ice do cc = 1, cat !.. Loop over all ice species has_ice(cc) = .false. if(qi(cc,k).gt.QSMALL) then !.. Check for ice !.. Get the number of species with ice has_ice(cc) = .true. numice = numice + 1 endif enddo !.. Determine the number of ice species to loop over !.. For vapor growth/riming, etc if(numice.gt.0.and..not.allocated(icearray)) then allocate(icearray(numice)) current_index = 1 do cc = 1, cat if(has_ice(cc)) then icearray(current_index) = cc current_index = current_index + 1 endif enddo endif !.. Initialize variables for all ice species and process rates do cc = 1, cat !.. Ice variables dry_growth(cc)=.true. qcrimefrac(cc)=0. rimetotal(cc)=0. nmult(cc)=0. qmult(cc)=0. ani(cc)=0. cni(cc)=0. rni(cc)=0. qagg(cc)=0. nagg(cc)=0. prd(cc) =0. nrd(cc) =0. ard(cc) =0. crd(cc) =0. prdr(cc)=0. ardr(cc)=0. crdr(cc)=0. qmlt(cc)=0. nmlt(cc)=0. amlt(cc)=0. cmlt(cc)=0. deltastr(cc)=1. dsold(cc)=1. dsdepout(cc)=1. rhobar(cc)=RHOI rhodepout(cc)=RHOI rhorimeout(cc)=RHOI if(cc.eq.3) rhobar(cc)=50. fv(cc)=1. fh(cc)=1. rimesum(cc)=0. rimesumr(cc)=0. numrateri(cc)=0. rainrateri(cc)=0. icerateri(cc)=0. dQRfzri(cc)=0. dQIfzri(cc)=0. dNfzri(cc)=0. dQImltri(cc)=0. dNmltri(cc)=0. gdenavg(cc)=400. gdenavgr(cc)=400. totalnuc(cc)=0. totalnucn(cc)=0. nibnuc(cc)=0. aniold(cc)=0. cniold(cc)=0. proc(cc,1)=0. proc(cc,2)=0. qi_qc_nrm(cc)=0. qi_qc_nrd(cc)=0. qi_qr_nrm(cc)=0. qi_qr_nrn(cc)=0. qi_qr_nrd(cc)=0. procr(cc,1)=0. procr(cc,2)=0. procr(cc,3)=0. procr(cc,4)=0. procr(cc,5)=0. procr(cc,6)=0. fluxqi(cc,k)=0. fluxni(cc,k)=0. fluxai(cc,k)=0. fluxci(cc,k)=0. falltndqi(cc)=0. falltndni(cc)=0. falltndai(cc)=0. falltndci(cc)=0. di1d(cc,k)=0. vmi1d(cc,k)=0. vtrmi1(cc,k)=0. vtrni1(cc,k)=0. vtrzi1(cc,k)=0. rhopo1d(cc,k)=0. phii1d(cc,k)=0. qipre1d(cc)=0. effi(cc,k)=0. icetype1d(cc,k)=0. enddo dbz1d(k)=-35. qrpre1d=0. fluxqr(k)=0. fluxnr(k)=0. falltndqr=0.; falltndnr=0. vtrm(k)=0.; vtrn(k)=0. effc1d(k) = 0. effi1d(k) = 0. nprc =0.; npre=0.; npra=0.; pra=0.; prc=0. nprc1=0.; nragg=0. nim =0.; mim =0. nimr=0.; mimr=0. mbiggr=0.; nbiggr=0. mnuccd=0.; nnuccd=0.; pcc=0.; fmult=0. gdentotal=400. rimedr=0.; rimedrr=0. dbzr=0.; dbzsum=0. !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Temperature and temperature-dependent variables i_temp= 1./temp rhoair(k)= pres_e(k)/(RD*temp) i_rhoair(k) = 1./rhoair(k) svpi = polysvp(temp,1) qvi = 0.622*svpi/(pres_e(k)-svpi) if(temp.gt.T0) qvi=qvs qs0 = 0.622*polysvp(T0,0)/(pres_e(k)-polysvp(T0,0)) sui = qv(k)/qvi-1. xxls = 3.15e6-2370.*temp+0.3337e6 xxlv = 3.1484E6-2370.*temp xxlf = xxls - xxlv ab = 1.+(xxlv*qvs/(RV*temp**2))*xxlv*i_cp abi = 1.+(xxls*qvi/(RV*temp**2))*xxls*i_cp mu = 1.496E-6*temp**1.5/(temp+120.) arn = ar*(R0*i_rhoair(k))**0.5 dv = 8.794E-5*temp**1.81/pres_e(k) kt = 2.3823e-2 +7.1177e-5*(temp-T0) nsch = mu*i_rhoair(k)/dv npr = mu*i_rhoair(k)/kt !.. Inherent growth ratio data from Lamb and Scott 1972 and Chen and Lamb 1994 !.. Assumes planar particles at temperatures below -20C (see Bailey and Hallett) !.. Function of temperature only (thanks to JHY and GHB for moving this fuction to where it belongs) igr = get_igr(igrdata, temp) if((temp-T0).lt.-20.) then igr=0.7 endif !.. Assume constant cloud drop number concentration (200 cm^-3) nc(k) = 200.e6*i_rhoair(k) !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Ice Microphysics if(numice.gt.0) then !.. If one or more ice species have ice tmpsum = 0. tmpsumr= 0. do ccvar = 1, numice cc = icearray(ccvar) !.. Variable check ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) !.. Set incoming volume limit after advection !.. This is the smallest bulk volume for which we care about shape !.. If smaller than this limit, assume spherical if(ai(cc,k).lt.1.e-12.or.ci(cc,k).lt.1.e-12) then ai(cc,k)=min(ai(cc,k),ci(cc,k)) ci(cc,k)=ai(cc,k) ani(cc)=(ai(cc,k)/ni(cc,k))**0.3333333333 cni(cc)=(ci(cc,k)/ni(cc,k))**0.3333333333 endif if(ani(cc).lt.2.e-6.or.cni(cc).lt.2.e-6) then ani(cc)=2.e-6 cni(cc)=2.e-6 endif ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k) ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k) deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) !.. Check to keep ice properties in bounds call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) dsold(cc) = deltastr(cc) !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Ice-cloud riming rates from lookup table if(RIMING.and.qc(k).gt.1.e-7) then !.. The lookup table is in log space. Do not change the following 4 lines without !.. changing the offline-built lookup table rrni = 13.498*log10(0.5e6*rni(cc)) rqci = 8.776*log10(1.0e7*(exp(qc(k))-1.)) rdsi = 50.*(deltastr(cc) - 0.5) rrho = 7.888*log10(0.02*rhobar(cc)) irni = int(rrni) iqci = int(rqci) idsi = int(rdsi) irho = int(rrho) !.. Small limit (stay in the bounds of the lookup table) rrni = max(rrni,1.) rqci = max(rqci,1.) rdsi = max(rdsi,1.) rrho = max(rrho,1.) irni = max(irni,1) iqci = max(iqci,1) idsi = max(idsi,1) irho = max(irho,1) !.. Large limit (stay in the bounds of the lookup table) rrni = min(rrni,real(size(itab,1))-1.) rqci = min(rqci,real(size(itab,2))-1.) rdsi = min(rdsi,real(size(itab,3))-1.) rrho = min(rrho,real(size(itab,4))-1.) irni = min(irni,size(itab,1)-1) iqci = min(iqci,size(itab,2)-1) idsi = min(idsi,size(itab,3)-1) irho = min(irho,size(itab,4)-1) do iti = 1, 2 call access_lookup_table(itab,irni,iqci,idsi,irho,iti, & rdsi,rrho,rqci,rrni,proc(cc,iti)) enddo !.. proc(cc,1) is the normalized riming rate !.. proc(cc,2) is related to the normalized rime density !.. rimesum is the riming rate in kg m^-3 s^-1 rimesum(cc) = max((proc(cc,1)*ni(cc,k)*nc(k)*rhoair(k)**2),0.) qi_qc_nrm(cc)= proc(cc,1) qi_qc_nrd(cc)= proc(cc,2) !.. Limit riming rate if too small if((rimesum(cc)*i_rhoair(k)).lt.QSMALL) then rimesum(cc) =0. qi_qc_nrm(cc)=0. qi_qc_nrd(cc)=0. endif endif !.. End ice-cloud riming !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Ice-rain riming rates from lookup table if(RAIN_ICE.and.qr(k).gt.QSMALL) then nr(k) = max(nr(k),QNSMALL) !.. Get rain size distribution properties lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333 n0rr = nr(k)*lamr if(lamr.LT.LAMMINR) then lamr = LAMMINR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr elseif(lamr.gt.LAMMAXR) then lamr = LAMMAXR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr endif rrr = 0.5 * (1./lamr) !.. The lookup table is in log space. Do not change the following 4 lines without !.. changing the offline-built lookup table rrni = 13.498*log10(0.5e6*rni(cc)) rrri= 23.273*log10(1.e5*rrr) rdsi = 50.*(deltastr(cc) - 0.5) rrho = 7.888*log10(0.02*rhobar(cc)) irni = int(rrni) irri = int(rrri) idsi = int(rdsi) irho = int(rrho) !.. Small lookup table limit rrni = max(rrni,1.) rrri = max(rrri,1.) rdsi = max(rdsi,1.) rrho = max(rrho,1.) irni = max(irni,1) irri = max(irri,1) idsi = max(idsi,1) irho = max(irho,1) !.. Large lookup table limit rrni = min(rrni,real(size(itabr,1))-1.) rrri = min(rrri,real(size(itabr,2))-1.) rdsi = min(rdsi,real(size(itabr,3))-1.) rrho = min(rrho,real(size(itabr,4))-1.) irni = min(irni,size(itabr,1)-1) irri = min(irri,size(itabr,2)-1) idsi = min(idsi,size(itabr,3)-1) irho = min(irho,size(itabr,4)-1) do iti = 1, 6 call access_lookup_table(itabr,irni,irri,idsi,irho,iti, & rdsi,rrho,rrri,rrni,procr(cc,iti)) enddo !.. procr(cc,4) is the normalized ice number tendency from ice-rain !.. procr(cc,5) is the normalized rain mass tendency from ice-rain !.. procr(cc,6) is the normalized ice mass tendency from ice-rain numrateri(cc) = max((procr(cc,4)*ni(cc,k)*nr(k)*rhoair(k)),0.) rainrateri(cc)= max((procr(cc,5)*ni(cc,k)*nr(k)*rhoair(k)),0.) icerateri(cc) = max((procr(cc,6)*ni(cc,k)*nr(k)*rhoair(k)),0.) !.. Limit ice-rain rates if to small if(rainrateri(cc).lt.QSMALL.or.icerateri(cc).lt.QSMALL) then rainrateri(cc)=0. icerateri(cc) =0. numrateri(cc) =0. endif if(temp.le.T0) then !.. For temperatures less than melting point ice-rain freezes rain if(qr(k).gt.0.1e-3.and.qi(cc,k).gt.0.1e-3) then dQRfzri(cc) = rainrateri(cc) dQIfzri(cc) = icerateri(cc) dNfzri(cc) = numrateri(cc) endif else !.. For temperatures greater than melting point ice-rain melts ice dQImltri(cc) = icerateri(cc) dNmltri(cc) = numrateri(cc) endif !.. procr(cc,1) is the normalized riming rate !.. procr(cc,2) is trelated to the he normalized rime density !.. procr(cc,3) is the normalized loss in rain number from riming !.. rimesumr is in kg m^-3 s^-1 !.. Ice-rain leads to riming when the fall speed of the collected rain is slower than !.. the fall speed of the ice collector. Otherwise, ice-rain freezes rain rimesumr(cc) = max((procr(cc,1)*ni(cc,k)*nr(k)*rhoair(k)**2),0.) qi_qr_nrm(cc)= procr(cc,1) qi_qr_nrd(cc)= procr(cc,2) qi_qr_nrn(cc)= procr(cc,3) !.. Limit ice-rain riming if small if((rimesumr(cc)*i_rhoair(k)).lt.QSMALL) then rimesumr(cc) =0. qi_qr_nrm(cc)=0. qi_qr_nrd(cc)=0. qi_qr_nrn(cc)=0. endif endif !.. End ice-rain collection !.. Do not over-deplete cloud water or rain from ice-cloud and ice-rain tmpsum = tmpsum + rimesum(cc) tmpsumr= tmpsumr+ rimesumr(cc) if(ccvar.eq.numice) then !.. Check at end !.. Do not let ice-cloud collection overdeplete cloud water sink = tmpsum*i_rhoair(k) source = qc(k)*i_dt if(sink.gt.source.and.qc(k).gt.QSMALL) then ratio = source / sink rimesum = rimesum*ratio qi_qc_nrm= qi_qc_nrm*ratio qi_qc_nrd= qi_qc_nrd*ratio endif !.. Do not let ice-rain collection overdeplete rainwater sink = tmpsum*i_rhoair(k) source = qr(k)*i_dt if(sink.gt.source.and.qr(k).gt.QSMALL) then ratio = source / sink rimesumr = rimesumr*ratio qi_qr_nrm = qi_qr_nrm*ratio qi_qr_nrd = qi_qr_nrd*ratio qi_qr_nrn = qi_qr_nrn*ratio endif endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Inital volume, iwc, #/m^3 for vapor growth gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) vi = fourthirdspi*rni(cc)**3*gamma_tab(gi)*i_gammnu iwci = ni(cc,k)*rhobar(cc)*vi*rhoair(k) nim3dum=ni(cc,k)*rhoair(k) !.. Total riming rate (the heating goes into vapor growth) rimetotal(cc) = rimesum(cc)+rimesumr(cc) alphstr=ao**(1.-deltastr(cc)) !.. Distribtuion-averaged capacitance (see Harrington et al. 2013) capgam = capacitance_gamma(ani(cc), deltastr(cc), NU, alphstr, i_gammnu) !.. This subroutine first calculates the bulk fall speeds, then the ventilation, !.. and finally the vapor growth (or sublimation) rate. (See Harrington et al. 2013) call vaporgrow(dt, ani(cc), cni(cc), rni(cc), igr, nim3dum, temp, rimetotal(cc), & pres_e(k), NU, alphstr, sui, sup, qvs, qvi, mu, iwci, rhoair(k), qi(cc,k), & dv, kt, ao, nsch, npr, gammnu, i_gammnu, fourthirdspi, svpi, xxls, xxlv, xxlf, & capgam, vtrni1(cc,k), vtrmi1(cc,k), vtrzi1(cc,k), anf, cnf, rnf, iwcf, fv(cc), & fh(cc), rhobar(cc), deltastr(cc),rhodepout(cc),dsdepout(cc)) !.. Vapor growth mass growth and axes rates prd(cc)=(iwcf-iwci)*i_rhoair(k)*i_dt prd(cc)=max(prd(cc),-qi(cc,k)*i_dt) ard(cc)=(2.*(anf-ani(cc))*cni(cc)+(cnf-cni(cc))*ani(cc))*ani(cc)*ni(cc,k)*i_dt crd(cc)=(2.*(cnf-cni(cc))*ani(cc)+(anf-ani(cc))*cni(cc))*cni(cc)*ni(cc,k)*i_dt !.. Sublimation number loss if(prd(cc).lt.0..and.qi(cc,k).gt.QSMALL) then nrd(cc)=prd(cc)*ni(cc,k)/qi(cc,k) endif !.. Limit vapor growth if values small if(abs(prd(cc)*dt).lt.(QSMALL*0.01)) then prd(cc)=0. nrd(cc)=0. ard(cc)=0. crd(cc)=0. endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Now calculated the mass growth and axes growth rates and density for riming !.. Check for wet-growth riming conditions !.. See Lamb and Verlinde (2011) if(WET_GROW) then rimetotal(cc)=rimesum(cc)+rimesumr(cc) call wet_growth_check(NU, temp, rhoair(k), xxlv, xxlf, qv(k), dv, kt, qs0, & fv(cc), fh(cc), rimetotal(cc), rni(cc), ni(cc,k), dry_growth(cc)) endif !.. If above the melting point, assume high-density riming if(temp.gt.T0) then dry_growth(cc) = .false. endif !.. Initial rime volume and r-axis length vir = vi rnfr= rni(cc) !.. Calculate the ice-cloud rime density !.. This is parameterized from Macklin (1962) and is described in !.. Jensen et al. 2017 (JAS) if(RIMING.and.qc(k).gt.QSMALL.and.qi_qc_nrm(cc).gt.0.) then rimec1 = 0.0066 !.. Value at -10C if((temp-T0).lt.-30.) then rimec1 = 0.0036 elseif((temp-T0).lt.-20.and.(temp-T0).ge.-30.) then dum = (abs((temp-T0))-20.) / 10. rimec1 = dum*(0.0036) + (1.-dum)*(0.004) elseif((temp-T0).lt.-15.and.(temp-T0).ge.-20.) then dum = (abs((temp-T0))-15.) / 5. rimec1 = dum*(0.004) + (1.-dum)*(0.005) elseif((temp-T0).lt.-10.and.(temp-T0).ge.-15.) then dum = (abs((temp-T0))-10.) / 5. rimec1 = dum*(0.005) + (1.-dum)*(0.0066) elseif((temp-T0).lt.-5.and.(temp-T0).ge.-10.) then dum = (abs((temp-T0))-5.) / 5. rimec1 = dum*(0.0066) + (1.-dum)*(0.012) elseif((temp-T0).ge.-5.) then rimec1 = 0.012 endif gdenavg(cc) = (1000.*(0.8*tanh(rimec1*qi_qc_nrd(cc)/qi_qc_nrm(cc))+0.1)) if((temp-T0).gt.-5..and.(temp-T0).le.0.) then dum = (abs((temp-T0))-0.) / 5. gdenavg(cc) = dum*gdenavg(cc) + (1.-dum)*900. endif if((temp-T0).gt.0..or..NOT.dry_growth(cc)) then gdenavg(cc)=900. endif !.. Keep rime density between 50 and 900 kg m^-3 gdenavg(cc) = max(gdenavg(cc),50.) gdenavg(cc) = min(gdenavg(cc),900.) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) !.. Change in r-axis from riming (based on the rime density) rimedr = max(((qi_qc_nrm(cc)/gdenavg(cc))/((gamma_tab(gi)*i_gammnu)* & 4.*PI*rni(cc)**2)),0.) rnfr = rnfr + max((rimedr*nc(k)*rhoair(k)),0.)*dt endif !.. End ice-cloud rime density calculation !.. Calculate the ice-rain rime density !.. This is parameterized from Macklin (1962) if(RAIN_ICE.and.qr(k).gt.QSMALL.and.qi_qr_nrm(cc).gt.0.) then rimec1 = 0.0066 !.. Value at -10C if((temp-T0).lt.-30.) then rimec1 = 0.0036 elseif((temp-T0).lt.-20.and.(temp-T0).ge.-30.) then dum = (abs((temp-T0))-20.) / 10. rimec1 = dum*(0.0036) + (1.-dum)*(0.004) elseif((temp-T0).lt.-15.and.(temp-T0).ge.-20.) then dum = (abs((temp-T0))-15.) / 5. rimec1 = dum*(0.004) + (1.-dum)*(0.005) elseif((temp-T0).lt.-10.and.(temp-T0).ge.-15.) then dum = (abs((temp-T0))-10.) / 5. rimec1 = dum*(0.005) + (1.-dum)*(0.0066) elseif((temp-T0).lt.-5.and.(temp-T0).ge.-10.) then dum = (abs((temp-T0))-5.) / 5. rimec1 = dum*(0.0066) + (1.-dum)*(0.012) elseif((temp-T0).ge.-5.) then rimec1 = 0.012 endif gdenavgr(cc) = (1000.*(0.8*tanh(rimec1*qi_qr_nrd(cc)/qi_qr_nrm(cc))+0.1)) if((temp-T0).gt.-5..and.(temp-T0).le.0.) then dum = (abs((temp-T0))-0.) / 5. gdenavgr(cc) = dum*gdenavgr(cc) + (1.-dum)*900. endif if((temp-T0).gt.0..or..NOT.dry_growth(cc)) then gdenavgr(cc)=900. endif !.. Keep rime density between 50 and 900 kg m^03 gdenavgr(cc) = max(gdenavgr(cc),50.) gdenavgr(cc) = min(gdenavgr(cc),900.) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) !.. Change in r-axs from riming (based on rime density) rimedrr = max(((qi_qr_nrm(cc)/gdenavgr(cc))/((gamma_tab(gi)*i_gammnu)* & 4.*PI*rni(cc)**2)),0.) rnfr = rnfr + max((rimedrr*nr(k)*rhoair(k)),0.)*dt endif !.. End ice-rain rime density calculation !.. Caculate the volume after riming gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) vfr = fourthirdspi*rnfr**3*gamma_tab(gi)*i_gammnu vfr = max(vfr,vi) rnfr= max(rnfr,rni(cc)) !.. Determine the fraction of rime mass from cloud versus small rain drops !.. Update the ice density and iwc if(rimetotal(cc).gt.0.and.vfr.gt.vi) then qcrimefrac(cc) = rimesum(cc)/rimetotal(cc) qcrimefrac(cc) = max(qcrimefrac(cc),0.) qcrimefrac(cc) = min(qcrimefrac(cc),1.) gdentotal = qcrimefrac(cc)*gdenavg(cc) + (1.-qcrimefrac(cc))*gdenavgr(cc) rhobarrime = rhobar(cc) rhorimeout(cc) = (rhobar(cc)*(vi/vfr))+(gdentotal*(1.-(vi/vfr))) rhorimeout(cc) = min(rhorimeout(cc),RHOI) iwcfr = rhorimeout(cc) * vfr * nim3dum rhorimeout(cc) = gdentotal rhorimeout(cc) = min(rhorimeout(cc),RHOI) else iwcfr= iwci rnfr = rni(cc) endif !.. For dry growth the drops freeze on contace and one axis (the axis perpendicular to the !.. fall direction) grows if(dry_growth(cc)) then cnfr =cni(cc) anfr =ani(cc) phibr= 1. if(rimetotal(cc).gt.0.and.vfr.gt.vi) then gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) phibr = cni(cc)/ani(cc)*gamma_tab(gi)*i_gammnu if(phibr.eq.1.) then !.. Spherical ice phifr= 1. anfr = rnfr cnfr = rnfr elseif(phibr.gt.1.25) then !.. Prolate ice phifr = phibr*(((rnfr/rni(cc))**3)**(-0.5)) anfr = (ani(cc)/rni(cc)**(1.5))*rnfr**(1.5) elseif(phibr.lt.0.8) then !.. Oblate ice phifr = phibr*(rnfr/rni(cc))**3 gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) cnfr = phifr*anfr*gammnu/gamma_tab(gi) else phifr = phibr !.. Constant aspect ratio if quasi-spherical gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+2.+deltastr(cc) gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) anfr = ((vfr*gamma_tab(gi))/(fourthirdspi*phifr*gamma_tab(gi2)))**0.3333333333 cnfr = phifr*anfr*gammnu/gamma_tab(gi) endif !.. For long time steps, riming can cause oblate ice to become prolate or vise versa !.. This prevents that from occuring if(phibr.le.1..and.phifr.gt.1.) then phifr = 0.99 gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+2.+deltastr(cc) gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) anfr = ((vfr*gamma_tab(gi))/(fourthirdspi*phifr*gamma_tab(gi2)))**0.3333333333 cnfr = phifr*anfr*gammnu/gamma_tab(gi) endif if(phibr.ge.1..and.phifr.lt.1.) then phifr = 1.01 gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+2.+deltastr(cc) gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) anfr = ((vfr*gamma_tab(gi))/(fourthirdspi*phifr*gamma_tab(gi2)))**0.3333333333 cnfr = phifr*anfr*gammnu/gamma_tab(gi) endif endif !.. End axis evolution from riming !.. Don't let the axes shrink from riming cnfr=max(cnfr,cni(cc)) anfr=max(anfr,ani(cc)) !.. Riming mass growth and axes rates prdr(cc)=(iwcfr-iwci)*i_rhoair(k)*i_dt ardr(cc)=(2.*(anfr-ani(cc))*cni(cc)+(cnfr-cni(cc))*ani(cc))*ani(cc)*ni(cc,k)*i_dt crdr(cc)=(2.*(cnfr-cni(cc))*ani(cc)+(anfr-ani(cc))*cni(cc))*cni(cc)*ni(cc,k)*i_dt !.. These can never be negative prdr(cc) = max(prdr(cc),0.) ardr(cc) = max(ardr(cc),0.) crdr(cc) = max(crdr(cc),0.) !.. If no rime mass, no axis growth if(prdr(cc).eq.0.0) then ardr(cc) = 0. crdr(cc) = 0. endif else !.. Wet Growth !.. For wet growth, mass is added but soaks into the particle and thus !.. the density increases but the size does not prdr(cc)=(iwcfr-iwci)*i_rhoair(k)*i_dt ardr(cc)=0. crdr(cc)=0. endif !.. End dry versus wet growth if(temp.gt.T0) then !.. Melting: ice shapes become more spherical qmlt(cc)=2.*PI*(kt*fh(cc)*(T0-temp)+rhoair(cc)*xxlv*dv*fv(cc)*(qs0-qv(k)))/ & xxlf*(ni(cc,k)*NU*max(ani(cc),cni(cc))) - & (CPW/xxlf*(temp-T0)*(rimetotal(cc)/rhoair(k) + dQImltri(cc))) qmlt(cc)=min(qmlt(cc),0.) qmlt(cc)=max(qmlt(cc),(-qi(cc,k)*i_dt)) !.. Don't let small very number mixing ratios (large sizes) cause ice precipitation if(qmlt(cc).lt.0.) then if(ai(cc,k).lt.1.e-12.or.ci(cc,k).lt.1.e-12) then qmlt(cc)=-qi(cc,k)*i_dt endif endif nmlt(cc) = max((-ni(cc,k)*i_dt),(ni(cc,k)*qmlt(cc)/qi(cc,k)-dNmltri(cc))) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) tmpmelt=fourthirdspi*alphstr*gamma_tab(gi)*i_gammnu amlt(cc) = (1./tmpmelt)*(1./rhobar(cc))*qmlt(cc) + ai(cc,k)*nmlt(cc)/ni(cc,k) cmlt(cc) = amlt(cc)*cni(cc)/ani(cc)*(1.+(2.*deltastr(cc)))/(2.+deltastr(cc)) endif enddo !.. End the loop over ice speices endif !.. End if ice !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Nucleation: this occurs outside of the ice species loop !.. Freeze evertying at temperatures below -35C !.. That is how it happens in nature, right? if(FREEZE_QC.and.qc(k).gt.QSMALL.and.temp.lt.(T0-35.)) then mim = qc(k)*i_dt nim = nc(k)*i_dt endif if(FREEZE_QC.and.qr(k).gt.QSMALL.and.temp.lt.(T0-35.)) then mimr = qr(k)*i_dt nimr = nr(k)*i_dt endif !.. Bigg (1953) freezing below -4C if(qr(k).gt.QSMALL.and.temp.lt.(T0-4.)) then nr(k) = max(nr(k),QNSMALL) lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333 n0rr = nr(k)*lamr if(lamr.LT.LAMMINR) then lamr = LAMMINR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr elseif(lamr.gt.LAMMAXR) then lamr = LAMMAXR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr endif mbiggr = 20.*PI*PI*RHOW*100.*nr(k)*(exp(0.66*(T0-temp))-1.)/lamr**3/lamr**3 nbiggr = PI*100.*nr(k)*(exp(0.66*(T0-temp))-1.)/lamr**3 mbiggr = min(mbiggr,qr(k)*i_dt) nbiggr = min(nbiggr,nr(k)*i_dt) endif !.. DeMott (2010) PNAS, Heterogeneous nuclation if(temp.lt.T0.and.sup.ge.0.) then a_demott = 0.0000594 b_demott = 3.33 c_demott = 0.0264 d_demott = 0.0033 !.. The 0.03 is the number of large aerosol (>0.5 microns) !.. See Chagnon and Junge 1962 inrate = a_demott * (273.16 - temp)**b_demott * & 0.03**((c_demott*(273.16 - temp)) + d_demott) inrate = inrate*1000.0 !.. #/L to #/m^3 inrate = inrate*i_rhoair(k) !.. # kg^-1 if((((ni(ICE1,k)+ni(ICE2,k)+inrate)*rhoair(k))/1000.).le.10000.) then mnuccd=(inrate)*(fourthirdspi*RHOI*(2.e-6)**3)*i_dt nnuccd=inrate*i_dt else curnum = (ni(ICE1,k)+ni(ICE2,k))*rhoair(k)/1000. ratel = max(0.,(10000.-curnum)) ratekg = ratel*1000.*i_rhoair(k) mnuccd=(ratekg)*(fourthirdspi*RHOI*(2.e-6)**3)*i_dt nnuccd=ratekg*i_dt endif endif !.. Rime splinter (Hallett and Mossop 1974, Nature) if(SPLINTERS) then if(temp.lt.270.16.and.temp.gt.265.16) then if (temp.GT.270.16) then fmult = 0. elseif(temp.le.270.16.and.temp.gt.268.16) then fmult = (270.16-temp)/2. elseif(temp.ge.265.16.and.temp.le.268.16) then fmult = (temp-265.16)/3. elseif(temp.lt.265.16) then fmult = 0. endif endif do cc = 1, cat if(prdr(cc).gt.0.) then nmult(cc) = 35.E4*prdr(cc)*fmult*1000. qmult(cc) = nmult(cc)*fourthirdspi*RHOI*(5.e-6)**3 qmult(cc) = min(qmult(cc),prdr(cc)) prdr(cc) = prdr(cc)-qmult(cc) endif enddo endif !.. Aggregation (currently only for T < T0) if(temp.le.T0) then do cc = 1, cat ni(cc,k) = max(ni(cc,k),QNSMALL) ai(cc,k) = max(ai(cc,k),QASMALL) ci(cc,k) = max(ci(cc,k),QASMALL) enddo !.. Get characteristic diameters for aggregation dn1 = 2.*((ai(ICE1,k)**2)/(ci(ICE1,k)*ni(ICE1,k)))**0.333333333333 dn2 = 2.*((ci(ICE2,k)**2)/(ai(ICE2,k)*ni(ICE2,k)))**0.333333333333 dn3 = 2.*((ai(ICE3,k)**2)/(ci(ICE3,k)*ni(ICE3,k)))**0.333333333333 dn1 = MIN(dn1,1.e-2) dn2 = MIN(dn2,1.e-2) dn3 = MIN(dn3,1.e-2) dn1 = MAX(dn1,1.e-6) dn2 = MAX(dn2,1.e-6) dn3 = MAX(dn3,1.e-6) !.. Get ice-one and ice-two aspect ratios !.. Spherical paritcles aggregate at a slower rate than less spherical ones gamma_arg = NU-1.+deltastr(1) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) phiagg1=ci(1,k)/ai(1,k)*gamma_tab(gi)*i_gammnu gamma_arg = NU-1.+deltastr(2) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) phiagg2=ci(2,k)/ai(2,k)*gamma_tab(gi)*i_gammnu phiagg1=MIN(phiagg1,100.) phiagg1=MAX(phiagg1,0.01) phiagg2=MIN(phiagg2,100.) phiagg2=MAX(phiagg2,0.01) !.. Ice-1 and ice-2 mass and number loss to aggregates occurs inside of aggregation !.. After aggregation, ice-1 and ice-2 bulk distribution properties are adjusted !.. assuming that aggregation did not cause the original shape and density of the !.. unaggregated particles to change if(qi(ICE1,k).gt.1.e-8.or.qi(ICE2,k).gt.1.e-8.or.qi(ICE3,k).gt.1.e-8) then call aggregation(dt,rhoair(k),temp,qi(ICE1,k),ni(ICE1,k),dn1,qi(ICE2,k),ni(ICE2,k),dn2, & qi(ICE3,k),ni(ICE3,k),dn3,rhobar(ICE1),rhobar(ICE2),phiagg1,phiagg2,coltab, & coltabn,qagg(ICE1),qagg(ICE2),qagg(ICE3),nagg(ICE1),nagg(ICE2),nagg(ICE3)) endif !.. If ice mass mixing ratios are high enough to aggregate endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Warm Microphysics !.. Get raindrop size distribution parameters if (qr(k).gt.QSMALL) then nr(k) = max(nr(k),QNSMALL) lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333 n0rr = nr(k)*lamr if(lamr.LT.LAMMINR) then lamr = LAMMINR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr elseif(lamr.gt.LAMMAXR) then lamr = LAMMAXR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr endif endif !.. note: nprc1 is change in Nr, !.. nprc is change in Nc !.. NOTE: NPRC IS CALCULATED BUT NOT CURRENTLY USED SINCE NC IS NOT PROGNOSED !.. Auto conversion of cloud water to rain (Khairoutdinov and Kogan 2000) if(qc(k).ge.1.E-6) then PRC=1350.*qc(k)**2.47* & (nc(k)/1.e6*rhoair(k))**(-1.79) NPRC1 = PRC/(4./3.*PI*RHOW*(25.E-6)**3) NPRC = PRC/(qc(k)/nc(k)) NPRC = min(NPRC,nc(k)*i_dt) NPRC1 = min(NPRC1,NPRC) endif !.. Accretion of cloud water by rain Khairoutdinov and Kogan 2000) !.. NOTE: NPRA CURRENTLY NOT USED SINCE NC IS NOT PROGNOSED if(qr(k).ge.1.e-8.and.qc(k).ge.1.e-8) then DUM=(qc(k)*qr(k)) PRA = 67.*(DUM)**1.15 NPRA = PRA/(qc(k)/nc(k)) endif !.. Self-collection (Beheng 1994), Breakup (Verlinde and Cotton 1993) if(qr(k).ge.1.e-8) then dum1=300.e-6 if(1./lamr.lt.dum1) then dum=1. elseif(1./lamr.ge.dum1) then dum=2.-exp(2300.*(1./lamr-dum1)) endif nragg = -5.78*dum*nr(k)*qr(k)*rhoair(k) endif !.. Evaporation of rain (Rutledge and Hobbs 1983) if(qr(k).gt.QSMALL) then epsr = 2.*PI*n0rr*rhoair(k)*dv* & (F1R/(LAMR*LAMR)+ & F2R*(arn*rhoair(k)/mu)**0.5* & nsch**0.333333333*(gamma(5./2.+BR/2.))/ & (lamr**(5./2.+BR/2.))) else epsr = 0. endif if(qv(k).lt.qvs) then pre = epsr*(qv(k)-qvs)/ab pre = MIN(PRE,0.) dumt = temp+xxlv*i_cp*pre*dt dumqv = qv(k)-pre*dt dumqvs = 0.622*polysvp(dumt,0)/ & (pres_e(k)-polysvp(dumt,0)) if(pre.lt.0..and.dumqv/dumqvs.gt.1.) then pre=(qv(k)-qvs)*i_dt/ab endif else pre = 0. endif !.. Rain number mixing ratio loss from evaporation if(pre.lt.0.) then dum = pre*dt/qr(k) dum = max(-1.,dum) npre = dum*nr(k)/dt endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Updates after process-rate calculations !.. For ice, updates to mass, number, volume, and volume !.. times aspect ratio occur first for vapor growth, !.. riming, and melting (shape is predicted to evolve). !.. Next updates the bulk distribution shapes are done for !.. nucleation and aggregation. For nucleation, ice particle shape !.. itself does not change but the distribution shape becomes !.. more spherical as nucleated particles are added. !.. Do not over-deplete !.. Cloud water !.. First sum the total loss from riming qcrimesum=0. do cc = 1, cat qcrimesum = qcrimesum + prdr(cc)*qcrimefrac(cc) enddo sink =(prc+pra+qcrimesum+mim)*dt source=qc(k) if(sink.gt.source.and.qc(k).gt.QSMALL) then ratio = source/sink prc = prc*ratio pra = pra*ratio mim = mim*ratio do cc = 1, cat prdr(cc) = prdr(cc)*ratio ardr(cc) = ardr(cc)*ratio crdr(cc) = crdr(cc)*ratio enddo endif !.. Rain !.. First sum the total loss from riming qcrimesumr=0. do cc = 1, cat qcrimesumr = qcrimesumr + prdr(cc)*(1.-qcrimefrac(cc)) enddo sink = (-pre+qcrimesumr+mimr+mbiggr)*dt source = qr(k) + (pra+prc)*dt do cc = 1, cat source = source + (-qmlt(cc)*dt) enddo do cc = 1, cat sink = sink + (dQRfzri(cc))*dt enddo if(sink.gt.source.and.qr(k).gt.QSMALL) then ratio= source/sink do cc = 1, cat prdr(cc) = prdr(cc)*ratio ardr(cc) = ardr(cc)*ratio crdr(cc) = crdr(cc)*ratio qi_qr_nrn(cc)=qi_qr_nrn(cc)*ratio dQRfzri(cc) = dQRfzri(cc)*ratio dQIfzri(cc) = dQIfzri(cc)*ratio dNfzri(cc) = dNfzri(cc)*ratio enddo mimr = mimr*ratio nimr = nimr*ratio mbiggr =mbiggr*ratio nbiggr =nbiggr*ratio pre = pre*ratio endif !.. Ice if(igr.le.1.) then !.. Nucleation initiation to ice 1 !.. Ice 1 sink =(-qmlt(ICE1)-qagg(ICE1))*dt source=qi(ICE1,k) + (prdr(ICE1))*dt if(prd(ICE1).gt.0.) then source = source + prd(ICE1)*dt else sink = sink + (-prd(ICE1))*dt endif !.. Added because of ice initiation to ice 1 source = source + (mnuccd+mim+mimr+mbiggr+dQIfzri(ICE2)+dQIfzri(ICE3))*dt do cc = 1, cat source = source + (qmult(cc)+dQRfzri(cc))*dt enddo if(sink.gt.source.and.qi(ICE1,k).gt.QSMALL) then ratio= source/sink qmlt(ICE1) = qmlt(ICE1)*ratio qagg(ICE1) = qagg(ICE1)*ratio nagg(ICE1) = nagg(ICE1)*ratio if(prd(ICE1).lt.0.) then prd(ICE1) = prd(ICE1)*ratio endif endif !.. Ice 2 sink =(dQIfzri(ICE2)-qmlt(ICE2)-qagg(ICE2))*dt source=qi(ICE2,k) + (prdr(ICE2))*dt if(prd(ICE2).gt.0.) then source = source + prd(ICE2)*dt else sink = sink + (-prd(ICE2))*dt endif if(sink.gt.source.and.qi(ICE2,k).gt.QSMALL) then ratio= source/sink dQIfzri(ICE2) = dQIfzri(ICE2)*ratio dQRfzri(ICE2) = dQRfzri(ICE2)*ratio dNfzri(ICE2) = dNfzri(ICE2)*ratio qmlt(ICE2) = qmlt(ICE2)*ratio qagg(ICE2) = qagg(ICE2)*ratio nagg(ICE2) = nagg(ICE2)*ratio if(prd(ICE2).lt.0.) then prd(ICE2) = prd(ICE2)*ratio endif endif !.. Ice 3 sink =(dQIfzri(ICE3)-qmlt(ICE3))*dt source=qi(ICE3,k) + (prdr(ICE3)+qagg(ICE3))*dt if(prd(ICE3).gt.0.) then source = source + prd(ICE3)*dt else sink = sink + (-prd(ICE3))*dt endif if(sink.gt.source.and.qi(ICE3,k).gt.QSMALL) then ratio= source/sink dQIfzri(ICE3) = dQIfzri(ICE3)*ratio dQRfzri(ICE3) = dQRfzri(ICE3)*ratio dNfzri(ICE3) = dNfzri(ICE3)*ratio qmlt(ICE3) = qmlt(ICE3)*ratio if(prd(ICE3).lt.0.) then prd(ICE3) = prd(ICE3)*ratio endif endif else !.. nucleation to ice 2 !.. Ice 2 sink =(-qmlt(ICE2)-qagg(ICE2))*dt source=qi(ICE2,k) + (prdr(ICE2))*dt if(prd(ICE2).gt.0.) then source = source + prd(ICE2)*dt else sink = sink + (-prd(ICE2))*dt endif !.. Added because of ice initiation to ice 2 source = source + (mnuccd+mim+mimr+mbiggr+dQIfzri(ICE1)+dQIfzri(ICE3))*dt do cc = 1, cat source = source + (qmult(cc)+dQRfzri(cc))*dt enddo if(sink.gt.source.and.qi(ICE2,k).gt.QSMALL) then ratio= source/sink qmlt(ICE2) = qmlt(ICE2)*ratio qagg(ICE2) = qagg(ICE2)*ratio nagg(ICE2) = nagg(ICE2)*ratio if(prd(ICE2).lt.0.) then prd(ICE2) = prd(ICE2)*ratio endif endif !.. Ice 1 sink =(dQIfzri(ICE1)-qmlt(ICE1)-qagg(ICE1))*dt source=qi(ICE1,k) + (prdr(ICE1))*dt if(prd(ICE1).gt.0.) then source = source + prd(ICE1)*dt else sink = sink + (-prd(ICE1))*dt endif if(sink.gt.source.and.qi(ICE1,k).gt.QSMALL) then ratio= source/sink dQIfzri(ICE1) = dQIfzri(ICE1)*ratio dQRfzri(ICE1) = dQRfzri(ICE1)*ratio dNfzri(ICE1) = dNfzri(ICE1)*ratio qmlt(ICE1) = qmlt(ICE1)*ratio qagg(ICE1) = qagg(ICE1)*ratio nagg(ICE1) = nagg(ICE1)*ratio if(prd(ICE1).lt.0.) then prd(ICE1) = prd(ICE1)*ratio endif endif !.. Ice 3 sink =(dQIfzri(ICE3)-qmlt(ICE3))*dt source=qi(ICE3,k) + (prdr(ICE3)+qagg(ICE3))*dt if(prd(ICE3).gt.0.) then source = source + prd(ICE3)*dt else sink = sink + (-prd(ICE3))*dt endif if(sink.gt.source.and.qi(ICE3,k).gt.QSMALL) then ratio= source/sink dQIfzri(ICE3) = dQIfzri(ICE3)*ratio dQRfzri(ICE3) = dQRfzri(ICE3)*ratio dNfzri(ICE3) = dNfzri(ICE3)*ratio qmlt(ICE3) = qmlt(ICE3)*ratio if(prd(ICE3).lt.0.) then prd(ICE3) = prd(ICE3)*ratio endif endif endif !.. Do not let vapor growth of ice create subsaturated conditions prdsum=0. do cc = 1, cat prdsum = prdsum + prd(cc) enddo vgflag = .false. sink = prdsum source=0.99*sui*qvi/abi*i_dt if(sink.gt.source.and.sui.gt.0.) then vgflag = .true. ratio = source/sink do cc = 1, cat prd(cc) = prd(cc)*ratio enddo endif !.. Update to ice do cc = 1, cat !.. Vapor growth update if((qi(cc,k)+(prd(cc))*dt).gt.QSMALL.and.qi(cc,k).gt.QSMALL) then rhobar(cc) = rhobar(cc)*(qi(cc,k)/(qi(cc,k)+(prd(cc))*dt)) + & rhodepout(cc)*(1.-(qi(cc,k)/(qi(cc,k)+(prd(cc))*dt))) else rhobar(cc) = RHOI endif qi(cc,k)=qi(cc,k)+(prd(cc))*dt ni(cc,k)=ni(cc,k)+(nrd(cc))*dt ai(cc,k)=ai(cc,k)+(ard(cc))*dt ci(cc,k)=ci(cc,k)+(crd(cc))*dt ni(cc,k) = max(ni(cc,k),QNSMALL) ai(cc,k) = max(ai(cc,k),QASMALL) ci(cc,k) = max(ci(cc,k),QASMALL) !.. If vapor growth of ice is going to cause !.. subsaturated condtions, reduce vapor growth !.. rates (done above) and rediagnose the axes if(vgflag.and.qi(cc,k).gt.QSMALL) then alphstr=ao**(1.-deltastr(cc)) alphv=fourthirdspi*alphstr betam=2.+deltastr(cc) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) ani(cc)=((qi(cc,k)*gammnu)/(rhobar(cc)*ni(cc,k)*alphv* & gamma_tab(gi)))**(1./betam) cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc) ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k) ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k) endif enddo !.. Update to ice do cc = 1, cat !.. Riming update if((qi(cc,k)+(prdr(cc))*dt).gt.QSMALL) then if(dry_growth(cc)) then rhobar(cc) = rhobar(cc)*(qi(cc,k)/(qi(cc,k)+(prdr(cc))*dt)) + & rhorimeout(cc)*(1.-(qi(cc,k)/(qi(cc,k)+(prdr(cc))*dt))) else rhobar(cc) = rhobar(cc)*(qi(cc,k)/(qi(cc,k)+(prdr(cc))*dt)) + & RHOW*(1.-(qi(cc,k)/(qi(cc,k)+(prdr(cc))*dt))) rhobar(cc) = max(rhobar(cc),RHOI) endif else rhobar(cc) = RHOI endif qi(cc,k)=qi(cc,k)+(prdr(cc))*dt ai(cc,k)=ai(cc,k)+(ardr(cc))*dt ci(cc,k)=ci(cc,k)+(crdr(cc))*dt enddo !.. Keep everything in bounds after vapor growth and riming do cc = 1, cat if(qi(cc,k).gt.QSMALL) then ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k) ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k) deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) dsold(cc) = deltastr(cc) endif enddo !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Update to ice after melting if(temp.gt.T0) then do cc = 1, cat if(qi(cc,k).gt.QSMALL) then rhobar(cc) = rhobar(cc)*((qi(cc,k)+(qmlt(cc))*dt)/qi(cc,k)) + & RHOW*((-qmlt(cc)*dt)/qi(cc,k)) rhobar(cc) = min(rhobar(cc),RHOI) else rhobar(cc) = RHOI endif qi(cc,k)=qi(cc,k)+qmlt(cc)*dt ni(cc,k)=ni(cc,k)+nmlt(cc)*dt ai(cc,k)=ai(cc,k)+amlt(cc)*dt ci(cc,k)=ci(cc,k)+cmlt(cc)*dt !.. Small masses of leftover ice can become rain if(qi(cc,k).le.QSMALL) then qr(k) = qr(k) + qi(cc,k) qi(cc,k)=0. endif !.. Keep everyting in bounds after melting if(qi(cc,k).gt.QSMALL) then ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k) ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k) deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) !.. Melting should only make the average distribution more spherical if(dsold(cc).le.1.) then if(deltastr(cc).lt.dsold(cc)) then deltastr(cc)=dsold(cc) endif if(deltastr(cc).gt.1.) then deltastr(cc)=1. endif else if(deltastr(cc).gt.dsold(cc)) then deltastr(cc)=dsold(cc) endif if(deltastr(cc).lt.1.) then deltastr(cc)=1. endif endif cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc) ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k) ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k) call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) dsold(cc) = deltastr(cc) endif !.. End if ice exists after melting enddo !.. Loop over all ices species endif !.. End temperature above melting !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Update to vapor qv(k)=qv(k)+(-pre-mnuccd)*dt do cc = 1, cat qv(k)=qv(k)+(-prd(cc))*dt enddo !.. Update to cloud qc(k)=qc(k)+(-prc-pra-mim)*dt do cc = 1, cat qc(k) = qc(k)-(prdr(cc)*qcrimefrac(cc)*dt) enddo !.. Update to rain qr(k)=qr(k)+(pre+prc+pra-mimr-mbiggr)*dt nr(k)=nr(k)+(nprc1+nragg+npre-nimr-nbiggr)*dt do cc = 1, cat qr(k)=qr(k)+((-prdr(cc)*(1.-qcrimefrac(cc)))-qmlt(cc)-dQRfzri(cc))*dt nr(k)=max((nr(k)+((-qi_qr_nrn(cc)*ni(cc,k)*nr(k)*rhoair(k))-nmlt(cc)-dNfzri(cc))*dt),0.) enddo !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Store the number before nucleation for later use do cc = 1, cat nibnuc(cc)=ni(cc,k) enddo if(igr.le.1.) then !.. Nucleation to ice1 (planar-nuclated) totalnuc(ICE1) =(mnuccd+mim+mimr+mbiggr+dQIfzri(ICE2)+dQIfzri(ICE3))*dt totalnucn(ICE1)=(nnuccd+nim+nimr+nbiggr+dNfzri(ICE2)+dNfzri(ICE3))*dt do cc = 1, cat totalnuc(ICE1) =totalnuc(ICE1) + (qmult(cc)+dQRfzri(cc))*dt totalnucn(ICE1)=totalnucn(ICE1) + (nmult(cc))*dt enddo qi(ICE1,k)=qi(ICE1,k)+totalnuc(ICE1) ni(ICE1,k)=ni(ICE1,k)+totalnucn(ICE1) !.. Loss of ice2 from ice2-rain producing ice1 qi(ICE2,k)=qi(ICE2,k)-(dQIfzri(ICE2))*dt ni(ICE2,k)=ni(ICE2,k)-(dNfzri(ICE2))*dt !.. Loss of ice3 from ice3-rain producing ice1 qi(ICE3,k)=qi(ICE3,k)-(dQIfzri(ICE3))*dt ni(ICE3,k)=ni(ICE3,k)-(dNfzri(ICE3))*dt else !.. Nucleation to ice2 (columnar-nucleated) totalnuc(ICE2) =(mnuccd+mim+mimr+mbiggr+dQIfzri(ICE1)+dQIfzri(ICE3))*dt totalnucn(ICE2)=(nnuccd+nim+nimr+nbiggr+dNfzri(ICE1)+dNfzri(ICE3))*dt do cc = 1, cat totalnuc(ICE2) =totalnuc(ICE2) + (qmult(cc)+dQRfzri(cc))*dt totalnucn(ICE2)=totalnucn(ICE2) + (nmult(cc))*dt enddo qi(ICE2,k)=qi(ICE2,k)+totalnuc(ICE2) ni(ICE2,k)=ni(ICE2,k)+totalnucn(ICE2) !.. Loss of ice1 from ice1-rain producing ice2 qi(ICE1,k)=qi(ICE1,k)-(dQIfzri(ICE1))*dt ni(ICE1,k)=ni(ICE1,k)-(dNfzri(ICE1))*dt !.. Loss of ice3 from ice3-rain producing ice2 qi(ICE3,k)=qi(ICE3,k)-(dQIfzri(ICE3))*dt ni(ICE3,k)=ni(ICE3,k)-(dNfzri(ICE3))*dt endif !.. Limit ice number concentrations to 1000/L do cc = 1, cat ni(cc,k)=min(ni(cc,k),(1000.*1000.*i_rhoair(k))) enddo !.. Final check on rain and fall speed calculation if (qr(k).gt.QSMALL) then nr(k) = max(nr(k),QNSMALL) lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333 n0rr = nr(k)*lamr if(lamr.LT.LAMMINR) then lamr = LAMMINR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr elseif(lamr.gt.LAMMAXR) then lamr = LAMMAXR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr endif vtrn(k) = ARN*(gamma(1.+BR))/LAMR**BR vtrm(k) = ARN*(gamma(4.+BR))/6./LAMR**BR vtrn(k) = min(vtrn(k),9.1) vtrm(k) = min(vtrm(k),9.1) dbzr = rhoair(k)*nr(k)/lamr**3/lamr**3*720. dbzr = 1.e18*dbzr else !.. If only small mass of rain, push to vapor qv(k)=qv(k)+qr(k) theta(k)=theta(k)+theta(k)*i_temp*(qr(k)*xxlv)*i_cp qr(k) =0. nr(k) =0. vtrm(k) =0. vtrn(k) =0. endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! do cc = 1, cat !.. Check for ice after nucleation and then update density and bulk distribution sizes if(qi(cc,k).ge.QSMALL) then ni(cc,k) =max(ni(cc,k),QNSMALL) nucfrac = totalnuc(cc)/qi(cc,k) nucfrac = min(nucfrac, 1.0) nucfrac = max(nucfrac, 0.0) !.. Do not nucleate to aggregates if(cc.eq.ICE3) then nucfrac=0. endif !.. Update density from nucleation rhobar(cc) = RHOI*nucfrac + (1.-nucfrac)*rhobar(cc) !.. Keep everything in bounds ai(cc,k) =max(ai(cc,k),QASMALL) ci(cc,k) =max(ci(cc,k),QASMALL) nibnuc(cc)=max(nibnuc(cc),QNSMALL) ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*nibnuc(cc)))**0.333333333333),2.e-6) cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*nibnuc(cc)))**0.333333333333),2.e-6) ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k) ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k) aniold(cc) = ani(cc) cniold(cc) =ao**(1.-dsold(cc))*aniold(cc)**dsold(cc) !.. Get the size distribution properties of nucleated particles if(nucfrac.gt.1.e-5.and.totalnucn(cc).gt.0.) then gamma_arg = NU+3. gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+4. gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) anuc=(((totalnuc(cc)*gammnu)/(RHOI*totalnucn(cc)*fourthirdspi* & gamma_tab(gi)))**0.3333333)*gamma_tab(gi2)/gamma_tab(gi) anuc = MAX(anuc,2.e-6) anuc = MIN(anuc,3.e-3) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+3.+deltastr(cc) gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+2.+(2.*deltastr(cc)) gi3=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) amass=ani(cc)*gamma_tab(gi2)/gamma_tab(gi) cmass=cni(cc)*gamma_tab(gi3)/gamma_tab(gi) !.. Mass-weighted size of old distribution and nucleated distribution amass = amass*(1.-nucfrac) + anuc*nucfrac cmass = cmass*(1.-nucfrac) + anuc*nucfrac !.. New ani and cni after nucleation ani(cc) = max((amass*gamma_tab(gi)/gamma_tab(gi2)),2.e-6) cni(cc) = max((cmass*gamma_tab(gi)/gamma_tab(gi3)),2.e-6) !.. New delta* (shape) after nucleation if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) else deltastr(cc)=1. endif !.. Nucleation should only make the average distribution more spherical if(dsold(cc).le.1.) then if(deltastr(cc).lt.dsold(cc)) then deltastr(cc)=dsold(cc) endif if(deltastr(cc).gt.1.) then deltastr(cc)=1. endif else if(deltastr(cc).gt.dsold(cc)) then deltastr(cc)=dsold(cc) endif if(deltastr(cc).lt.1.) then deltastr(cc)=1. endif endif endif !.. End nucleation update distribution properties !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. This needs to be done because of transfer from e.g. !.. ice1 to ice2 during nucleation. !.. Assume constant density/delta*, loss of mass and number !..With the new delta* and density after nucleation re-diagnose a and c !..Use the old density and deltastr for loss values e.g. q1 = q1 - (transfered to q2) alphstr=ao**(1.-deltastr(cc)) alphv=fourthirdspi*alphstr betam=2.+deltastr(cc) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) ani(cc)=((qi(cc,k)*gammnu)/(rhobar(cc)*ni(cc,k)*alphv* & gamma_tab(gi)))**(1./betam) cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc) !.. volume and volume*aspect ratio process rates ai(cc,k) = ai(cc,k) + cniold(cc)*nibnuc(cc)*2.*aniold(cc)*(ani(cc)-aniold(cc)) + & nibnuc(cc)*aniold(cc)**2*(cni(cc)-cniold(cc)) + & aniold(cc)**2*cniold(cc)*(ni(cc,k)-nibnuc(cc)) ci(cc,k) = ci(cc,k) + aniold(cc)*nibnuc(cc)*2.*cniold(cc)*(cni(cc)-cniold(cc)) + & nibnuc(cc)*cniold(cc)**2*(ani(cc)-aniold(cc)) + & cniold(cc)**2*aniold(cc)*(ni(cc,k)-nibnuc(cc)) !.. Keep everything in bounds ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) else deltastr(cc)=1. endif call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) dsold(cc) = deltastr(cc) !.. Limit the fall speeds here to 25 m s^-1 vtrzi1(cc,k)=min(vtrzi1(cc,k),25.) vtrmi1(cc,k)=min(vtrmi1(cc,k),25.) vtrni1(cc,k)=min(vtrni1(cc,k),25.) else !.. If there is no ice even after nucleation, !.. add remainder to vapor and update temperature qv(k)=qv(k)+qi(cc,k) theta(k)=theta(k)+theta(k)*i_temp*(qi(cc,k)*xxls)*i_cp qi(cc,k)=0. ai(cc,k)=0. ni(cc,k)=0. ci(cc,k)=0. vtrzi1(cc,k)=0. vtrmi1(cc,k)=0. vtrni1(cc,k)=0. endif !.. End check for ice enddo !.. Loop over all ice species !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! if(temp.le.T0) then !.. Aggregation occurs similarly to nulceation !.. We assume that the loss of ice1 and ice2 to aggregates does !.. not significantly change the shape and density of the !.. remaining particles. Both ani and cni are updated !.. with no change to density and delta*. !.. Aggregation (see Jensen et al. 2017, 2018a, 2018b) do cc = 1, cat dsold(cc) = deltastr(cc) nibnuc(cc) = max(ni(cc,k),QNSMALL) qi(cc,k) = qi(cc,k) + qagg(cc) ni(cc,k) = ni(cc,k) + nagg(cc) if(qi(cc,k).gt.QSMALL) then !.. Aggregate (unmelted) density and aspect ratio are assumed !.. Force the shape of aggregates if(cc.eq.ICE3) then rhobar(cc)=50. phiagg3 = 0.2 ani(cc)=.5*dn3 gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) cni(cc)=phiagg3*ani(cc)*gammnu/gamma_tab(gi) ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k) ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k) ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) else deltastr(cc)=1. endif call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) else !.. Variable check on ice-1 and ice-2 after aggregation ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*nibnuc(cc)))**0.333333333333),2.e-6) cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*nibnuc(cc)))**0.333333333333),2.e-6) aniold(cc) = ani(cc) cniold(cc) = cni(cc) deltastr(cc)= dsold(cc) !.. Update ice1 and ice2 sizes after aggregation from !.. updated mass and number mixiing ratios (transfer to aggregates) and an !.. assumed constant deltastar alphstr=ao**(1.-deltastr(cc)) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) ani(cc) = ((qi(cc,k))/(ni(cc,k)*rhobar(cc)*fourthirdspi*alphstr*gamma_tab(gi)*i_gammnu))**& (1./(2.+deltastr(cc))) ani(cc)=max(ani(cc),2.e-6) cni(cc)=alphstr*ani(cc)**deltastr(cc) !.. volume and volume*aspect ratio updates ai(cc,k) = ai(cc,k) + cniold(cc)*nibnuc(cc)*2.*aniold(cc)*(ani(cc)-aniold(cc)) + & nibnuc(cc)*aniold(cc)**2*(cni(cc)-cniold(cc)) + & aniold(cc)**2*cniold(cc)*(ni(cc,k)-nibnuc(cc)) ci(cc,k) = ci(cc,k) + aniold(cc)*nibnuc(cc)*2.*cniold(cc)*(cni(cc)-cniold(cc)) + & nibnuc(cc)*cniold(cc)**2*(ani(cc)-aniold(cc)) + & cniold(cc)**2*aniold(cc)*(ni(cc,k)-nibnuc(cc)) !.. Keep everything in bounds ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) else deltastr(cc)=1. endif call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) endif endif !.. If ice exist enddo !.. Loop over all species endif !.. If temperature < 273.15 K !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Update the thermodynamics theta(k)=theta(k)+theta(k)*i_temp*((mnuccd)*xxls*i_cp*dt) do cc = 1, cat theta(k)=theta(k)+theta(k)*i_temp*((prd(cc))*xxls*i_cp*dt) theta(k)=theta(k)+theta(k)*i_temp*((qmlt(cc))*xxlf*i_cp*dt) enddo theta(k)=theta(k)+theta(k)*i_temp*(pre*xxlv)*i_cp*dt if(temp.le.T0) then theta(k)=theta(k)+theta(k)*i_temp*(mim+mimr+mbiggr)*xxlf*i_cp*dt do cc = 1, cat theta(k)=theta(k)+theta(k)*i_temp*(prdr(cc)+dQRfzri(cc))*xxlf*i_cp*dt enddo endif temp = theta(k)/(100000./pres_e(k))**(RCP) i_temp= 1./temp qvs = 0.622*polysvp(temp,0)/(pres_e(k)-polysvp(temp,0)) xxlv = 3.1484E6-2370.*temp rhoair(k)= pres_e(k)/(RD*temp) i_rhoair(k) = 1./rhoair(k) !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Saturation adjustment for qc pcc = (qv(k)-qvs) / (1.0 + xxlv**2*qvs/(CP*rv*temp**2))*i_dt if (pcc*dt+qc(k).lt.0.) then pcc=-qc(k)/dt end if theta(k)=theta(k) + theta(k)*i_temp*pcc*xxlv*i_cp*dt qv(k) = qv(k) - pcc*dt qc(k) = qc(k) + pcc*dt if(qc(k).le.QSMALL) then qv(k)=qv(k)+qc(k) theta(k)=theta(k)+theta(k)*i_temp*(qc(k)*xxlv)*i_cp qc(k)=0. nc(k)=0. endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Cloud water effective radius lrsig = 1.3 sig = log(lrsig) !.. Distribution standard deviation lwc = qc(k)*rhoair(k) !.. liquid water content ncm3dum = nc(k)*rhoair(k) if(ncm3dum .gt. 0.) then r_n = ((lwc)/(fourthirdspi*ncm3dum*RHOW*exp(4.5*(sig**2))))**(0.33333333333) else r_n = 0. endif effc1d(k) = r_n * (5.*EXP((sig**2)/2.)) effc1d(k) = MIN(effc1d(k),50.e-6) effc1d(k) = MAX(effc1d(k),2.51e-6) ! For now, assume cloud droplet fallspeed is zero vtrmc(k) = 0.0 vtrnc(k) = 0.0 !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Final check on ice before sedimentation do cc = 1, cat if(qi(cc,k).gt.QSMALL) then ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k) ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k) deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) !.. Final check on aggregates if(cc.eq.ICE3) then if(temp.le.T0) then rhobar(cc)=50. phiagg3 = 0.2 gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) cni(cc)=phiagg3*ani(cc)*gammnu/gamma_tab(gi) ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k) ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k) endif if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) else deltastr(cc)=1. endif !.. Keep aggregates spherical or oblate (for now) deltastr(cc) = min(deltastr(cc),1.) alphstr=ao**(1.-deltastr(cc)) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) ani(cc) = ((qi(cc,k))/(ni(cc,k)*rhobar(cc)*fourthirdspi*alphstr*& gamma_tab(gi)/gammnu))**(1./(2.+deltastr(cc))) ani(cc) = max(ani(cc),2.e-6) cni(cc)=alphstr*ani(cc)**deltastr(cc) !.. Final check on aggregate size (implicit breakup) if(ani(cc).gt.0.5e-3) then ani(cc)=0.5e-3 cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) ni(cc,k) = (qi(cc,k))/(rhobar(cc)*fourthirdspi*ao**(1.-deltastr(cc))* & ani(cc)**(2.+deltastr(cc))*gamma_tab(gi)/gammnu) ni(cc,k)=max(ni(cc,k),QNSMALL) endif ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k) ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k) call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) endif ! Check on aggregates !.. Fall speed moved to here !.. See Harrington et al. 2013 and Mitchell and Heymsfield 2005 fv(cc) = 1.0 fh(cc) = 1.0 gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) phiivt = cni(cc)/ani(cc)*gamma_tab(gi)*i_gammnu alphstr=ao**(1.-deltastr(cc)) !.. Determine coefficeients for the Best number for fall speeds if(phiivt.lt.1.0)then bl = 1. al = 2. aa = PI ba = 2. qe = (1.-phiivt)*(rhobar(cc)/RHOI) + phiivt else if(phiivt .gt. 1.0)then al = 2.0 bl = 1. aa = PI*alphstr ba = deltastr(cc) + 1. qe = 1.0 else if(phiivt .eq. 1.0)then bl = 1. al = 2. aa = PI ba = 2. qe = 1.0 endif qe = min(qe, 1.0) !.. Fall speed and ventilation (Best number formulation) xn = 2./rhobar(cc)*(rhobar(cc)-rhoair(k))*G_HOME*rhoair(k)/mu**2 * & (fourthirdspi*rhobar(cc))*alphstr*al**2/(aa*qe) * qe**(3./4.) bx = deltastr(cc)+2.+2.*bl-ba !.. Number-average Best Number xm = xn*ani(cc)**bx * (gamma(NU+bx))*i_gammnu !.. The following fall speed coefficients are from !.. Mitchell and Heymsfiled 2005 f_c1 = 4.0 / (5.83 * 5.83 * SQRT(0.6)) f_c2 = (5.83 * 5.83) / 4.0 bm = ((f_c1 * SQRT(xm)) / & (2.0 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0) * & (SQRT(1.0 + f_c1 * SQRT(xm))))) - & (1.0e-5 * xm) / & (f_c2 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0)**2.0) am = ((f_c2 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0)**2.0) - & (1.0e-5 * xm)) / (xm**bm) if(xm.gt.1.e8)then am = 1.0865 bm = 0.499 endif !.. Reynolds Number Nre = am*xm**bm !.. Number-averaged ice fall speed vtrni1(cc,k) = min((mu/rhoair(k)*0.5 * am*(xn)**bm*ani(cc)**(bx*bm-1.) * & (gamma(NU+bx*bm-1.))*i_gammnu),25.) !.. Mass-averaged ice fall speed vtrmi1(cc,k) = min((mu/rhoair(k)*0.5 * am*(xn)**bm*ani(cc)**(bx*bm-1.) * & (gamma(NU+bx*bm-1.+2.+deltastr(cc)))/(gamma(NU+2.+deltastr(cc)))),25.) !.. Remove size sorting when melting if(temp.gt.T0.and.qmlt(cc).lt.0.) then vtrni1(cc,k) = vtrmi1(cc,k) endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Determine the ice radiative properties !.. For simlicity this is !.. Planar ice: integral of a^2c / integral of a^2 !.. Columnar ice: integral of a^2c / integral of a*c gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) phirad = cni(cc)/ani(cc)*gamma_tab(gi)*i_gammnu if(phirad.le.1.) then gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+2. gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) effi(cc,k) = cni(cc)*gamma_tab(gi)/gamma_tab(gi2) else gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+1.+deltastr(cc) gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) effi(cc,k) = ani(cc)*gamma_tab(gi)/gamma_tab(gi2) endif !.. For advection use volume and volume times aspect ratio mixing ratios ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k) ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) ni(cc,k)=max(ni(cc,k),QNSMALL) else !.. If not enough ice convert to vapor qv(k)=qv(k)+qi(cc,k) theta(k)=theta(k)+theta(k)*i_temp*(qi(cc,k)*xxls)*i_cp qi(cc,k)=0. ai(cc,k)=0. ni(cc,k)=0. ci(cc,k)=0. vtrzi1(cc,k)=0. vtrmi1(cc,k)=0. vtrni1(cc,k)=0. endif !.. If ice exists !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Output diagnostics gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) !.. For aspect ratio diagnostic set lower size limit if(ai(cc,k).gt.1.e-12.and.ci(cc,k).gt.1.e-12.and. & ani(cc).gt.2.e-6.and.cni(cc).gt.2.e-6.and.qi(cc,k).gt.1.e-9) then phii1d(cc,k) = cni(cc)/ani(cc)*gamma_tab(gi)*i_gammnu rhopo1d(cc,k) = rhobar(cc) endif gamma_arg = NU+3.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+2.+deltastr(cc) gi2=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gamma_arg = NU+2.+(2.*deltastr(cc)) gi3=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) if(ai(cc,k).gt.1.e-12.and.ci(cc,k).gt.1.e-12.and. & ani(cc).gt.2.e-6.and.cni(cc).gt.2.e-6.and.qi(cc,k).gt.1.e-9) then di1d(cc,k) = MAX((2.*ani(cc)*gamma_tab(gi)/gamma_tab(gi2)), & (2.*cni(cc)*gamma_tab(gi3)/gamma_tab(gi2))) vmi1d(cc,k) = vtrmi1(cc,k) endif if(ni(cc,k).gt.0..and.ani(cc).gt.0.) then vtrzi1(cc,k) = 0.176/0.93*(6./PI)**2*((fourthirdspi*rhobar(cc)* & (ao)**(1.-deltastr(cc)))/(2.**(2.+deltastr(cc))))**2 * & rhoair(k)/900./900.*ni(cc,k)*(2.*ani(cc))**(2.*(2.+deltastr(cc))) * & gamma(NU+2.*(2.+deltastr(cc)))*i_gammnu vtrzi1(cc,k) = 1.e18*vtrzi1(cc,k) else vtrzi1(cc,k) = 0. endif enddo !.. Loop over all ice species dbzsum = dbzr + vtrzi1(ICE1,k) + vtrzi1(ICE2,k) + vtrzi1(ICE3,k) if(dbzsum.ne.0.) then dbz1d(k) = 10.*Log10(dbzsum) dbz1d(k) = MAX(-35.,dbz1d(k)) endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Ice effective radius is a mass-weighting of all ice species totaliwc = 0. do cc = 1, cat totaliwc = totaliwc + qi(cc,k) enddo if(totaliwc.gt.0.) then radw1 = MAX((MIN((qi(ICE1,k)/totaliwc),1.)),0.) radw2 = MAX((MIN((qi(ICE2,k)/totaliwc),1.)),0.) radw3 = MAX((MIN((qi(ICE3,k)/totaliwc),1.)),0.) effi1d(k) = MAX((MIN((effi(ICE1,k)*radw1 + effi(ICE2,k)*radw2 + & effi(ICE3,k)*radw3),999.e-6)),5.e-6) else effi1d(k)=0. endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Check for rain or ice to do sedimentation if(qr(k).gt.QSMALL) then sedi=.true. endif do cc = 1, cat if(qi(cc,k).gt.QSMALL) then sedi=.true. endif enddo !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. else on domicro else !.. If we don't do all of the above !.. Push to vapor and get density for seimentation in the column xxls = 3.15e6-2370.*temp+0.3337e6 xxlv = 3.1484E6-2370.*temp i_temp= 1./temp rhoair(k)= pres_e(k)/(RD*temp) i_rhoair(k) = 1./rhoair(k) qv(k)=qv(k)+qc(k)+qr(k) theta(k)=theta(k)+theta(k)*i_temp*((qc(k)+qr(k))*xxlv)*i_cp do cc = 1, cat qv(k)=qv(k)+qi(cc,k) theta(k)=theta(k)+theta(k)*i_temp*((qi(cc,k))*xxls)*i_cp qi(cc,k)=0. ni(cc,k)=0. ai(cc,k)=0. ci(cc,k)=0. vtrmi1(cc,k)=0. vtrni1(cc,k)=0. vtrzi1(cc,k)=0. vmi1d(cc,k)=0. di1d(cc,k)=0. phii1d(cc,k)=0. rhopo1d(cc,k)=0. qipre1d(cc) =0. icetype1d(cc,k)=0. enddo effc1d(k)=0. effi1d(k)=0. dbz1d(k)=-35. qc(k)=0. qr(k)=0. nr(k)=0. vtrm(k)=0. vtrn(k)=0. qrpre1d=0. endif !.. If domicro enddo !.. do k = kts, kte !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Sedimentation (Reisner et al. 1998) !.. Fallout terms are calculated on split time steps to ensure numerical stability (Courant < 1) nstep = 1 maxfall=0. if(sedi) then do k = kte,kts,-1 !.. Calculate number of split timesteps maxfall = MAX(vtrmi1(ICE1,k),vtrni1(ICE1,k),vtrmi1(ICE2,k),vtrni1(ICE2,k), & vtrmi1(ICE3,k),vtrni1(ICE3,k),vtrm(k),vtrn(k)) nstep = MAX(INT(maxfall*dt/dzmic(k)+1.),nstep) do cc = 1, cat qi(cc,k) = qi(cc,k)*rhoair(k) ni(cc,k) = ni(cc,k)*rhoair(k) ai(cc,k) = ai(cc,k)*rhoair(k) ci(cc,k) = ci(cc,k)*rhoair(k) enddo qr(k) = qr(k)*rhoair(k) nr(k) = nr(k)*rhoair(k) enddo do nn = 1, nstep do k = kts, kte do cc = 1, cat fluxqi(cc,k) = vtrmi1(cc,k)*qi(cc,k) fluxni(cc,k) = vtrni1(cc,k)*ni(cc,k) fluxai(cc,k) = vtrmi1(cc,k)*ai(cc,k) fluxci(cc,k) = vtrmi1(cc,k)*ci(cc,k) enddo fluxqr(k) = vtrm(k)*qr(k) fluxnr(k) = vtrn(k)*nr(k) enddo do cc = 1, cat falltndqi(cc) = fluxqi(cc,KTE)/dzmic(KTE) falltndni(cc) = fluxni(cc,KTE)/dzmic(KTE) falltndai(cc) = fluxai(cc,KTE)/dzmic(KTE) falltndci(cc) = fluxci(cc,KTE)/dzmic(KTE) enddo falltndqr = fluxqr(KTE)/dzmic(KTE) falltndnr = fluxnr(KTE)/dzmic(KTE) do cc = 1, cat qi(cc,KTE) = qi(cc,KTE)-falltndqi(cc)*dt/nstep ni(cc,KTE) = ni(cc,KTE)-falltndni(cc)*dt/nstep ai(cc,KTE) = ai(cc,KTE)-falltndai(cc)*dt/nstep ci(cc,KTE) = ci(cc,KTE)-falltndci(cc)*dt/nstep enddo qr(KTE) = qr(KTE)-falltndqr*dt/nstep nr(KTE) = nr(KTE)-falltndnr*dt/nstep do k = kte-1,kts,-1 do cc = 1, cat falltndqi(cc) = (fluxqi(cc,k+1)-fluxqi(cc,k))/dzmic(k) falltndni(cc) = (fluxni(cc,k+1)-fluxni(cc,k))/dzmic(k) falltndai(cc) = (fluxai(cc,k+1)-fluxai(cc,k))/dzmic(k) falltndci(cc) = (fluxci(cc,k+1)-fluxci(cc,k))/dzmic(k) enddo falltndqr = (fluxqr(k+1)-fluxqr(k))/dzmic(k) falltndnr = (fluxnr(k+1)-fluxnr(k))/dzmic(k) do cc = 1, cat qi(cc,k) = qi(cc,k)+falltndqi(cc)*dt/nstep ni(cc,k) = ni(cc,k)+falltndni(cc)*dt/nstep ai(cc,k) = ai(cc,k)+falltndai(cc)*dt/nstep ci(cc,k) = ci(cc,k)+falltndci(cc)*dt/nstep !.. GET PRECIPITATION AND SNOWFALL ACCUMULATION DURING THE TIME STEP !.. FACTOR OF 1000 CONVERTS FROM M TO MM, BUT DIVISION BY DENSITY !.. OF LIQUID WATER CANCELS THIS FACTOR OF 1000 !.. Ice precipitation rate if(k.eq.kts.and.qi(cc,k).gt.1.e-9) then qipre1d(cc) = qipre1d(cc)+fluxqi(cc,k)*dt/nstep endif enddo qr(k) = qr(k)+falltndqr*dt/nstep nr(k) = nr(k)+falltndnr*dt/nstep !.. Rain precipitation rate if(k.eq.kts.and.qr(k).gt.1.e-9) then qrpre1d = qrpre1d + fluxqr(k)*dt/nstep endif enddo enddo do k = kts, kte do cc = 1, cat qi(cc,k) = qi(cc,k)*i_rhoair(k) ni(cc,k) = ni(cc,k)*i_rhoair(k) ai(cc,k) = ai(cc,k)*i_rhoair(k) ci(cc,k) = ci(cc,k)*i_rhoair(k) enddo qr(k) = qr(k)*i_rhoair(k) nr(k) = nr(k)*i_rhoair(k) enddo endif !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Check after sedimentation do k = kts, kte temp = theta(k)/(100000./pres_e(k))**(RCP) i_temp= 1./temp xxls = 3.15e6-2370.*temp+0.3337e6 xxlv = 3.1484E6-2370.*temp do cc = 1, CAT if(qi(cc,k).gt.QSMALL) then !.. Keep variables in bounds ni(cc,k)=max(ni(cc,k),QNSMALL) ai(cc,k)=max(ai(cc,k),QASMALL) ci(cc,k)=max(ci(cc,k),QASMALL) ani(cc) =max((((ai(cc,k)**2)/(ci(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) cni(cc) =max((((ci(cc,k)**2)/(ai(cc,k)*ni(cc,k)))**0.333333333333),2.e-6) !.. Set incoming limit after sedimentation !.. smallest values for which we care about shape if(ai(cc,k).lt.1.e-12.or.ci(cc,k).lt.1.e-12) then ai(cc,k)=min(ai(cc,k),ci(cc,k)) ci(cc,k)=ai(cc,k) ani(cc)=(ai(cc,k)/ni(cc,k))**0.3333333333 cni(cc)=(ci(cc,k)/ni(cc,k))**0.3333333333 endif if(ani(cc).lt.2.e-6.or.cni(cc).lt.2.e-6) then ani(cc)=2.e-6 cni(cc)=2.e-6 endif ci(cc,k) = cni(cc)**2*ani(cc)*ni(cc,k) ai(cc,k) = ani(cc)**2*cni(cc)*ni(cc,k) deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) !.. Final check on aggregates if(cc.eq.ICE3) then if(temp.le.T0) then rhobar(cc)=50. phiagg3 = 0.2 gamma_arg = NU-1.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) cni(cc)=phiagg3*ani(cc)*gammnu/gamma_tab(gi) endif ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k) ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k) if(ani(cc).gt.(1.1*ao).and.cni(cc).gt.(1.1*ao)) then deltastr(cc)=(log(cni(cc))-log(ao))/(log(ani(cc))-log(ao)) else deltastr(cc)=1. endif !.. Keep aggregates spherical or oblate (for now) deltastr(cc) = min(deltastr(cc),1.) alphstr=ao**(1.-deltastr(cc)) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) ani(cc) = ((qi(cc,k))/(ni(cc,k)*rhobar(cc)*fourthirdspi*alphstr*& gamma_tab(gi)/gammnu))**(1./(2.+deltastr(cc))) ani(cc) = max(ani(cc),2.e-6) cni(cc)=alphstr*ani(cc)**deltastr(cc) !.. Final check on aggregate size (implicit breakup) if(ani(cc).gt.0.5e-3) then ani(cc)=0.5e-3 cni(cc)=ao**(1.-deltastr(cc))*ani(cc)**deltastr(cc) gamma_arg = NU+2.+deltastr(cc) gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) ni(cc,k) = (qi(cc,k))/(rhobar(cc)*fourthirdspi*ao**(1.-deltastr(cc))* & ani(cc)**(2.+deltastr(cc))*gamma_tab(gi)/gammnu) ni(cc,k)=max(ni(cc,k),QNSMALL) endif ai(cc,k)=ani(cc)**2*cni(cc)*ni(cc,k) ci(cc,k)=cni(cc)**2*ani(cc)*ni(cc,k) call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) endif !.. Check on aggregates call var_check(NU, ao, fourthirdspi, gammnu, qi(cc,k), deltastr(cc), & ani(cc), cni(cc), rni(cc), rhobar(cc), ni(cc,k), ai(cc,k), ci(cc,k), & alphstr, alphv, betam) else qv(k)=qv(k)+qi(cc,k) theta(k)=theta(k)+theta(k)*i_temp*(qi(cc,k)*xxls)*i_cp qi(cc,k)=0. ni(cc,k)=0. ai(cc,k)=0. ci(cc,k)=0. endif enddo if(qr(k).gt.QSMALL) then nr(k) = max(nr(k),QNSMALL) lamr = (PI*RHOW*nr(k)/qr(k))**0.333333333 n0rr = nr(k)*lamr if(lamr.LT.LAMMINR) then lamr = LAMMINR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr elseif(lamr.gt.LAMMAXR) then lamr = LAMMAXR n0rr = lamr**4*qr(k)/(PI*RHOW) nr(k) = n0rr/lamr endif else qv(k)=qv(k)+qr(k) theta(k)=theta(k)+theta(k)*i_temp*(qr(k)*xxlv)*i_cp qr(k) =0. nr(k) =0. endif enddo end subroutine me_ishmael !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. This subroutine keeps the ice properties in check subroutine var_check(NU, ao, fourthirdspi, gammnu, qidum, dsdum, ani, & cni, rni, rbdum, nidum, aidum, cidum, alphstr, alphv, betam) implicit none REAL, INTENT(IN) :: NU, ao, fourthirdspi, gammnu, qidum REAL :: voltmp, maxsize, gamma_arg INTEGER :: gi REAL, INTENT(INOUT) :: dsdum, ani, cni, rbdum, nidum, aidum, cidum REAL, INTENT(OUT) :: rni, alphstr, alphv, betam !.. Deltastr check if(dsdum.lt.0.55) then dsdum=0.55 ani=(cni/(ao**(1.-dsdum)))**(1./dsdum) aidum=ani**2*cni*nidum cidum=cni**2*ani*nidum else if (dsdum.gt.1.3) then dsdum=1.3 cni=ao**(1.-dsdum)*ani**dsdum aidum=ani**2*cni*nidum cidum=cni**2*ani*nidum endif alphstr=ao**(1.-dsdum) alphv = fourthirdspi*alphstr betam = 2.0 + dsdum gamma_arg = NU+2.+dsdum gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) !.. Ice density check !.. Keep ice density between 50 and RHOI=920 kg m^-3 if(ani.gt.2.e-6) then rbdum = qidum*gammnu/(nidum*alphv* & ani**betam*gamma_tab(gi)) else rbdum = RHOI endif if(rbdum.gt.RHOI) then rbdum=RHOI ani=((qidum*gammnu)/(rbdum*nidum*alphv* & gamma_tab(gi)))**(1./betam) cni=ao**(1.-dsdum)*ani**dsdum aidum=ani**2*cni*nidum cidum=cni**2*ani*nidum elseif(rbdum.lt.50.) then rbdum=50. ani=((qidum*gammnu)/(rbdum*nidum*alphv* & gamma_tab(gi)))**(1./betam) cni=ao**(1.-dsdum)*ani**dsdum aidum=ani**2*cni*nidum cidum=cni**2*ani*nidum endif !.. Small ice limit !.. Keep rni > 2 microns rni= (qidum*3./(nidum*rbdum*4.*PI* & (gamma_tab(gi)/gammnu)))**0.333333333333 if(rni.lt.2.e-6) then rni=2.e-6 nidum=3.*qidum*gammnu/(4.*pi*rbdum*rni**3* & (gamma_tab(gi))) ani=((qidum*gammnu)/(rbdum*nidum*alphv* & gamma_tab(gi)))**(1./betam) cni=ao**(1.-dsdum)*ani**dsdum aidum=ani**2*cni*nidum cidum=cni**2*ani*nidum endif !.. Large ice limit !.. This is a number weighted diameter of 8 mm and !.. a spherical mass weighted diameter of 14 mm maxsize = max(ani,cni) if (maxsize.gt.1.e-3) then if(ani.ge.cni) then ani = 1.e-3 nidum=qidum*gammnu/(fourthirdspi*rbdum*ao**(1.-dsdum)*ani**(2.+dsdum)* & (gamma_tab(gi))) cni=ao**(1.-dsdum)*ani**dsdum aidum=ani**2*cni*nidum cidum=cni**2*ani*nidum else cni = 1.e-3 ani = (cni/(ao**(1.-dsdum)))**(1./dsdum) nidum=qidum*gammnu/(fourthirdspi*rbdum*ao**(1.-dsdum)*ani**(2.+dsdum)* & (gamma_tab(gi))) aidum=ani**2*cni*nidum cidum=cni**2*ani*nidum endif rni= (qidum/(nidum*rbdum*fourthirdspi* & (gamma_tab(gi)/gammnu)))**0.333333333333 end if return end subroutine var_check !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. This is the ventilation, fall speed, vapor growth routine !.. See Harrington et al. 2013 subroutine vaporgrow(dt, ani, cni, rni, igr, nidum, temp, rimesum, presdum, & NU, alphstr, sui, sup, qvs, qvi, mu, iwci, rhodum, qidum, dv, kt, ao, & nsch, npr, gammnu, i_gammnu, fourthirdspi, svpi, xxls, xxlv, xxlf, & capgam, vtbarb, vtbarbm, vtbarbz, anf, cnf, rnf, iwcf, fvdum, fhdum, & rbdum, dsdum, rdout, dsdumout) implicit none real, parameter :: QASMALL= 1.e-19 ! Smallest ice size for sublimation (squared) REAL, INTENT(IN) :: dt, ani, cni, rni, igr, nidum, temp, presdum REAL, INTENT(IN) :: NU, alphstr, mu, iwci, rhodum, qidum REAL, INTENT(IN) :: dv, kt, ao, nsch, npr, gammnu, i_gammnu, fourthirdspi, svpi REAL, INTENT(IN) :: xxls, xxlv, xxlf, capgam, rimesum, sui, sup, qvs, qvi, dsdum REAL :: gammnubet, phii, fs, alphanr, bl, al, aa, ba REAL :: xn, bx, xm, am, bm, f_c1, f_c2, del1, del2, alpha, afn, del REAL :: xvent, ntherm, bv1, bv2, gv, bt1, bt2, gt, vtbranch, rhodep REAL :: videp, vmin, betavol, vi, vf, cf, maxsui, qe, Nre, phif, gamma_arg, rbdumtmp INTEGER :: gi, gi2 REAL, INTENT(INOUT) :: vtbarb, vtbarbm, vtbarbz REAL, INTENT(OUT) :: anf, cnf, rnf, iwcf, rdout REAL, INTENT(INOUT) :: fvdum, fhdum, dsdumout, rbdum fvdum = 1.0 fhdum = 1.0 dsdumout = dsdum gamma_arg = NU-1.+dsdum gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) phii = cni/ani*gamma_tab(gi)*i_gammnu fs = capgam/rni alphanr = ani/rni**(3./(2.+igr)) !.. Determine coefficeients for the Best number for fall speeds if(phii.lt.1.0)then bl = 1. al = 2. aa = PI ba = 2. qe = (1.-phii)*(rbdum/RHOI) + phii else if(phii .gt. 1.0)then al = 2.0 bl = 1. aa = PI*alphstr ba = dsdum + 1. qe = 1.0 else if(phii .eq. 1.0)then bl = 1. al = 2. aa = pi ba = 2. qe = 1.0 endif qe = min(qe, 1.0) !.. Fall speed and ventilation (Best number formulation) xn = 2./rbdum*(rbdum-rhodum)*G_HOME*rhodum/mu**2 * & (fourthirdspi*rbdum)*alphstr*al**2/(aa*qe) * qe**(3./4.) bx = dsdum+2.+2.*bl-ba !.. Number-average Best Number xm = xn*ani**bx * (gamma(nu+bx))*i_gammnu !.. The following fall speed coefficients are from !.. Mitchell and Heymsfiled 2005 f_c1 = 4.0 / (5.83 * 5.83 * SQRT(0.6)) f_c2 = (5.83 * 5.83) / 4.0 bm = ((f_c1 * SQRT(xm)) / & (2.0 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0) * & (SQRT(1.0 + f_c1 * SQRT(xm))))) - & (1.0e-5 * xm) / & (f_c2 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0)**2.0) am = ((f_c2 * (SQRT(1.0 + f_c1 * SQRT(xm)) - 1.0)**2.0) - & (1.0e-5 * xm)) / (xm**bm) if(xm.gt.1.e8)then am = 1.0865 bm = 0.499 endif !.. Reynolds Number Nre = am*xm**bm !.. Number-averaged ice fall speed vtbarb = mu/rhodum*0.5 * am*(xn)**bm*ani**(bx*bm-1.) * & (gamma(nu+bx*bm-1.))*i_gammnu !.. Mass-averaged ice fall speed vtbarbm = mu/rhodum*0.5 * am*(xn)**bm*ani**(bx*bm-1.) * & (gamma(nu+bx*bm-1.+2.+dsdum))/(gamma(nu+2.+dsdum)) !.. Reflectivity-weighted fall speed (if needed) vtbarbz = mu/rhodum*0.5 * am*(xn)**bm*ani**(bx*bm-1.) * & exp(gammln(nu+bx*bm-1.+4.+2.*dsdum))/exp(gammln(nu+4.+2.*dsdum)) !.. Calculate Ventilation xvent = nsch**0.333333333*Nre**0.5 ntherm = Nre**0.5*npr**0.333333333 if(xvent.le.1.0)then bv1 = 1.0 bv2 = 0.14 gv = 2. else bv1 = 0.86 bv2 = 0.28 gv = 1. endif if(ntherm.lt.1.4)then bt1 = 1.0 bt2 = 0.108 gt = 2.0 else bt1 = 0.78 bt2 = 0.308 gt = 1.0 endif fvdum = bv1 + bv2*xvent**gv fhdum = bt1 + bt2*ntherm**gt if(temp.le.T0) then !.. If T < T0 do vapor growth/sublimation (otherwise assume water on surface) !.. Fall speed needed to determine when branching occurs vtbranch = vtbarbm if(sup.ge.0.) then maxsui = 1. elseif(sui.ge.0..and.qvi.lt.qvs) then maxsui = ((sui+1.)*qvi)/(qvs-qvi) - qvi/(qvs-qvi) maxsui = min(maxsui,1.) maxsui = max(maxsui,0.) else maxsui = 0. endif !.. Vapor growth density if(igr.le.1.) then !.. Planar if(vtbranch.gt.0.) then if(ani.gt.SQRT((dv*PI*2.*cni)/(vtbranch*nu))) then rhodep = (RHOI*igr)*maxsui + RHOI*(1.-maxsui) else rhodep = RHOI endif else rhodep = RHOI endif else !.. Columnar rhodep = (RHOI/igr)*maxsui + RHOI*(1.-maxsui) endif !.. high limit on rhodep for vapor growth = 700 kg m^-3 !.. Cotton et al. 2012 (QJRMS) rhodep = min(rhodep,700.) !.. Vapor growth solution which includes heating from rime mass !.. See J.P. Chen's thesis alpha = (dv*fvdum*svpi*xxls)/(RV*kt*fhdum*temp) if(nidum.gt.0.0.and.rimesum.gt.0.0) then del1 = (xxlf*(rimesum/nidum)/(4.*PI*kt*fhdum*capgam)) * & ((temp + alpha*((xxls/(RV*temp))-1.))**(-1.0)) else del1 = 0.0 endif del2 = sui * & (((temp/alpha) + ((xxls/(RV*temp))-1.))**(-1.0)) del = del1 + del2 afn = ((dv*fvdum*polysvp(temp,1))/(RV*temp)) * & (sui - del*((xxls/(RV*temp))-1.)) !.. During sublimation using polynomial removal of density if(afn.lt.0.0) then rhodep = rbdum videp = rni**3 vmin = (10.e-6)**3 if(Vmin.lt.videp)then betavol = log(RHOI/rbdum)*1./(log(Vmin/videp)) rhodep = rbdum*(1.+betavol) else rhodep = rbdum endif endif rhodep=max(rhodep,50.) rhodep=min(rhodep,RHOI) gamma_arg = NU+2.+dsdum gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) gammnubet = gamma_tab(gi) !.. characteristic r-axis and a-axis after growth timestep rnf = (max((rni**2 + 2.*afn*fs/rhodep*gammnu/gammnubet*dt),QASMALL))**(0.5) anf = alphanr*rnf**(3./(2.+igr)) !.. Do not sublimation change the shape of ice from !.. prolate to oblate or vice versa phif = phii*(rnf**3/rni**3)**((igr-1.)/(igr+2.)) if(sui.lt.0.0.or.afn.lt.0.0) then if(phii.gt.1.0.and.phif.lt.1.0) then phif = phii alphanr = ani/rni anf = alphanr*rnf endif if(phii.lt.1.0.and.phif.gt.1.0) then phif = phii alphanr = ani/rni anf = alphanr*rnf endif !.. Do not let sublimation create extreme shapes if(phii.gt.1.) then if(phif.gt.phii) then phif = phii alphanr = ani/rni anf = alphanr*rnf endif else if(phif.lt.phii) then phif = phii alphanr = ani/rni anf = alphanr*rnf endif endif endif vi = fourthirdspi*rni**3*gamma_tab(gi)*i_gammnu vf = fourthirdspi*rnf**3*gamma_tab(gi)*i_gammnu rdout = rhodep rbdumtmp = rbdum*(vi/vf) + rhodep*(1.-vi/vf) rbdumtmp = min(rbdumtmp,RHOI) iwcf = nidum*rbdumtmp*vf !.. Update delta* from vapor growth/sublimation if(igr.ne.1.0)then if(anf.gt.(1.1*ao)) then dsdumout = (3.*log(rnf)-2.*log(anf)-log(ao))/ & (log(anf)-log(ao)) else dsdumout=1. endif endif !.. Do not let particles sublimate to sizes that are too small if(afn.lt.0..and.rnf.lt.1.e-6) then rbdum =RHOI rdout =RHOI dsdumout =1. phif =1. alphanr = ani/rni anf = alphanr*rnf endif if(afn.lt.0..and.anf.lt.1.e-6) then rbdum =RHOI rdout =RHOI dsdumout =1. phif =1. alphanr = ani/rni anf = alphanr*rnf endif !.. Sublimation check if(afn.lt.0..and.dsdumout.le.0.) then dsdumout=1. anf =rnf endif !.. C-axis after vapor growth gamma_arg = NU-1.+dsdumout gi=MIN((MAX((NINT((gamma_arg*100000.)-355000+1)),1)),505001) cnf = phif*anf*gammnu/gamma_tab(gi) else !.. T > T0 anf =ani cnf =cni rnf =rni iwcf=iwci endif end subroutine vaporgrow !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. This routine is used to calculate the distribution-averaged !.. capacitance (see Harrington et al. 2013) real function capacitance_gamma(ani,dsdum,NU,alphstr,i_gammnu) implicit none REAL, INTENT(IN) :: ani, dsdum, NU, alphstr, i_gammnu REAL :: a1, a2, b1, b2, c1, c2, d1, d2, gammad1, gammad2 !.. Oblate Spheroid if (dsdum.le.1.0) then a1 = 0.6369427 a2 = 0.57*a1 b1 = 0.0 b2 = 0.95 c1 = a1*alphstr**b1 c2 = a2*alphstr**b2 d1 = b1*(dsdum - 1.0) + 1.0 d2 = b2*(dsdum - 1.0) + 1.0 !.. Prolate Spheroid elseif (dsdum.gt.1.) then a1 = 0.5714285 a2 = 0.75*a1 b1 = -1.0 b2 = -0.18 c1 = a1*alphstr**(b1+1.0) c2 = a2*alphstr**(b2+1.0) d1 = b1*(dsdum - 1.0) + dsdum d2 = b2*(dsdum - 1.0) + dsdum endif if(dsdum.eq.1..and.nu.eq.1.) then gammad1 = (gamma(nu+1.0)) capacitance_gamma = ani*gammad1*i_gammnu elseif(dsdum.le.1.)then gammad1 = (gamma(nu+d1)) gammad2 = (gamma(nu+d2)) capacitance_gamma = c1*ani**d1 * gammad1*i_gammnu + & c2*ani**d2 * gammad2*i_gammnu elseif(dsdum.gt.1.) then gammad1 = (gamma(nu+d1)) gammad2 = (gamma(nu+d2)) capacitance_gamma = c1*ani**d1 * gammad1*i_gammnu + & c2*ani**d2 * gammad2*i_gammnu endif return end function capacitance_gamma !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! SUBROUTINE wet_growth_check(NU, temp, rhodum, xxlv, xxlf, qvdum, dv, kt, qs0, & fvdum, fhdum, rimedum, rni, nidum, dgflag) IMPLICIT NONE !.. Check for wet growth conditions !.. See Lamb and Verlinde (2011) REAL, INTENT(IN) :: NU, temp, rhodum, xxlv, xxlf, qvdum, dv, kt REAL, INTENT(IN) :: fvdum, fhdum, rni, nidum REAL, INTENT(IN) :: qs0, rimedum REAL :: wetg, dum LOGICAL, INTENT(INOUT) :: dgflag dum=nidum*NU*2.*rni wetg=2.*PI*dum*(kt*fhdum*(T0-temp)+rhodum*xxlv*dv*fvdum*(qs0-qvdum))/(xxlf + (CPW*(temp-T0))) if(rimedum/rhodum.gt.wetg) then dgflag = .false. endif end subroutine wet_growth_check !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! real function polysvp(T,TYPE) !-------------------------------------------------------------- !.. Taken from 'module_mp_morr_two_moment.F' (WRFV3.4) !.. COMPUTE SATURATION VAPOR PRESSURE !.. POLYSVP RETURNED IN UNITS OF PA. !.. T IS INPUT IN UNITS OF K. !.. TYPE REFERS TO SATURATION WITH RESPECT TO LIQUID (0) OR ICE (1) !.. REPLACE GOFF-GRATCH WITH FASTER FORMULATION FROM FLATAU ET AL. 1992, !.. TABLE 4 (RIGHT-HAND COLUMN) !-------------------------------------------------------------- IMPLICIT NONE REAL DUM REAL T INTEGER TYPE !.. ice real a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i data a0i,a1i,a2i,a3i,a4i,a5i,a6i,a7i,a8i /& 6.11147274, 0.503160820, 0.188439774e-1, & 0.420895665e-3, 0.615021634e-5,0.602588177e-7, & 0.385852041e-9, 0.146898966e-11, 0.252751365e-14/ !.. liquid real a0,a1,a2,a3,a4,a5,a6,a7,a8 !.. V1.7 data a0,a1,a2,a3,a4,a5,a6,a7,a8 /& 6.11239921, 0.443987641, 0.142986287e-1, & 0.264847430e-3, 0.302950461e-5, 0.206739458e-7, & 0.640689451e-10,-0.952447341e-13,-0.976195544e-15/ real dt !.. ICE IF (TYPE.EQ.1) THEN dt = max(-80.,t-273.16) polysvp = a0i + dt*(a1i+dt*(a2i+dt*(a3i+dt*(a4i+dt*(a5i+dt*(a6i+dt*(a7i+a8i*dt))))))) polysvp = polysvp*100. END IF !.. LIQUID IF (TYPE.EQ.0) THEN dt = max(-80.,t-273.16) polysvp = a0 + dt*(a1+dt*(a2+dt*(a3+dt*(a4+dt*(a5+dt*(a6+dt*(a7+a8*dt))))))) polysvp = polysvp*100. END IF end function polysvp !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Access and interpolate lookup table values (from Hugh Morrison) subroutine access_lookup_table(itabdum,dumjj,dumii,dumi,dumk,index, & dum1,dum2,dum4,dum5,proc) implicit none real :: itabdum(:,:,:,:,:), dum1,dum2,dum4,dum5,proc real :: dproc1,dproc2,iproc1,gproc1,tmp1,tmp2 integer :: dumjj,dumii,dumi,dumk,index !.. get value at current density index !.. first interpolate for current rimed fraction index dproc1 =itabdum(dumjj,dumii,dumi,dumk,index)+(dum1-real(dumi))* & (itabdum(dumjj,dumii,dumi+1,dumk,index)-itabdum(dumjj,dumii,dumi,dumk,index)) dproc2 =itabdum(dumjj,dumii,dumi,dumk+1,index)+(dum1-real(dumi))* & (itabdum(dumjj,dumii,dumi+1,dumk+1,index)-itabdum(dumjj,dumii,dumi,dumk+1,index)) iproc1=dproc1+(dum2-real(dumk))*(dproc2-dproc1) !.. linearly interpolate to get process rates for rimed fraction index + 1 dproc1 =itabdum(dumjj,dumii+1,dumi,dumk,index)+(dum1-real(dumi))* & (itabdum(dumjj,dumii+1,dumi+1,dumk,index)- & itabdum(dumjj,dumii+1,dumi,dumk,index)) dproc2 =itabdum(dumjj,dumii+1,dumi,dumk+1,index)+(dum1-real(dumi))* & (itabdum(dumjj,dumii+1,dumi+1,dumk+1,index)- & itabdum(dumjj,dumii+1,dumi,dumk+1,index)) gproc1=dproc1+(dum2-real(dumk))*(dproc2-dproc1) tmp1=iproc1+(dum4-real(dumii))*(gproc1-iproc1) !.. get value at density index + 1 !.. first interpolate for current rimed fraction index dproc1 =itabdum(dumjj+1,dumii,dumi,dumk,index)+(dum1-real(dumi))* & (itabdum(dumjj+1,dumii,dumi+1,dumk,index)- & itabdum(dumjj+1,dumii,dumi,dumk,index)) dproc2 =itabdum(dumjj+1,dumii,dumi,dumk+1,index)+ & (dum1-real(dumi))*(itabdum(dumjj+1,dumii,dumi+1,dumk+1,index)- & itabdum(dumjj+1,dumii,dumi,dumk+1,index)) iproc1=dproc1+(dum2-real(dumk))*(dproc2-dproc1) !.. linearly interpolate to get process rates for rimed fraction index + 1 dproc1 =itabdum(dumjj+1,dumii+1,dumi,dumk,index)+(dum1-real(dumi))* & (itabdum(dumjj+1,dumii+1,dumi+1,dumk,index)- & itabdum(dumjj+1,dumii+1,dumi,dumk,index)) dproc2 =itabdum(dumjj+1,dumii+1,dumi,dumk+1,index)+(dum1-real(dumi))* & (itabdum(dumjj+1,dumii+1,dumi+1,dumk+1,index)- & itabdum(dumjj+1,dumii+1,dumi,dumk+1,index)) gproc1=dproc1+(dum2-real(dumk))*(dproc2-dproc1) tmp2=iproc1+(dum4-real(dumii))*(gproc1-iproc1) !.. get final process rate proc=tmp1+(dum5-real(dumjj))*(tmp2-tmp1) end subroutine access_lookup_table !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Computes ln(gamma(xx)) using Lanczos real function gammln(xx) implicit none REAL :: xx REAL*8 :: cof(6), stp, x, y, tmp, ser INTEGER :: j DATA cof, stp/76.18009172947146d0, & -86.50532032941677d0, & 24.01409824083091d0, & -1.231739572450155d0, & 0.1208650973866179d-2, & -0.5395239384953d-5, & 2.5066282746310005d0/ x = xx y = x tmp = x + 5.5d0 tmp = (x + 0.5d0) * log(tmp) - tmp ser = 1.000000000190015d0 do j = 1, 6 y = y + 1.0d0 ser = ser + cof(j) / y enddo gammln = tmp + log(stp * ser / x) return end function gammln !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! REAL FUNCTION GAMMA(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 ! !---------------------------------------------------------------------- implicit none INTEGER I,N LOGICAL PARITY REAL & CONV,EPS,FACT,HALF,ONE,RES,SUM,TWELVE, & TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO REAL, DIMENSION(7) :: C REAL, DIMENSION(8) :: P REAL, DIMENSION(8) :: Q REAL, PARAMETER :: xxx = 0.9189385332046727417803297 !---------------------------------------------------------------------- ! MATHEMATICAL CONSTANTS !---------------------------------------------------------------------- DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/ !---------------------------------------------------------------------- ! MACHINE DEPENDENT PARAMETERS !---------------------------------------------------------------------- DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/,XINF/3.4E38/ !---------------------------------------------------------------------- ! 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/ !---------------------------------------------------------------------- ! 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/ !---------------------------------------------------------------------- ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT !---------------------------------------------------------------------- CONV(I) = REAL(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 I=1,8 XNUM=(XNUM+P(I))*Z XDEN=XDEN*Z+Q(I) END DO 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 I=1,N RES=RES*Y Y=Y+ONE END DO ENDIF ELSE !---------------------------------------------------------------------- ! EVALUATE FOR ARGUMENT .GE. 12.0, !---------------------------------------------------------------------- IF(Y.LE.XBIG)THEN YSQ=Y*Y SUM=C(7) DO I=1,6 SUM=SUM/YSQ+C(I) END DO SUM=SUM/Y-Y+xxx 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 RETURN ! ---------- LAST LINE OF GAMMA ---------- END FUNCTION GAMMA !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Get the inherent growth ratio data as a function of temperature real function get_igr(igrdatadum, temp) implicit none real :: igrdatadum(:), dum, igr1, igr2, temp !.. Inherent growth ratio (IGR) from data !.. See Lamb and Scott (1972) and Chen and Lamb (1994) if((temp-T0).ge.-59..and.(temp-T0).le.-1.) then dum = (abs(real(int(temp-T0))) + 1.) - abs(temp-T0) igr1 = igrdatadum(max((int(temp-T0)*(-1)),1)) igr2 = igrdatadum(min(((int(temp-T0)*(-1))+1),60)) get_igr = dum*igr1 + (1.-dum)*igr2 elseif((temp-T0).gt.-1..and.(temp-T0).le.0.) then dum = 1.0 - abs(temp-T0) igr1 = 1.0 igr2 = igrdatadum(1) get_igr = dum*igr1 + (1.-dum)*igr2 elseif((temp-T0).ge.-60..and.(temp-T0).lt.-59.) then dum = 60.0 - abs(temp-T0) igr1 = igrdatadum(59) igr2 = igrdatadum(60) get_igr = dum*igr1 + (1.-dum)*igr2 elseif((temp-T0).lt.-60.) then get_igr = igrdatadum(60) else get_igr = 1. endif end function get_igr !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Aggregation routine (see Jensen et al. 2017, implentation by JYH and AAJ) !.. See Walko et al. 1995 and Meyers et al. 1997 for more details subroutine aggregation(dtdum,rhoairdum,tempdum,qdum1,ndum1,ddum1,qdum2,ndum2,ddum2, & qdum3,ndum3,ddum3,rhodum1,rhodum2,phidum1,phidum2,coltab,coltabn, & qagg1,qagg2,qagg3,nagg1,nagg2,nagg3) implicit none real, parameter :: QNSMALL= 1.25e-7 ! Smallest ice number integer icat,k,kp1,kp2,isnow,ipris,it,ihab,jhab integer, parameter :: ncat=7,npair=35,ndn=60,nftab=1000,ngam=5000,nz=1 real, parameter :: rsmall=1.e-12 real gnu(ncat),shape(ncat),cfmas(ncat),pwmas(ncat),& cfden(ncat),pwden(ncat)& ,cfvt(ncat),pwvt(ncat),tablo(ncat),tabhi(ncat),tabloln(ncat),& tabhiln(ncat)& ,dnmin(ncat),dnmax(ncat),vtfac(ncat),pow1(ncat),pow2(ncat),& frefac1(ncat)& ,frefac2(ncat),cfmasft(ncat),dict(ncat),dbmi(ncat),& gamm(ncat),gamn1(ncat)& ,gam(ngam,ncat),ncrossloss,rcrossloss,nselfloss,& rselfloss,ncrossgain,rcrossgain,nselfgain,rselfgain real rictmin,rictmax,rict,rictmm,dtlt,totalc(nz),totaldc(nz) real dn(nz,ncat),en(nz,ncat),qr(nz,ncat),r(nz,ncat),qq(nz,ncat) real wct1(nz,ncat),wct2(nz,ncat),rxfer(nz,ncat,ncat),& qrxfer(nz,ncat,ncat),enxfer(nz,ncat,ncat),& parxfer(nz,ncat,ncat), paqrxfer(nz,ncat,ncat),& paenxfer(nz,ncat,ncat), pprxfer(nz,ncat,ncat),& ppqrxfer(nz,ncat,ncat), ppenxfer(nz,ncat,ncat) real qs(nz),ns(nz),dns(nz),qa(nz),na(nz),dna(nz) real t(nz,ncat),tc(nz),rhoa(nz),tk(nz),denfac(nz),eff(nz) real ftable(nftab,ncat) real, intent(IN) :: coltab(ndn,ndn,npair),coltabn(ndn,ndn,npair) integer ipair(ncat,ncat),ict(nz,ncat) real, dimension(ncat,9) :: dstprms real, intent(IN) :: qdum1,ndum1,qdum2,ndum2,qdum3,ndum3,ddum1,ddum2 real, intent(INOUT) :: qagg1,qagg2,qagg3,nagg1,nagg2,nagg3,ddum3 real, intent(IN) :: rhodum1, rhodum2, phidum1, phidum2 real :: rhoeff, phieff, rhoeffmax, phieffmax, efffact real :: source3, sink3, source4, sink4, ratioagg real, intent(IN) :: dtdum, rhoairdum, tempdum dstprms=& !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! cloud rain planar columnar aggreg graup hail ! 1 2 3 4 5 6 7 ---> ncat !----------------------------------------------------------------------- reshape((/.5, .5, .318, .318, .5, .5, .5, &!shape (capacitance shape) 524., 524., .333, .333, .496, 157., 471., &!cfmas (m=cfmas*D**pwmas) 3., 3., 2.4, 2.4, 2.4, 3., 3., &!pwmas 3173., 149., 4.836, 4.836, 3.084, 93.3, 161., &!cfvt (v = cfvt*D**pwvt) 2., .5, 0.25, .25, .2, .5, .5, &!pwvt 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, 1.e-6, &!tablo (size limits for table) 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, 1.e-2, &!tabhi .1e-5, .1e-4, .1e-5, .1e-5, .1e-4, .1e-4, .1e-4, &!dnmin (characteristic diam .1e-2, .1e-1, .1e-1, .1e-1, .1e-1, .1e-1, .1e-1/), &!dnmax limits) (/ncat,9/)) !----------------------------------------------------------------------- ipair=& reshape((/0,0,3,7,11,15,19 ,1,0,5,9,13,17,21 & ! collection pairs ,2,4,22,24,0,0,0 ,6,8,23,28,0,0,0 & ,10,12,25,29,35,0,0 ,14,16,26,30,32,0,0 & ,18,20,27,31,33,34,0/), & (/ncat,ncat/)) do icat=1,ncat shape(icat)=dstprms(icat,1) cfmas(icat)=dstprms(icat,2) pwmas(icat)=dstprms(icat,3) cfvt (icat)=dstprms(icat,4) pwvt (icat)=dstprms(icat,5) tablo(icat)=dstprms(icat,6) tabhi(icat)=dstprms(icat,7) dnmin(icat)=dstprms(icat,8) dnmax(icat)=dstprms(icat,9) gnu(icat) = 4.0 ! shape distribution tabloln(icat)=log(tablo(icat)) tabhiln(icat)=log(tabhi(icat)) dict(icat) = float(ndn - 1) / (tabhiln(icat) - tabloln(icat)) enddo rictmin=1.0001 rictmax=0.9999*float(ndn) !.. Collection integral multipliers for number and mass mixing ratio loss !.. due to self collection from a cateogry. Multipliers for the gain term !.. for aggregates given as well. Below this are the multipliers for !.. cross collection. These should be correct but a more detailed !.. mathematical treatment is warranted. JYH Sept 2016. nselfloss = 1. rselfloss = 1. nselfgain = 1./2. rselfgain = 1.0 ncrossloss = 1. rcrossloss = 1.0 ncrossgain = 1./2. rcrossgain = 1.0 !.. Assumed mass-size coefficients pwmas(3) = 3.0 cfmas(3) = 3.14159/6.*250.*0.05 pwmas(4) = 3.0 cfmas(4) = 3.14159/6.*250.*0.05 pwmas(5) = 3.0 cfmas(5) = 3.14159/6.*50.*0.2 dtlt = dtdum ! time step kp1 = 1 kp2 = nz ipris= 5 ! double moment !.. collection routine !.. nz = # vertical levels !.. kp1 = lowest level with cloud !.. kp2 = higest level with cloud !.. 3 => the first is the collector !.. 3 => the second is the collectee !.. 5 => where the collection goes (aggregates) !.. dn(nz,ncat) = characteristic diameter !.. en(nz,ncat) = number concentration !.. r(nz,ncat) = mass mixing ratio !.. qr(nz,ncat) = hydrometeor internal enthalpy !.. qq(nz,ncat) = hydrometeor temperature (internal energy) !.. t(nz,ncat) = hydrometeor temperature (air temperature) !.. ict,wct1,wct2 -> table interpolation !.. rhoa(nz) = air density !.. denfac(nz) = sqrt(1/rhoa) !.. eff(nz) = collection efficiency !.. ipris = 5 double moment !.. rxfer(nz,ncat,ncat) = change in mass mixing ratio due to collection !.. enxfer(nz,ncat,ncat) = change in number concentration due to collection !.. qrxfer(nz,ncat,ncat) = change in internal energy due to collection !.. ipair,2 -> 2 is for self-collection, 1 if cross collection do k = 1,nz !.. nz = 1 currently (the k loop is done outsize of this routine) rhoa(k) = rhoairdum tk(k) = tempdum denfac(k) = sqrt(1./rhoa(k)) dn(k,3) = ddum1 en(k,3) = ndum1*rhoa(k) r(k,3) = qdum1 dn(k,4) = ddum2 en(k,4) = ndum2*rhoa(k) r(k,4) = qdum2 qr(k,3) = 0.0 !.. internal energy... not needed so pass in zero qq(k,3) = t(k,3) t(k,3) = tk(k) - 273.15 qr(k,4) = 0.0 !.. internal energy... not needed so pass in zero qq(k,4) = t(k,3) t(k,4) = t(k,3) qr(k,5) = 0.0 qq(k,5) = t(k,5) t(k,5) = t(k,3) dn(k,5) = ddum3 en(k,5) = ndum3*rhoa(k) r(k,5) = qdum3 totaldc(k)=0. enddo do it = 1,1 !.. One timestep at a time do k = 1,nz !.. One vertical grid at a time do ihab = 1,7 do jhab = 1,7 !.. planar = qice, columnar=qice2, aggregates = qice3 !.. planar + columnar = aggregates (qice + qice2 -> qice3) rxfer(k,ihab,jhab) = 0.0 enxfer(k,ihab,jhab) = 0.0 qrxfer(k,ihab,jhab) = 0.0 rxfer(k,ihab,jhab) = 0.0 enxfer(k,ihab,jhab) = 0.0 qrxfer(k,ihab,jhab) = 0.0 !.. planar + aggregates = aggregates parxfer(k,ihab,jhab) = 0.0 paenxfer(k,ihab,jhab) = 0.0 paqrxfer(k,ihab,jhab) = 0.0 !.. columnar + aggregates = aggregates parxfer(k,ihab,jhab) = 0.0 paenxfer(k,ihab,jhab) = 0.0 paqrxfer(k,ihab,jhab) = 0.0 !.. Self collection (aggregates) pprxfer(k,ihab,jhab) = 0.0 ppenxfer(k,ihab,jhab) = 0.0 ppqrxfer(k,ihab,jhab) = 0.0 pprxfer(k,ihab,jhab) = 0.0 ppenxfer(k,ihab,jhab) = 0.0 ppqrxfer(k,ihab,jhab) = 0.0 enddo enddo !.. needed for table interpolation rict=dict(3)*(log(max(1.e-10,dn(k,3)))-tabloln(3))+1. rictmm=max(rictmin,min(rictmax,rict)) ict(k,3)=int(rictmm) wct2(k,3)=rictmm-float(ict(k,3)) wct1(k,3)=1.0-wct2(k,3) !.. needed for table interpolation rict=dict(4)*(log(max(1.e-10,dn(k,4)))-tabloln(4))+1. rictmm=max(rictmin,min(rictmax,rict)) ict(k,4)=int(rictmm) wct2(k,4)=rictmm-float(ict(k,4)) wct1(k,4)=1.0-wct2(k,4) !.. needed for table interpolation rict=dict(5)*(log(max(1.e-10,dn(k,5)))-tabloln(5))+1. rictmm=max(rictmin,min(rictmax,rict)) ict(k,5)=int(rictmm) wct2(k,5)=rictmm-float(ict(k,5)) wct1(k,5)=1.0-wct2(k,5) totalc(k)=0. enddo !.. planar + columnar = aggregates !.. let high density reduce the collection efficiency rhoeffmax = max(rhodum1,rhodum2) phieffmax = max(min(phidum1,1./phidum1),min(phidum2,1./phidum2)) if(rhoeffmax.le.400.) then rhoeff = 1. else rhoeff = -0.001923*rhoeffmax + 1.76916 endif rhoeff = max(rhoeff,0.) rhoeff = min(rhoeff,1.) !.. Spherical particles collect less efficiently than more extreme !.. particle shapes !.. Based on Connolly et al 2012 !.. See Jensen et al. 2018b if(phieffmax.le.0.03) then phieff = 1. elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then phieff = 0.0001*((phieffmax+0.07)**(-4.)) else phieff = 0. endif if(phieff.le.0.001) then phieff = 0. endif phieff = max(phieff,0.01) !.. Low limit currently set to 0.01 AAJ phieff = min(phieff,1.) efffact = phieff * rhoeff call col1 & (nz,kp1,kp2,ndn,ncat,3,4,5 & ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab & ,coltabn,ipair,1.,eff,efffact,ipris,rxfer,qrxfer,enxfer) call col1 & (nz,kp1,kp2,ndn,ncat,4,3,5 & ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab & ,coltabn,ipair,1.,eff,efffact,ipris,rxfer,qrxfer,enxfer) !.. plates + aggregates = aggregates rhoeffmax = rhodum1 phieffmax = min(phidum1,1./phidum1) if(rhoeffmax.le.400.) then rhoeff = 1. else rhoeff = -0.001923*rhoeffmax + 1.76916 endif rhoeff = max(rhoeff,0.) rhoeff = min(rhoeff,1.) !.. Based on Connolly et al 2012 !.. See Jensen et al. 2018b if(phieffmax.le.0.03) then phieff = 1. elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then phieff = 0.0001*((phieffmax+0.07)**(-4.)) else phieff = 0. endif if(phieff.le.0.001) then phieff = 0. endif phieff = max(phieff,0.01) phieff = min(phieff,1.) efffact = phieff * rhoeff call col1 & (nz,kp1,kp2,ndn,ncat,3,5,5 & ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab & ,coltabn,ipair,1.,eff,efffact,ipris,parxfer,paqrxfer,paenxfer) !.. Columnar + Aggregates = Aggregates rhoeffmax = rhodum2 phieffmax = min(phidum2,1./phidum2) if(rhoeffmax.le.400.) then rhoeff = 1. else rhoeff = -0.001923*rhoeffmax + 1.76916 endif rhoeff = max(rhoeff,0.) rhoeff = min(rhoeff,1.) !.. Based on Connolly et al 2012 !.. See Jensen et al. 2018b if(phieffmax.le.0.03) then phieff = 1. elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then phieff = 0.0001*((phieffmax+0.07)**(-4.)) else phieff = 0. endif if(phieff.le.0.001) then phieff = 0. endif phieff = max(phieff,0.01) phieff = min(phieff,1.) efffact = phieff * rhoeff call col1 & (nz,kp1,kp2,ndn,ncat,4,5,5 & ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab & ,coltabn,ipair,1.,eff,efffact,ipris,parxfer,paqrxfer,paenxfer) !.. Plates + Plates = Aggregates rhoeffmax = rhodum1 phieffmax = min(phidum1,1./phidum1) if(rhoeffmax.le.400.) then rhoeff = 1. else rhoeff = -0.001923*rhoeffmax + 1.76916 endif rhoeff = max(rhoeff,0.) rhoeff = min(rhoeff,1.) !.. Based on Connolly et al 2012 !.. See Jensen et al. 2018b if(phieffmax.le.0.03) then phieff = 1. elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then phieff = 0.0001*((phieffmax+0.07)**(-4.)) else phieff = 0. endif if(phieff.le.0.001) then phieff = 0. endif phieff = max(phieff,0.01) phieff = min(phieff,1.) efffact = phieff * rhoeff call col1 & (nz,kp1,kp2,ndn,ncat,3,3,5 & ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab & ,coltabn,ipair,1.,eff,efffact,ipris,pprxfer,ppqrxfer,ppenxfer) !.. Columns + Columns = Aggregates rhoeffmax = rhodum2 phieffmax = min(phidum2,1./phidum2) if(rhoeffmax.le.400.) then rhoeff = 1. else rhoeff = -0.001923*rhoeffmax + 1.76916 endif rhoeff = max(rhoeff,0.) rhoeff = min(rhoeff,1.) !.. Based on Connolly et al 2012 !.. See Jensen et al. 2018b if(phieffmax.le.0.03) then phieff = 1. elseif(phieffmax.gt.0.03.and.phieffmax.lt.0.5) then phieff = 0.0001*((phieffmax+0.07)**(-4.)) else phieff = 0. endif if(phieff.le.0.001) then phieff = 0. endif phieff = max(phieff,0.01) phieff = min(phieff,1.) efffact = phieff * rhoeff call col1 & (nz,kp1,kp2,ndn,ncat,4,4,5 & ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab & ,coltabn,ipair,1.,eff,efffact,ipris,pprxfer,ppqrxfer,ppenxfer) efffact=1. !.. Aggregate self collection call col1 & (nz,kp1,kp2,ndn,ncat,5,5,5 & ,dtlt,dn,en,r,qr,qq,t,ict,wct1,wct2,rhoa,denfac,coltab & ,coltabn,ipair,1.,eff,efffact,ipris,pprxfer,ppqrxfer,ppenxfer) do k = 1,nz !.. Do not over-deplete from aggregation sink3 = rcrossloss*rxfer(k,3,5) + & rcrossloss*parxfer(k,3,5) + & rselfloss*pprxfer(k,3,5) source3 = r(k,3) sink4 = rcrossloss*rxfer(k,4,5) + & rcrossloss*parxfer(k,4,5) + & rselfloss*pprxfer(k,4,5) source4 = r(k,4) if(sink3.gt.source3.and.source3.gt.1.e-8) then ratioagg= source3/sink3 rxfer(k,3,5) = rxfer(k,3,5)*ratioagg parxfer(k,3,5) = parxfer(k,3,5)*ratioagg pprxfer(k,3,5) = pprxfer(k,3,5)*ratioagg endif if(sink4.gt.source4.and.source4.gt.1.e-8) then ratioagg= source4/sink4 rxfer(k,4,5) = rxfer(k,4,5)*ratioagg parxfer(k,4,5) = parxfer(k,4,5)*ratioagg pprxfer(k,4,5) = pprxfer(k,4,5)*ratioagg endif totaldc(k)=totaldc(k) + & ncrossloss*enxfer(k,3,5)+ncrossloss*paenxfer(k,3,5)+ & nselfloss*ppenxfer(k,3,5)+ & ncrossloss*enxfer(k,4,5)+ncrossloss*paenxfer(k,4,5)+ & nselfloss*ppenxfer(k,4,5) !.. qagg is the mass tendency, nagg is the number tendency qagg1 = -rcrossloss*rxfer(k,3,5) - rcrossloss*parxfer(k,3,5) - rselfloss*pprxfer(k,3,5) nagg1 = -ncrossloss*enxfer(k,3,5) - ncrossloss*paenxfer(k,3,5) - nselfloss*ppenxfer(k,3,5) qagg2 = -rcrossloss*rxfer(k,4,5) - rcrossloss*parxfer(k,4,5) - rselfloss*pprxfer(k,4,5) nagg2 = -ncrossloss*enxfer(k,4,5) - ncrossloss*paenxfer(k,4,5) - nselfloss*ppenxfer(k,4,5) qagg3 = rcrossgain*rxfer(k,3,5) + rcrossgain*rxfer(k,4,5) + rcrossgain*parxfer(k,3,5) + & rcrossgain*parxfer(k,4,5) + rselfgain*pprxfer(k,3,5) + rselfgain*pprxfer(k,4,5) nagg3 = ncrossgain*enxfer(k,3,5) + ncrossgain*enxfer(k,4,5) + nselfgain*ppenxfer(k,3,5) + & nselfgain*ppenxfer(k,4,5) - 0.5*ppenxfer(k,5,5) qagg1 = min(qagg1,0.) qagg2 = min(qagg2,0.) qagg3 = max(qagg3,0.) !.. Update to characteristic diameters after aggregation !.. Note that for ice-one and ice-two, the axes will be re-diagnosed !.. after aggregation based on delta*, so only output aggregation size dn(k,5) = ((r(k,5)+qagg3)*rhoa(k)/(MAX((en(k,5)+nagg3),QNSMALL)*cfmas(5) & *(gamma(gnu(5)+ & pwmas(5)))/(gamma(gnu(5))) ) )**(1./pwmas(5)) ddum3=dn(k,5) nagg1 = nagg1/rhoa(k) nagg2 = nagg2/rhoa(k) nagg3 = nagg3/rhoa(k) totalc(k)=0. totalc(k) = en(k,3)+en(k,4) + 0. rxfer(k,3,5) = 0.0 enxfer(k,3,5) = 0.0 qrxfer(k,3,5) = 0.0 rxfer(k,4,5) = 0.0 enxfer(k,4,5) = 0.0 qrxfer(k,4,5) = 0.0 enddo enddo end subroutine aggregation !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Needed for the aggregation routine subroutine col1(n1,k1,k2,ndn,ncat,mx,my,mz,dtlt,dn,en,r,qr,qq,t & ,ict,wct1,wct2,dn0,denfac,coltab,coltabn,ipair,cr,eff,efdum,icf & ,rxfer,qrxfer,enxfer) integer ict(n1,*),ipair(ncat,ncat), ndn, ncat,n1 real dn(n1,*),en(n1,*),r(n1,*),qr(n1,*),qq(n1,*),t(n1,*) & ,wct1(n1,*),wct2(n1,*),dn0(*),denfac(*) & ,coltab(ndn,ndn,*),eff(*),coltabn(ndn,ndn,*) & ,rxfer(n1,ncat,ncat),qrxfer(n1,ncat,ncat),enxfer(n1,ncat,ncat) integer :: k1, k2, mx, my, mz, icf, k, ip real :: dtlt, cr, csizei, tmax, colamt, colamtn, wght, deltan,efdum ip=ipair(mx,my) csizei=7.e4 do k=k1,k2 tmax=max(t(k,mx),t(k,my)) !.. efdum added by AAJ from ramsagg !.. if takes into account the shape and density !.. dependence on the collection efficiency eff(k)=min(0.2,10.**(0.035*tmax-0.7))*efdum if(my.eq.6.and.qr(k,6).gt.0)eff(k)=1. if(my.eq.7.and.qr(k,7).gt.0)eff(k)=1. !.. see note above on efdum if(mz.eq.5.and.abs(t(k,mx)+14.).le.2.)eff(k)=1.4*efdum if(mx.eq.1)eff(k)=min(1.,dn(k,1)*csizei) colamt=dtlt*cr*0.785*eff(k)*denfac(k)/dn0(k) & *en(k,mx)*en(k,my)*( & wct1(k,mx)*wct1(k,my)*coltab(ict(k,mx) ,ict(k,my) ,ip) & +wct2(k,mx)*wct1(k,my)*coltab(ict(k,mx)+1,ict(k,my) ,ip) & +wct1(k,mx)*wct2(k,my)*coltab(ict(k,mx) ,ict(k,my)+1,ip) & +wct2(k,mx)*wct2(k,my)*coltab(ict(k,mx)+1,ict(k,my)+1,ip)) colamt=min(colamt,r(k,mx)) rxfer(k,mx,mz) = rxfer(k,mx,mz) + colamt qrxfer(k,mx,mz) = qrxfer(k,mx,mz) + colamt * qq(k,mx) colamtn=dtlt*cr*0.785*eff(k)*denfac(k)/dn0(k) & *en(k,mx)*en(k,my)*( & wct1(k,mx)*wct1(k,my)*coltabn(ict(k,mx) ,ict(k,my) ,ip) & +wct2(k,mx)*wct1(k,my)*coltabn(ict(k,mx)+1,ict(k,my) ,ip) & +wct1(k,mx)*wct2(k,my)*coltabn(ict(k,mx) ,ict(k,my)+1,ip) & +wct2(k,mx)*wct2(k,my)*coltabn(ict(k,mx)+1,ict(k,my)+1,ip)) colamtn=min(colamtn,en(k,mx)) if(icf.eq.5)then wght = min(max(0.,100.*abs(colamt)/max(1.e-20,r(k,mx))),1.) deltan = colamtn*(1.-wght) + & wght*en(k,mx) * colamt / max(1.e-20,r(k,mx)) enxfer(k,mx,mz) = enxfer(k,mx,mz) + deltan endif enddo return end subroutine col1 !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Needed for the aggregation routine subroutine mkcoltb(ndn,ncat,coltab,coltabn,ipair,gnu,tablo,tabhi & ,cfmas,pwmas,cfvt,pwvt) integer :: ncat real :: coltab(ndn,ndn,35),gnu(ncat) & ,tablo(ncat), tabhi(ncat), cfmas(ncat),pwmas(ncat) & ,cfvt(ncat), pwvt(ncat), coltabn(ndn,ndn,35) integer :: ipair(ncat,ncat), ix, iy, idny, idnx, idx, ndn integer, parameter :: ndx=20 real :: dx(ndx),gx(ndx),fx(20), gyn,gyn1,gyn2,gynp1,gynp2,gynp real :: dxhi, dny, vny, dxlo, dnx, ans do ix=1,ncat do iy=1,ncat if(ipair(ix,iy).gt.0)then gyn=(gamma(gnu(iy))) gyn1=(gamma(gnu(iy)+1.))/gyn gyn2=(gamma(gnu(iy)+2.))/gyn gynp=(gamma(gnu(iy)+pwvt(iy)))/gyn gynp1=(gamma(gnu(iy)+pwvt(iy)+1.))/gyn gynp2=(gamma(gnu(iy)+pwvt(iy)+2.))/gyn dxlo=tablo(ix)*0.01 dxhi=tabhi(ix)*10.0 do idny=1,ndn dny=tablo(iy)*(tabhi(iy)/tablo(iy))**(float(idny-1) & /float(ndn-1)) vny=cfvt(iy)*dny**pwvt(iy) do idnx=1,ndn dnx=tablo(ix)*(tabhi(ix)/tablo(ix))**(float(idnx-1) & /float(ndn-1)) do idx=1,ndx dx(idx)=dxlo*(dxhi/dxlo) & **(float(idx-1)/float(ndx-1)) fx(idx) = xjnum(dx(idx),cfvt(ix),pwvt(ix),cfvt(iy) & ,pwvt(iy),vny,dnx,dny,gnu(ix),gnu(iy) & ,gyn1,gyn2,gynp,gynp1,gynp2) if(fx(idx).lt.1.e-15) then !.. low limit added by AAJ fx(idx) = 0. endif gx(idx) = fx(idx) * cfmas(ix) & * dx(idx) ** pwmas(ix) enddo call avint(dx,gx,ndx,dxlo,dxhi,ans) coltab(idnx,idny,ipair(ix,iy))=max(0.,ans) call avint(dx,fx,ndx,dxlo,dxhi,ans) coltabn(idnx,idny,ipair(ix,iy))=max(0.,ans) enddo enddo endif enddo enddo return end subroutine mkcoltb !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Needed for the aggregation routine function xjnum(dx,cvx,pvx,cvy,pvy,vny,dnx,dny,xnu,ynu & ,gyn1,gyn2,gynp,gynp1,gynp2) implicit none real :: dx,cvx,pvx,cvy,pvy,vny,dnx,dny,xnu,ynu,gyn1,gyn2,gynp, & gynp1,gynp2,xjnum & ,dnxi,rdx,vx,dxy,ynup dnxi = 1. / dnx rdx = dx * dnxi vx = max((cvx * dx ** pvx),1.e-6) !.. low limit added by AAJ dxy = (vx / cvy) ** (1. / pvy) / dny dxy = max(dxy,1.e-5) dxy = min(dxy,70.) ynup = ynu + pvy if (rdx .lt. 38.) then xjnum=exp(-rdx-gammln(xnu)-gammln(ynu))*rdx**(xnu-1.)*dnxi*( & vx*(dx*dx*(gammap(ynu,dxy)-gammq(ynu,dxy)) & +2.*dx*dny*gyn1*(gammap(ynu+1.,dxy)-gammq(ynu+1.,dxy)) & +dny*dny*gyn2*(gammap(ynu+2.,dxy)-gammq(ynu+2.,dxy))) & -vny*(dx*dx*gynp*(gammap(ynup,dxy)-gammq(ynup,dxy)) & +2.*dx*dny*gynp1*(gammap(ynup+1.,dxy)-gammq(ynup+1.,dxy)) & +dny*dny*gynp2*(gammap(ynup+2.,dxy)-gammq(ynup+2.,dxy)))) else xjnum=0. endif return end function xjnum !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Used by aggregation REAL FUNCTION GAMMQ(A,X) real :: x, a, gln, gamser IF (X .LT. A+1.) THEN CALL LOWGSERIES(GAMSER,A,X,GLN) GAMMQ = 1. - GAMSER ELSE CALL HIGHGCONTFRAC(GAMMQ,A,X,GLN) END IF RETURN END FUNCTION GAMMQ !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Computes P(a, x) = lower_gamma(a, x) / gamma(a) !.. Using either series solution or continued fractions real function gammap(a, x) implicit none REAL :: a, x, gammaser, gammacf, gln if(x.lt.0.0.or.a.lt.0.0) then write(*,*) 'either x or a is less than zero' write(*,*) x, a return endif if(x.lt.(a + 1.0)) then call lowgseries(gammaser, a, x, gln) gammap = gammaser else call highgcontfrac(gammacf, a, x, gln) gammap = 1.0 - gammacf endif return end function gammap !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Returns: gln = ln(gamma(a)) !.. Returns: gammaser = lower_gamma(a, x) as a series solution subroutine lowgseries(gammaser, a, x, gln) implicit none INTEGER, PARAMETER :: ITMAX = 100 INTEGER :: n REAL :: gammaser, a, x, gln, sum, del, ap REAL, PARAMETER :: EPS = 3.0e-7 gln = gammln(a) if(x.le.0) then gammaser = 0.0 return endif 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) exit enddo gammaser = sum*exp(-x+a*log(x)-gln) return end subroutine lowgseries !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Returns: gln = ln(gamma(a)) !.. Returns: gammacf = upper_gamma(a, x) as a continued fractions solution subroutine highgcontfrac(gammacf, a, x, gln) implicit none INTEGER, PARAMETER :: ITMAX = 100 INTEGER :: n REAL :: gammacf, a, x, gln, b, c, d, h, del, an REAL, PARAMETER :: EPS = 3.0e-7, FPMIN = 1.0e-30 gln = gammln(a) b=x+1.-a c=1./FPMIN d=1./b h=d do n = 1, ITMAX an=-n*(n-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) exit enddo gammacf = exp(-x+a*log(x)-gln)*h return end subroutine highgcontfrac !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! !.. Used by aggregation SUBROUTINE AVINT(X,Y,N,XLO,XUP,ANS) DOUBLE PRECISION :: R3,RP5,SUM,SYL,SYL2,SYL3,SYU,SYU2,SYU3,X1,X2,X3 & ,X12,X13,X23,TERM1,TERM2,TERM3,A,B,C,CA,CB,CC REAL :: XLO,XUP,ANS,FL,FR,SLOPE INTEGER :: N, I, INLFT,INRT,ISTART,ISTOP REAL, DIMENSION(N) :: X,Y ANS =0.0 IF (XLO-XUP) 3,100,200 3 IF (N.LT.2) GO TO 215 DO 5 I=2,N IF (X(I).LE.X(I-1)) GO TO 210 IF (X(I).GT.XUP) GO TO 6 5 CONTINUE 6 CONTINUE IF (N.GE.3) GO TO 9 !.. special n=2 case SLOPE = (Y(2)-Y(1))/(X(2)-X(1)) FL = Y(1) + SLOPE*(XLO-X(1)) FR = Y(2) + SLOPE*(XUP-X(2)) ANS = 0.5*(FL+FR)*(XUP-XLO) RETURN 9 CONTINUE IF (X(N-2).LT.XLO) GO TO 205 IF (X(3).GT.XUP) GO TO 205 I = 1 10 IF (X(I).GE.XLO) GO TO 15 I = I+1 GO TO 10 15 INLFT = I I = N 20 IF (X(I).LE.XUP) GO TO 25 I = I-1 GO TO 20 25 INRT = I IF ((INRT-INLFT).LT.2) GO TO 205 ISTART = INLFT IF (INLFT.EQ.1) ISTART = 2 ISTOP = INRT IF (INRT.EQ.N) ISTOP = N-1 R3 = 3.0D0 RP5= 0.5D0 SUM = 0.0 SYL = XLO SYL2= SYL*SYL SYL3= SYL2*SYL DO 50 I=ISTART,ISTOP X1 = X(I-1) X2 = X(I) X3 = X(I+1) X12 = X1-X2 X13 = X1-X3 X23 = X2-X3 TERM1 = DBLE(Y(I-1))/(X12*X13) TERM2 =-DBLE(Y(I)) /(X12*X23) TERM3 = DBLE(Y(I+1))/(X13*X23) A = TERM1+TERM2+TERM3 B = -(X2+X3)*TERM1 - (X1+X3)*TERM2 - (X1+X2)*TERM3 C = X2*X3*TERM1 + X1*X3*TERM2 + X1*X2*TERM3 IF (I-ISTART) 30,30,35 30 CA = A CB = B CC = C GO TO 40 35 CA = 0.5*(A+CA) CB = 0.5*(B+CB) CC = 0.5*(C+CC) 40 SYU = X2 SYU2= SYU*SYU SYU3= SYU2*SYU SUM = SUM + CA*(SYU3-SYL3)/R3 + CB*RP5*(SYU2-SYL2) + CC*(SYU-SYL) CA = A CB = B CC = C SYL = SYU SYL2= SYU2 SYL3= SYU3 50 CONTINUE SYU = XUP ANS = SUM + CA*(SYU**3-SYL3)/R3 + CB*RP5*(SYU**2-SYL2) & + CC*(SYU-SYL) 100 RETURN 200 PRINT*, 'Upper limit of integration not greater than lower limit.' STOP 'AVINT2' 205 PRINT*, 'Less than 3 function values between integration limits.' STOP 'AVINT3' 210 PRINT*, 'Abscissas not strictly increasing.' STOP 'AVINT4' 215 PRINT*, 'Less than 2 function values were supplied.' STOP 'AVINT5' END SUBROUTINE !--------------------------------------------------------------------------------------------------------------! !--------------------------------------------------------------------------------------------------------------! end module module_mp_jensen_ishmael !.. Full stop