! The National Taiwan University (NTU) scheme is a multi-moment bulk scheme (two moments for the liquid-phase and 3 moments ! for the ice-phase hydrometeors) with consideration for ice crystal shape and density variations. This scheme is developed ! by Dr. Tzu-Chin Tsai and Prof. Jen-Ping Chen at the National Taiwan University. The NTU multi-moment bulk microphysics scheme ! default version implemented into the WRF model by Tzu-Chin Tsai at NTU (d97229002@ntu.edu.tw; tzuchin12@gmail.com) ! Realsed on 2021/1/6 !------------------------------------------------------------------------------------------------------------------------------- ! References : ! 1. Chen, J.-P. and S.-T. Liu, 2004: Physically based two-moment bulkwater parametrization for warm-cloud microphysics. ! Quart. J. Roy. Meteor. Soc., 130, 51-78. ! 2. Cheng, C.-T., W.-C. Wang, and J.-P. Chen, 2007: A modeling study of aerosol impacts on cloud microphysics and radiative ! properties. Quart. J. Roy. Meteor. Soc., 133, 283-297. ! 3. Cheng, C.-T., W.-C. Wang, and J.-P. Chen, 2010: Simulation of the effects of increasing cloud condensation nuclei on ! mixed-phase clouds and precipitation of a front system. Atmos. Res., 96, 461-476. ! 4. Chen, J.-P. and T.-C. Tsai, 2016: Triple-moment modal parameterization for the adaptive growth habit of pristine ice ! crystals. J. Atmos. Sci., 73, 2105-2122. ! 5. Tsai, T.-C., and J.-P. Chen, 2020: Multi-moment ice bulk microphysics scheme with consideration for particle shape ! and apparent density. Part I: Methodology and idealized simulation. J. Atmos. Sci., 77, 1821-1850. ! !-------------------------------------------------------------------------------------------------------------------------------- ! The notation for terms involving two interacting categories are denoted by B(x,y)MPxy, where B is the prognostic variable ! MP represents the microphysical processes, and the subscripts x and y indicate the quantity B being transferred from the ! hydrometeor category x to y. If terms involve three interacting categories, then a three-letter subscript of xyz would ! follow the microphysical process with the destination of hydrometeor category z. Another notation denoted by QHomp is ! only used for the change rate of the hailstone mixing ratio, and o indicates whether the dry (denote d) or wet (denote w) ! growth mode occurs. Please referred to supplemental information related to the Tsai and Chen (2020) is available at the ! Journals Online website: https://doi.org/10.1175/JAS-D-19-0125.s1. ! !================================================================================================================================= MODULE module_mp_ntu USE module_wrf_error IMPLICIT NONE PUBLIC :: MP_NTU PRIVATE :: GAMMA,GAMLN,GAMIN,GAMMP,GSER,CFG,GUESS_RC,YEQU,DYEQU, & PDF,DPDF,DLNX,POLYSVP INTEGER, PRIVATE, PARAMETER :: ID_NH42SO4 = 1, &! 1. (NH4)2SO4 ID_DUST = 0, &! 2. DUST ID_IN = 2 ! 3. IN INTEGER, PRIVATE, PARAMETER :: INSPEC = 2 ! 1. SOOT; 2. SD; 3. AD; 4. BT INTEGER, PRIVATE, PARAMETER :: ICE_SHAPE = 1 ! 0. SPHERICAL; 1. ASPHERICAL; PRISTINE ICE (HABIT) INTEGER, PRIVATE, PARAMETER :: AGG_SHAPE = 1 ! 0. SPHERICAL; 1. CL94; ASPHERICAL AGGREGATES INTEGER, PRIVATE, PARAMETER :: ICE_RHOI = 1 ! 0. FIXED; 1. PREDICT; 2. FIXED; PRISTINE ICE DENSITY INTEGER, PRIVATE, PARAMETER :: ICE_RHOS = 1 ! 0. FIXED; 1. PREDICT; 2. MAPPING; AGGREGATES DENSITY INTEGER, PRIVATE, PARAMETER :: ICE_RHOG = 1 ! 0. FIXED; 1. CL1993; 2. RH1985; PREDICT GRAUPEL DENSITY INTEGER, PRIVATE, PARAMETER :: ICE_VENT = 2 ! 0. OFF; 1. HP76; 2. CL94; 3. JW99; VENTILATION EFFECT INTEGER, PRIVATE, PARAMETER :: HAIL_VENT = 1 ! 0. DEFAULT; 1. CHENG ET AL. (2014); VENTILATION EFFECT INTEGER, PRIVATE, PARAMETER :: HWET_MODE = 1 ! 0. OFF; 1. ON; HAILSTONE WET GROWTH INTEGER, PRIVATE, PARAMETER :: LIQ_VTC = 1 ! 0. RAMS/FIXED; 1.,2. CL04; CLOUD DROPS/FLUX FALL SPEED INTEGER, PRIVATE, PARAMETER :: LIQ_VTR = 1 ! 0. RAMS/FIXED; 1.,2. CL04; RAIN DROPS/FLUX FALL SPEED INTEGER, PRIVATE, PARAMETER :: ICE_VTI = 1 ! 0. MORR/FIXED; 1. MH05; PRISTINE ICE FALL SPEED INTEGER, PRIVATE, PARAMETER :: ICE_VTS = 1 ! 0. RAMS/FIXED; 1. MH05/EM16 SYNOPTIC; ! 2. SH09 CRYSTAL-FACE; AGGREGATES FALL SPEED INTEGER, PRIVATE, PARAMETER :: ICE_VTG = 1 ! 0. RAMS/FIXED; 1. MH05; GRAUPEL FALL SPEED INTEGER, PRIVATE, PARAMETER :: ICE_VTH = 1 ! 0. RAMS/FIXED; 1. MH05; HAIL FALL SPEED INTEGER, PRIVATE, PARAMETER :: AFAC_3M = 1 ! 0. FIXED; 1.,2. VARIABLE AFAC PARAMETER INTEGER, PRIVATE, PARAMETER :: AFAR_3M = 1 ! 0. FIXED; 1.,2. VARIABLE AFAR PARAMETER INTEGER, PRIVATE, PARAMETER :: AFAI_3M = 1 ! 0. FIXED; 1.,2. VARIABLE AFAI PARAMETER INTEGER, PRIVATE, PARAMETER :: AFAS_3M = 1 ! 0. FIXED; 1.,2. VARIABLE AFAS PARAMETER INTEGER, PRIVATE, PARAMETER :: AFAG_3M = 1 ! 0. FIXED; 1.,2. VARIABLE AFAG PARAMETER INTEGER, PRIVATE, PARAMETER :: AFAH_3M = 1 ! 0. FIXED; 1.,2. VARIABLE AFAH PARAMETER INTEGER, PRIVATE, PARAMETER :: SAT_ADJ = 0 ! 0. OFF (SUPER-SAT); 1. WITH SATURATION ADJUSTMENT INTEGER, PRIVATE, PARAMETER :: NCCN = 3 ! for aerosol modes INTEGER, PRIVATE, PARAMETER :: NAER = 2 ! for aerosol categories of CCN and IN INTEGER, PRIVATE, PARAMETER :: MAER = 4, NAER1 = 4, NAER2 = 1 INTEGER, PRIVATE, PARAMETER :: NAERT = NAER1+NAER2 ! 5 INTEGER, PRIVATE, PARAMETER :: NTBXA = 25 ! size of tables for the first guess of Rc INTEGER, PRIVATE, DIMENSION(NAER) :: IBAER INTEGER, PRIVATE, DIMENSION(NAER) :: NAERN(1:NAER)=(/NAER1,NAER2/) DOUBLE PRECISION, PRIVATE, SAVE, DIMENSION(NTBXA) :: TBLRC ! table of curresponding R-cut-off DOUBLE PRECISION, PRIVATE, SAVE, DIMENSION(NTBXA,NAER) :: TBLXF !----------------- FOR MICROPHYSICAL PROCESSES ------------------------- REAL, PRIVATE, PARAMETER :: DTMIN = 0.01 ! [sec] mininum dt for cond./evap. calculation REAL, PRIVATE, PARAMETER :: DT20S = 20. ! [sec] dt for coalscence and accretion calculation REAL, PRIVATE, PARAMETER :: PI = 3.1415926535897932384626434 REAL, PRIVATE, PARAMETER :: SQRTPI = 9.189385332046727417803297E-1 REAL, PRIVATE, PARAMETER :: SQRT2 = 1.4142135623730950488016887 REAL, PRIVATE, PARAMETER :: THRD = 1./3., C4PI3 = 4.*PI/3. REAL, PRIVATE, PARAMETER :: CP = 1.00546E3, TK0C = 2.7315E2 REAL, PRIVATE, PARAMETER :: R = 2.87058E2, RV = 4.61495E2 REAL, PRIVATE, PARAMETER :: CPI = 2.093E3, CPW = 4.218E3 ! SPECIFIC HEAT REAL, PRIVATE, PARAMETER :: CMW = 1.8015E-2 REAL, PRIVATE, PARAMETER :: RHOSU = 8.5E4/(2.8715E2*2.7315E2) REAL, PRIVATE, PARAMETER :: RHOW = 9.97E2, RHOG1 = 4.E2 ! BULK DENSITY REAL, PRIVATE, PARAMETER :: RHOI0 = 9.1E2, iRHOI0 = 1./RHOI0 ! SOLID ICE DENSITY REAL, PRIVATE, PARAMETER :: RHOS0 = 1.E2, RHOI1 = 5.E2 REAL, PRIVATE, PARAMETER :: RHOG0 = (0.078+0.184*6.-0.015*36.)*1000 REAL, PRIVATE, PARAMETER :: RHOH = 9.E2, iRHOH = 1./RHOH REAL, PRIVATE, PARAMETER :: RHOIMIN = 5.E1, RHOIMAX = RHOI0 REAL, PRIVATE, PARAMETER :: C4PI3W = 4.*PI*RHOW/3. REAL, PRIVATE, PARAMETER :: iAMI0 = 6./(PI*RHOI0) REAL, PRIVATE, PARAMETER :: G = 9.806, iAPW = 4./(PI*RHOW) REAL, PRIVATE, PARAMETER :: AAW = PI/4., V2M3 = 6./PI REAL, PRIVATE, PARAMETER :: BMW = 3., AMW = PI*RHOW/6. ! MASS-DIAMETER PARAMETERS REAL, PRIVATE, PARAMETER :: BMI0 = 3., AMI0 = PI*RHOI0/6. ! REAL, PRIVATE, PARAMETER :: BMS = 3., AMS0 = PI*RHOS0/6. ! REAL, PRIVATE, PARAMETER :: BMG = 3., AMG0 = PI*RHOG1/6. ! REAL, PRIVATE, PARAMETER :: BMH = 3., AMH = PI*RHOH/6. ! REAL, PRIVATE, PARAMETER :: iAMW = 1./AMW, iAMH = 1./AMH ! REAL, PRIVATE, PARAMETER :: AIMM = 6.6E-1, BIMM = 1.E2 ! BIGG 1953 PARAMETERS FOR IMMERSION FREEZING REAL, PRIVATE, PARAMETER :: AVSG = 8.6E-1, BVSG = 2.8E-1 ! VENTILATION COEFFICIENTS REAL, PRIVATE, PARAMETER :: AVRH = 7.8E-1, BVRH = 3.08E-1 ! REAL, PRIVATE, PARAMETER :: AVIS = 1., BVIS = 1.4E-1 ! REAL, PRIVATE, PARAMETER :: VENC1 = 3.09E-2, VENC2 = 1.447E-1 ! VENTILATION COEFFICIENTS REAL, PRIVATE, PARAMETER :: VENP1 = 1.05E-2, VENP2 = 2.28E-2 ! FOR PRISITNE ICE CRYSTAL REAL, PRIVATE, PARAMETER :: VENH1 = 0.22385, VENH2 = 0.00101 ! FOR HAILSTONE (Cheng et al., 2014) REAL, PRIVATE, PARAMETER :: NSMALL = 1.E-2, NSMAL1 = 1.E2 ! MOMENTS LIMITS REAL, PRIVATE, PARAMETER :: RLIMIT = 1.E-32, SLIMIT = 1.E-2 ! REAL, PRIVATE, PARAMETER :: QSMAL1 = 1.E-9, QSMALL = 1.E-14 ! REAL, PRIVATE, PARAMETER :: QLIMIT = 1.E-6, RSMALL = 1.E-20 ! REAL, PRIVATE, PARAMETER :: ASMALL = 1.E-12, ISMALL = 1.E-17 ! REAL, PRIVATE, PARAMETER :: BOLTZ = 1.38E-23, MLIMIT = 1.E-2 ! REAL, PRIVATE, PARAMETER :: VTZ0 = 5.83, VTC0 = 0.6 ! SURFACE ROUGHNESS PARAMETERS REAL, PRIVATE, PARAMETER :: VTC1 = 0.151931, VTC2 = VTZ0**2./4. ! VTC1=4./(VTZ0**2.*VTC0**5.E-1) REAL, PRIVATE, PARAMETER :: VTA0 = 1.7E-3, VTB0 = 0. ! VTB0 = 0.8 REAL, PRIVATE, PARAMETER :: AVC0 = 3.E7, BVC0 = 2. ! Morrison scheme REAL, PRIVATE, PARAMETER :: AVR0 = 841.997, BVR0 = 0.8 ! Morrison scheme REAL, PRIVATE, PARAMETER :: AVI0 = 700., BVI0 = 1. ! Morrison scheme REAL, PRIVATE, PARAMETER :: AVS0 = 11.72, BVS0 = 0.41 ! Locatelli and Hobbs (1974) REAL, PRIVATE, PARAMETER :: AVG0 = 19.3, BVG0 = 0.37 ! Ferrier (1994) REAL, PRIVATE, PARAMETER :: AVH0 = 206.89, BVH0 = 0.6384 ! Ferrier (1994) REAL, PRIVATE, PARAMETER :: VTCMAX = 1., VTSMAX = 10. ! UPPER LIMIT FALL SPEED CHECK REAL, PRIVATE, PARAMETER :: VTIMAX = 10., VTRMAX = 15. ! REAL, PRIVATE, PARAMETER :: VTGMAX = 20., VTHMAX = 25. ! REAL, PRIVATE, PARAMETER :: AFAC0 = 0., AFAR0 = 0. ! INITIAL SPECTRAL INDEX REAL, PRIVATE, PARAMETER :: AFAI0 = 3., AFAS0 = 0. ! REAL, PRIVATE, PARAMETER :: AFAG0 = 0., AFAH0 = 0. ! REAL, PRIVATE, PARAMETER :: SASMAX = 1., SASMIN = 1.E-3 ! ASPECT RATIO CHECK FOR AGGREGATES SHAPE REAL, PRIVATE, PARAMETER :: KCCMIN = 0.223, KCCMAX = 0.999 ! MOMENT RATIO CHECK REAL, PRIVATE, PARAMETER :: KCRMIN = 0.223, KCRMAX = 0.999 ! REAL, PRIVATE, PARAMETER :: KCIMIN = 0.556, KCIMAX = 0.999 ! KC = 0.999, AFAMAX = 2996 REAL, PRIVATE, PARAMETER :: KCSMIN = 0.223, KCSMAX = 0.999 ! KC = 0.995, AFAMAX = 596 REAL, PRIVATE, PARAMETER :: KCGMIN = 0.223, KCGMAX = 0.999 ! KC = 0.99, AFAMAX = 296 REAL, PRIVATE, PARAMETER :: KCHMIN = 0.223, KCHMAX = 0.999 ! KC = 0.95, AFAMAX = 56 REAL, PRIVATE, PARAMETER :: AFAMAX = 3.E4, AFAMIN = 0. ! SPECTRAL INDEX CHECK REAL, PRIVATE, PARAMETER :: AFU = 3.125E-1, BFU = 1.0552E-3 ! CL04 FALL SPEED ADJUSTMENT REAL, PRIVATE, PARAMETER :: CFU = -2.4023 REAL, PRIVATE, PARAMETER :: DCMIN = 1.E-7, DCMAX = 1.E-4 ! DIAMETER CHECK REAL, PRIVATE, PARAMETER :: DIMIN = 1.E-6, DIMAX = 5.E-3 ! REAL, PRIVATE, PARAMETER :: DRMIN = 3.E-5, DRMAX = 6.E-3 ! REAL, PRIVATE, PARAMETER :: DSMIN = 2.E-5, DSMAX = 1.E-2 ! REAL, PRIVATE, PARAMETER :: DGMIN = 5.E-5, DGMAX = 2.E-2 ! REAL, PRIVATE, PARAMETER :: DHMIN = 1.E-3, DHMAX = 4.E-2 ! REAL, PRIVATE, PARAMETER :: RCMIN = DCMIN/2., RCMAX = DCMAX/2. ! RADIUS CHECK REAL, PRIVATE, PARAMETER :: RRMIN = DRMIN/2., RRMAX = DRMAX/2. ! REAL, PRIVATE, PARAMETER :: DI0 = 6.E-6, DCR = 100.E-6 ! ISOMETRIC LIMIT FOR PI; THRESHOLD FOR RAINDROPS REAL, PRIVATE, PARAMETER :: MI0 = AMI0*DI0**BMI0 ! INITIAL PI MASS REAL, PRIVATE, PARAMETER :: SIG1 = -6.2685, SIG2 = -2.7312E-1 ! CL04 CLOUD DROPLETS SPECTRUM REAL, PRIVATE, PARAMETER :: SIG3 = 2.2606E-1, MNR1 = -2.2920 ! REAL, PRIVATE, PARAMETER :: MNR2 = -3.5158E-1, MNR3 = 3.4708E-1 ! REAL, PRIVATE, PARAMETER :: EFC1 = -1.2560, EFC2 = -1.7904E-02 ! CL04 CLOUD DROPLETS EFFECTIVE RADIUS REAL, PRIVATE, PARAMETER :: EFC3 = 8.5536E-01, EFR1 = -9.9216E-02 ! CL04 RAINDROPS EFFECTIVE RADIUS REAL, PRIVATE, PARAMETER :: EFR2 = 2.9490E-02, EFR3 = 9.9238E-01 ! REAL, DIMENSION(0:120) :: ITBLE ! DEPOSITION GROWTH COEFFICIENTS DATA ITBLE /1.000000,0.979490,0.959401,0.939723,0.920450, & 0.899498,0.879023,0.857038,0.833681,0.810961, & 0.783430,0.755092,0.703072,0.537032,0.467735, & 0.524807,0.630957,0.812831,1.096478,1.479108, & 1.905461,2.089296,2.290868,2.398833,2.454709, & 2.426610,2.371374,2.290868,2.137962,1.995262, & 1.862087,1.737801,1.621810,1.513561,1.396368, & 1.288250,1.188502,1.096478,1.000000,0.922571, & 0.851138,0.785236,0.724436,0.668344,0.616595, & 0.575440,0.537032,0.501187,0.467735,0.436516, & 0.407380,0.380189,0.354813,0.331131,0.316228, & 0.301995,0.291743,0.285102,0.281838,0.278612, & 0.275423,0.278612,0.281838,0.285102,0.291743, & 0.298538,0.309030,0.319890,0.331131,0.346737, & 0.367282,0.393550,0.426580,0.457088,0.489779, & 0.524807,0.562341,0.609537,0.660693,0.716143, & 0.785236,0.860994,0.954993,1.047129,1.148154, & 1.258925,1.380384,1.496236,1.603245,1.698244, & 1.778279,1.840772,1.883649,1.905461,1.905461, & 1.883649,1.862087,1.840772,1.798871,1.737801, & 1.698244,1.640590,1.584893,1.548817,1.513561, & 1.475707,1.452112,1.428894,1.412538,1.393157, & 1.377209,1.361445,1.348963,1.336596,1.327394, & 1.318257,1.309182,1.303167,1.294196,1.288250,1.279381/ REAL, DIMENSION(0:167) :: IECC ! JW99 RIMING COLLECTION EFFICIENCY COLUMNS DATA IECC /0.00,0.03,0.29,0.25,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.10,0.42,0.50,0.47,0.21,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.15,0.49,0.56,0.55,0.46,0.10,0.00, & 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.20,0.52,0.61,0.62,0.59,0.49, & 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.25,0.62,0.72,0.75,0.74, & 0.71,0.68,0.57,0.00,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.00,0.30,0.70,0.80,0.84, & 0.85,0.85,0.84,0.83,0.81,0.77,0.69,0.10,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.35,0.75,0.85, & 0.89,0.91,0.92,0.92,0.91,0.90,0.89,0.87,0.85,0.82, & 0.71,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.40,0.80, & 0.89,0.93,0.94,0.95,0.95,0.95,0.95,0.95,0.95,0.94, & 0.93,0.92,0.91,0.88,0.80,0.15,0.00,0.00/ REAL, DIMENSION(0:167) :: IEPC ! JW99 RIMING COLLECTION EFFICIENCY PLATES DATA IEPC /0.00,0.00,0.13,0.41,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.25,0.54,0.56,0.39,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.08,0.57,0.75,0.83,0.85,0.86,0.84, & 0.78,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.10,0.65,0.80,0.87,0.90,0.91, & 0.91,0.91,0.88,0.78,0.00,0.00,0.00,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.12,0.67,0.81,0.88,0.91, & 0.92,0.93,0.93,0.92,0.91,0.89,0.81,0.00,0.00,0.00, & 0.00,0.00,0.00,0.00,0.00,0.00,0.15,0.68,0.82,0.89, & 0.91,0.93,0.94,0.94,0.95,0.95,0.95,0.94,0.91,0.83, & 0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.18,0.69,0.82, & 0.89,0.92,0.94,0.95,0.95,0.96,0.96,0.96,0.96,0.95, & 0.94,0.91,0.80,0.00,0.00,0.00,0.00,0.00,0.22,0.70, & 0.83,0.90,0.93,0.95,0.96,0.97,0.98,0.98,0.98,0.98, & 0.97,0.96,0.94,0.92,0.84,0.00,0.00,0.00/ REAL, DIMENSION(0:44) :: AMS1 ! EM16, MASS-DIAMETER PARAMETERS SYNOPTIC CLOUDS DATA AMS1 /2.728E+01,1.258E+00,6.452E-02,6.259E-03,8.658E-04, & 4.657E+01,1.193E+00,3.790E-02,2.688E-03,3.101E-04, & 2.260E+01,1.000E+00,5.059E-02,4.889E-03,6.826E-04, & 6.176E+00,6.002E-01,5.889E-02,8.892E-03,1.647E-03, & 2.982E+00,1.795E-01,1.261E-02,1.647E-03,3.128E-04, & 2.783E+00,2.179E-01,1.901E-02,2.825E-03,5.699E-04, & 1.630E+00,3.267E-01,6.223E-02,1.493E-02,3.902E-03, & 1.145E+01,3.346E-02,1.178E-04,7.701E-05,2.463E-05, & 6.606E+01,1.799E-03,1.178E-04,7.701E-05,2.463E-05/ REAL, DIMENSION(0:44) :: BMS1 ! SYNOPTIC CLOUDS DATA BMS1 /2.792,2.455,2.085,1.748,1.411,2.846,2.449,2.015,1.618,& 1.221,2.773,2.429,2.053,1.710,1.367,2.642,2.371,2.073,& 1.802,1.530,2.556,2.254,1.923,1.621,1.320,2.549,2.276,& 1.977,1.704,1.431,2.495,2.322,2.133,1.960,1.787,2.686,& 2.064,1.382,1.382,1.382,2.863,1.732,1.382,1.382,1.382/ REAL, DIMENSION(0:44) :: AAS1 ! EM16, AREA-DIAMETER PARAMETERS SYNOPTIC CLOUDS DATA AAS1 /1.782E+00,2.778E-01,4.794E-02,1.241E-02,4.069E-03, & 2.571E+00,2.635E-01,3.180E-02,6.517E-03,1.879E-03, & 1.910E+00,2.522E-01,3.759E-02,8.905E-03,2.787E-03, & 6.238E-01,1.697E-01,4.686E-02,1.662E-02,6.761E-03, & 5.604E-01,9.194E-02,1.726E-02,5.030E-03,1.927E-03, & 2.802E-01,8.345E-02,2.602E-02,1.042E-02,4.854E-03, & 1.269E-01,1.106E-01,8.453E-02,5.913E-02,3.763E-02, & 3.569E-01,2.314E-02,1.664E-03,1.397E-03,8.521E-04, & 6.339E-01,5.656E-03,1.664E-03,1.397E-03,8.521E-04/ REAL, DIMENSION(0:44) :: BAS1 ! SYNOPTIC CLOUDS DATA BAS1 /2.133,1.938,1.725,1.531,1.337,2.170,1.932,1.671,1.432,& 1.194,2.140,1.927,1.693,1.480,1.267,2.027,1.882,1.722,& 1.576,1.431,2.011,1.821,1.612,1.422,1.232,1.941,1.810,& 1.666,1.534,1.403,1.861,1.842,1.821,1.801,1.782,1.960,& 1.669,1.350,1.350,1.350,2.018,1.509,1.350,1.350,1.350/ !----------------- FOR AEROSOL PARAMETER ------------------------------- REAL, SAVE :: DNC0,DNC1,DNC2,DNC3 REAL, SAVE, DIMENSION(NAER) :: SENS,ASH,CMAS,DNAS,AVAN,BETA1,RXMIN REAL, SAVE, DIMENSION(NCCN,NAER) :: ZCCN,CNMOD,CNSTD,WMAS,RFACT CONTAINS !====================================================================== SUBROUTINE NTU_INIT(PHB,PH,P,PB,ALT,QV,QDCN,QTCN,QCCN,QRCN,QNIN, & XLAND,CCNTY,RESTART,IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME, & JMS,JME,KMS,KME,ITS,ITE,JTS,JTE,KTS,KTE) !====================================================================== IMPLICIT NONE LOGICAL, INTENT(IN) :: RESTART INTEGER, INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME, & KMS,KME,ITS,ITE,JTS,JTE,KTS,KTE,CCNTY REAL, INTENT(IN), DIMENSION(IMS:IME,JMS:JME) :: XLAND ! 1:land ; 2:ocean REAL, INTENT(IN), DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: PHB,PH,P,& PB,ALT,QV REAL, INTENT(INOUT), DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: QDCN, & QTCN,QCCN,QRCN,QNIN REAL, DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE) :: DZ,RHO,DZ8W,P_PHY REAL, DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE,NAERT) :: QAERO INTEGER :: I,J,K,NK,ITF,JTF ITF = MIN(ITE,IDE-1) JTF = MIN(JTE,JDE-1) CALL AERO_CONST(CCNTY) DO J = JTS,JTF DO K = KTS,KTE DO I = ITS,ITF DZ(I,K,J) = (PHB(I,K,J)+PH(I,K,J))/G ENDDO ENDDO ENDDO DO J = JTS,JTF DO K = KTS,KTE-1 DO I = ITS,ITF NK = KTE-K DZ8W(I,K,J) = DZ(I,NK+1,J)-DZ(I,NK,J) RHO(I,K,J) = 1./ALT(I,NK,J)*(1.+QV(I,NK,J)) P_PHY(I,K,J) = P(I,NK,J)+PB(I,NK,J) QAERO(I,K,J,1) = QDCN(I,NK,J) QAERO(I,K,J,2) = QTCN(I,NK,J) QAERO(I,K,J,3) = QCCN(I,NK,J) QAERO(I,K,J,4) = QRCN(I,NK,J) QAERO(I,K,J,5) = QNIN(I,NK,J) ENDDO ENDDO ENDDO IF (.NOT.RESTART) THEN CALL INIT_AEROSOL(P_PHY,RHO,DZ8W,XLAND,QAERO,IDS,IDE,JDS,JDE, & KDS,KDE,IMS,IME,JMS,JME,KMS,KME,ITS,ITE,JTS,JTE,KTS,KTE) ENDIF DO J = JTS,JTF DO K = KTS,KTE-1 DO I = ITS,ITF NK = KTE-K QDCN(I,NK,J) = QAERO(I,K,J,1) QTCN(I,NK,J) = QAERO(I,K,J,2) QCCN(I,NK,J) = QAERO(I,K,J,3) QRCN(I,NK,J) = QAERO(I,K,J,4) QNIN(I,NK,J) = QAERO(I,K,J,5) ENDDO ENDDO ENDDO END SUBROUTINE NTU_INIT !====================================================================== ! !====================================================================== SUBROUTINE FIND_RC0(XAFRC,CMODE,CSTDV,WMAS,RC,TBLXA,TBLRC) !====================================================================== IMPLICIT NONE INTEGER, PARAMETER :: ITERMAX = 50 INTEGER :: IM,I !----------------- AEROSOL TRI-MODAL SIZE DISTRIBUTION ---------------- REAL :: RC ! cut-off radius [m] REAL, DIMENSION(NCCN) :: WMAS, &! mass weighting for each mode CMODE, &! modal value (ln R) of each mode [m] CSTDV ! standard distribution (in logarithm) of the mode REAL, PARAMETER :: RC_MAX = 99., RC_MIN = 1.E-9, XTOR = 1.E-5 DOUBLE PRECISION, DIMENSION(NCCN) :: DWMAS,DMODE,DSTDV !----------------- GET THE FIRST GUESS OF RC --------------------------- DOUBLE PRECISION :: DRC,UL,X1MAFRC,Y0,Y,DY,DX,DERF,DLOG,DEXP, & DSQRT,XAFRC ! ratio of dry/total aerosol mass DOUBLE PRECISION, DIMENSION(NTBXA) :: TBLRC,TBLXA ! table of dry aerosol mass fraction DO IM = 1,NCCN DWMAS(IM) = DBLE(WMAS(IM)) DSTDV(IM) = DBLE(CSTDV(IM)) DMODE(IM) = DLOG(DBLE(CMODE(IM)))+DBLE(3.*CSTDV(IM)*CSTDV(IM)) ! for mass, the 3rd momentum DBLE(CMODE (IM))+ ENDDO UL = (1.-DERF(5.D+0/DSQRT(2.D+0)))*0.5 ! after 5 sigma, x1mfrc will < 2.86E-7 X1MAFRC = 1.D+0-XAFRC IF (XAFRC.LT.1.D-10) THEN RC = RC_MIN RETURN ENDIF IF (X1MAFRC.LT.UL) THEN RC = RC_MAX RETURN ENDIF DRC = GUESS_RC(XAFRC,TBLXA,TBLRC) ! cut-off radius [m] Y0 = YEQU(DRC,XAFRC,DWMAS,DMODE,DSTDV) DO I = 1,ITERMAX Y = Y0 DY = DYEQU(DRC,DWMAS,DMODE,DSTDV) IF (DY.LE.1.D-50) THEN PRINT *,'IN FIND_RC0.F DY IS',DY PRINT *,I,DX,Y,DY,DEXP(DRC),DRC,XAFRC,X1MAFRC,DMODE,DSTDV, & CMODE,CSTDV STOP ENDIF DX = -Y/DY IF (DABS(DX).GT.3.D+0) DX = DX/(DABS(DX/0.4)) IF (DABS(DX).LT.XTOR) THEN RC = REAL(DEXP(DRC+DX)) RC = MIN(RC_MAX,MAX(RC,RC_MIN)) RETURN ENDIF Y0 = YEQU(DRC+DX,XAFRC,DWMAS,DMODE,DSTDV) IF (Y*Y0.LT.0.) THEN DX = DX*5.D-1 Y0 = YEQU(DRC+DX,XAFRC,DWMAS,DMODE,DSTDV) ENDIF DRC = DRC+DX ENDDO RC = REAL(DEXP(DRC)) END SUBROUTINE FIND_RC0 !====================================================================== ! !====================================================================== FUNCTION GUESS_RC(XAFRC,TBLXA,TBLRC) !====================================================================== DOUBLE PRECISION :: GUESS_RC,XAFRC,TBLXA(NTBXA),TBLRC(NTBXA) INTEGER :: IBNG,IEND,I IF (XAFRC.GE.TBLXA(NTBXA)) THEN GUESS_RC = TBLRC(NTBXA) RETURN ELSEIF (XAFRC.LT.TBLXA(1)) THEN GUESS_RC = TBLRC(1) RETURN ENDIF IF (XAFRC.LT.TBLXA(NTBXA/4)) THEN IBNG = 1 IEND = NTBXA/4-1 ELSEIF (XAFRC.LT.TBLXA(NTBXA*2/4)) THEN IBNG = NTBXA/4 IEND = NTBXA*2/4-1 ELSEIF (XAFRC.LT.TBLXA(NTBXA*3/4)) THEN IBNG = NTBXA*2/4 IEND = NTBXA*3/4-1 ELSE IBNG = NTBXA*3/4 IEND = NTBXA-1 ENDIF DO I = IEND,IBNG,-1 IF (XAFRC.GE.TBLXA(I)) GOTO 111 ENDDO 111 GUESS_RC = TBLRC(I) END FUNCTION GUESS_RC !====================================================================== ! !====================================================================== FUNCTION YEQU(X,XAFRC,WEGHT,XLBAR,SIGMA) !====================================================================== INTEGER :: I DOUBLE PRECISION :: YEQU,DX,XAFRC,X,DSQRT2,DSQRT DOUBLE PRECISION, DIMENSION(NCCN) :: WEGHT,XLBAR,SIGMA DSQRT2 = DSQRT(2.D+0) YEQU = 0.D+0 DO I = 1,NCCN DX = (X-XLBAR(I))/SIGMA(I)/DSQRT2 YEQU = YEQU+WEGHT(I)*PDF(DX) ENDDO YEQU = YEQU-XAFRC END FUNCTION YEQU !---------------------------------------------------------------------- FUNCTION DYEQU(X,WEGHT,XLBAR,SIGMA) !====================================================================== INTEGER :: I DOUBLE PRECISION :: DYEQU,X,DX,DSQRT2,DSQRT DOUBLE PRECISION, DIMENSION(NCCN) :: WEGHT,XLBAR,SIGMA DSQRT2 = DSQRT(2.D+0) DYEQU = 0.D+0 DO I = 1,NCCN DX = (X-XLBAR(I))/SIGMA(I)/DSQRT2 DYEQU = DYEQU+WEGHT(I)*DPDF(DX)/SIGMA(I)/DSQRT2 ENDDO END FUNCTION DYEQU !====================================================================== FUNCTION PDF(X) DOUBLE PRECISION :: PDF,X,DERF PDF = (1.D+0+DERF(X))*5.D-1 END FUNCTION PDF !---------------------------------------------------------------------- FUNCTION DPDF(X) DOUBLE PRECISION :: DPDF,X,DPI,DACOS,DEXP,DSQRT DPI = DACOS(-1.D+0) DPDF = DEXP(-X*X)/DSQRT(DPI) END FUNCTION DPDF !====================================================================== ! !====================================================================== REAL FUNCTION GAMMA(X) ! IMPLEMETED FROM MORRISON SCHEME !====================================================================== 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 !---------------------------------------------------------------------- ! 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.71618513886549492533811E0,2.47656508055759199108314E1,& -3.79804256470945635097577E2,6.29331155312818442661052E2,& 8.66966202790413211295064E2,-3.14512729688483675254357E4,& -3.61444134186911729807069E4,6.64561438202405440627855E4/ DATA Q /-3.08402300119738975254353E1,3.15350626979604161529144E2,& -1.01515636749021914166146E3,-3.10777167157231109440444E3,& 2.25381184209801510330112E4,4.75584627752788110767815E3, & -1.34659959864969306392456E5,-1.15132259675553483497211E5/ !---------------------------------------------------------------------- ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). !---------------------------------------------------------------------- DATA C /-1.910444077728E-3,8.4171387781295E-4, & -5.952379913043012E-4,7.93650793500350248E-4, & -2.777777777777681622553E-3,8.333333333333333331554247E-2,& 5.7083835261E-3/ !---------------------------------------------------------------------- ! 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+SQRTPI SUM = SUM+(Y-HALF)*LOG(Y) RES = EXP(SUM) ELSE RES = XINF GOTO 900 ENDIF ENDIF !---------------------------------------------------------------------- ! FINAL ADJUSTMENTS AND RETURN !---------------------------------------------------------------------- IF (PARITY) RES = -RES IF (FACT.NE.ONE) RES = FACT/RES 900 GAMMA = RES RETURN END FUNCTION GAMMA !====================================================================== !====================================================================== REAL FUNCTION POLYSVP(T,TYPE) ! IMPLEMETED FROM MORRISON SCHEME !====================================================================== ! 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,T,DT 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 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/ IF (TYPE.EQ.1) THEN ! ICE 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. ENDIF IF (TYPE.EQ.0) THEN ! LIQUID 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. ENDIF END FUNCTION POLYSVP !====================================================================== ! !====================================================================== REAL FUNCTION GAMLN(XX) ! Referred to MY2 scheme !====================================================================== ! Returns value of ln(GAMMA(XX)) for XX>0 (modified from "Numerical Recipes") IMPLICIT NONE REAL, INTENT(IN) :: XX INTEGER :: J DOUBLE PRECISION :: ser,stp,TMP,X,y,cof(6) SAVE cof,stp DATA cof,stp /76.18009172947146d0,-86.50532032941677d0, & 24.01409824083091d0,-1.231739572450155d0, & .1208650973866179d-2,-.5395239384953d-5, & 2.5066282746310005d0/ X = DBLE(XX) y = X TMP = X+5.5D0 TMP = (X+0.5D0)*LOG(TMP)-TMP ser = 1.000000000190015d0 DO J = 1,6 !original y = y+1.D0 ser = ser+cof(J)/y ENDDO #ifdef DOUBLE_PRECISION GAMLN = TMP+LOG(stp*ser/X) #else GAMLN = SNGL(TMP+LOG(stp*ser/X)) #endif END FUNCTION GAMLN !====================================================================== ! !====================================================================== REAL FUNCTION GAMIN(P,XMAX) !====================================================================== ! The incomplete gamma function below the limit of P = AFA+1 and XMAX = DI0*LAMDA REAL :: P,XMAX GAMIN = GAMMP(P,XMAX)*EXP(GAMLN(P)) END FUNCTION GAMIN !======================================================================= REAL FUNCTION GAMMP(A,X) !======================================================================= ! The fraction of distribution below the limit of A = AFA+1 and X = Db*LAMDA IMPLICIT NONE REAL :: A,X,GAMMCF,GAMSER,GLN IF (X.LT.0..OR.A.LE.0.) & CALL wrf_error_fatal('warning : bad arguments in gammq') IF (X.LT.A+1.) THEN CALL GSER(GAMSER,A,X,GLN) GAMMP = GAMSER ELSE CALL CFG(GAMMCF,A,X,GLN) GAMMP = 1.-GAMMCF ENDIF RETURN END FUNCTION GAMMP !====================================================================== SUBROUTINE GSER(GAMSER,A,X,GLN) ! Referred to MY2 scheme !====================================================================== ! USES GAMLN, Returns the incomplete gamma function P(A,X) evaluated by its series ! representation as GAMSER. Also returns GAMMA(A) as GLN. IMPLICIT NONE INTEGER :: N REAL :: A,GAMSER,GLN,X,AP,de1,summ INTEGER, PARAMETER :: ITMAX = 500 REAL, PARAMETER :: EPS = 3.E-7 GLN = GAMLN(A) IF (X.LE.0.) THEN IF (X.LT.0.) CALL wrf_error_fatal ( 'WARNING: X <0 in GSER' ) GAMSER = 0. RETURN ENDIF AP = A summ = 1./A de1 = summ DO N = 1,ITMAX AP = AP+1. de1 = de1*X/AP summ = summ+de1 IF (ABS(de1).LT.ABS(summ)*EPS) GOTO 777 ENDDO CALL wrf_error_fatal ('Warning : ITMAX too small in GSER') 777 GAMSER = summ*EXP(-X+A*LOG(X)-GLN) RETURN END SUBROUTINE GSER !======================================================================= SUBROUTINE CFG(GAMMCF,A,X,GLN) ! Referred to MY2 scheme !======================================================================= ! USES GAMLN, Returns the incomplete gamma function (Q(A,X) evaluated by tis continued fraction ! representation as GAMMCF. Also returns ln(GAMMA(A)) as GLN. ITMAX is the maximum allowed number of iterations; ! EPS is the relative accuracy; FPMIN is a number near the smallest representable floating-point number. IMPLICIT NONE INTEGER :: I REAL :: A,GAMMCF,GLN,X,AN,b,c,d,de1,h INTEGER, PARAMETER :: ITMAX = 500 REAL, PARAMETER :: EPS = 3.E-7 REAL, PARAMETER :: fpmin = 1.E-30 GLN = GAMLN(A) b = X+1.-A c = 1./fpmin d = 1./b h = d DO I = 1,ITMAX AN = -I*(I-A) b = b+2. d = AN*d+b IF (ABS(d).LT.fpmin) d = fpmin c = b+AN/c IF (ABS(c).LT.fpmin) c = fpmin d = 1./d de1 = d*c h = h*de1 IF (ABS(de1-1.).LT.EPS) GOTO 888 ENDDO CALL wrf_error_fatal ('Warning : ITMAX too small in gcf') 888 GAMMCF = EXP(-X+A*LOG(X)-GLN)*h RETURN END SUBROUTINE CFG !======================================================================= ! !====================================================================== SUBROUTINE SOLVE_AFAC(TK1D,QC1D,NC1D,LAMC,MVDC,AFAC) !====================================================================== IMPLICIT NONE REAL :: TK1D,QC1D,NC1D,LAMC,GC1,MVDC,AFAC,LAMCMIN,LAMCMAX,C3M1D, & SIGC,MNRC,MVRC,EFRC,KDX,LTK,LQC IF (NC1D.LT.NSMALL) THEN LTK = LOG(TK1D) LQC = -1.*LOG(QC1D) NC1D = EXP(DNC0+DNC1*LTK+DNC2*LTK**2.+DNC3*LTK**3.-0.25*LQC) ENDIF IF (QC1D.GE.QSMAL1.AND.NC1D.LT.NSMAL1) THEN LTK = LOG(TK1D) LQC = -1.*LOG(QC1D) NC1D = EXP(DNC0+DNC1*LTK+DNC2*LTK**2.+DNC3*LTK**3.-0.25*LQC) ENDIF C3M1D = QC1D*V2M3/RHOW IF (AFAC_3M.EQ.0) THEN AFAC = AFAC0 ELSEIF (AFAC_3M.EQ.1) THEN MVRC = (QC1D/NC1D/C4PI3W)**THRD MVRC = MIN(MAX(MVRC,RCMIN),RCMAX) EFRC = EXP(EFC1+EFC2*LOG(NC1D)+EFC3*LOG(MVRC)) ! KDX = MIN(MAX(KCCMIN,C3M1D/(8.*EFRC**3.*NC1D)),KCCMAX) KDX = MAX(KCCMIN,MIN(KCCMAX,(MVRC/EFRC)**3.)) AFAC = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAC = MIN(MAX(AFAC,AFAMIN),AFAMAX) ELSEIF (AFAC_3M.EQ.2) THEN SIGC = EXP(SIG1+SIG2*LOG(NC1D)+SIG3*LOG(QC1D)) MNRC = EXP(MNR1+MNR2*LOG(NC1D)+MNR3*LOG(QC1D)) AFAC = MIN(MAX(SIGC/MNRC,AFAMIN),AFAMAX) ENDIF GC1 = GAMLN(AFAC+1.) LAMC = (EXP(GAMLN(AFAC+4.)-GC1)*NC1D/C3M1D)**THRD LAMCMIN = (EXP(GAMLN(AFAC+4.)-GC1))**THRD/DCMAX LAMCMAX = (EXP(GAMLN(AFAC+4.)-GC1))**THRD/DCMIN IF (LAMC.LT.LAMCMIN) THEN LAMC = LAMCMIN NC1D = C3M1D*EXP(GAMLN(AFAC+1.)-GAMLN(AFAC+4.)+3.*LOG(LAMC)) ELSEIF (LAMC.GT.LAMCMAX) THEN LAMC = LAMCMAX NC1D = C3M1D*EXP(GAMLN(AFAC+1.)-GAMLN(AFAC+4.)+3.*LOG(LAMC)) ENDIF MVDC = (EXP(GAMLN(AFAC+4.)-GAMLN(AFAC+1.)))**THRD/LAMC END SUBROUTINE SOLVE_AFAC !====================================================================== ! !====================================================================== SUBROUTINE SOLVE_AFAR(TK1D,QR1D,NR1D,LAMR,MVDR,AFAR) !====================================================================== IMPLICIT NONE REAL :: TK1D,QR1D,NR1D,LAMR,GR1,MVDR,LAMRMIN,LAMRMAX,R3M1D,AFAR, & BDR,MVRR,EFRR,KDX,LTK,LQR IF (NR1D.LT.NSMALL) THEN LTK = LOG(TK1D) LQR = -1.*LOG(QR1D) NR1D = EXP(-5793.7852+3191.1171*LTK-582.73279*LTK**2.+ & 35.346854*LTK**3.-0.25*LQR) ENDIF R3M1D = QR1D*V2M3/RHOW IF (AFAR_3M.EQ.0) THEN AFAR = AFAR0 ELSEIF (AFAR_3M.EQ.1) THEN MVRR = MIN(RRMAX,MAX(RRMIN,(QR1D/NR1D/C4PI3W)**THRD)) EFRR = EXP(EFR1+EFR2*LOG(NR1D)+EFR3*LOG(MVRR)) ! KDX = MIN(MAX(KCRMIN,R3M1D/(8.*EFRR**3.*NR1D)),KCRMAX) KDX = MAX(KCRMIN,MIN(KCRMAX,(MVRR/EFRR)**3.)) AFAR = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAR = MIN(MAX(AFAR,AFAMIN),AFAMAX) ELSEIF (AFAR_3M.EQ.2) THEN BDR = MIN(MAX((R3M1D/NR1D)**THRD,DRMIN),DRMAX) AFAR = MAX(AFAMIN,19.*TANH(0.6*(1.E3*BDR-1.8))+17.) ENDIF GR1 = GAMLN(AFAR+1.) LAMR = (EXP(GAMLN(AFAR+4.)-GR1)*NR1D/R3M1D)**THRD LAMRMIN = (EXP(GAMLN(AFAR+4.)-GR1))**THRD/DRMAX LAMRMAX = (EXP(GAMLN(AFAR+4.)-GR1))**THRD/DRMIN IF (LAMR.LT.LAMRMIN) THEN LAMR = LAMRMIN NR1D = R3M1D*EXP(GAMLN(AFAR+1.)-GAMLN(AFAR+4.)+3.*LOG(LAMR)) ELSEIF (LAMR.GT.LAMRMAX) THEN LAMR = LAMRMAX NR1D = R3M1D*EXP(GAMLN(AFAR+1.)-GAMLN(AFAR+4.)+3.*LOG(LAMR)) ENDIF MVDR = (EXP(GAMLN(AFAR+4.)-GAMLN(AFAR+1.)))**THRD/LAMR END SUBROUTINE SOLVE_AFAR !====================================================================== !====================================================================== SUBROUTINE SOLVE_AFAI(TK1D,P1D,RHO,QV1D,QI1D,NI1D,VI1D,FI1D, & I2M1D,I3M1D,ADAGR,ZETA,LAMI,AFAI,MVDI,RHOI,AMI,BMI, & AVI,BVI,BEST) !====================================================================== IMPLICIT NONE INTEGER :: HID REAL :: TK1D,P1D,RHO,QV1D,QI1D,NI1D,VI1D,FI1D,I2M1D,I3M1D,ADAGR, & ZETA,RHOI,LAMI,AFAI,AMI,BMI,AVI,BVI,MVDI,GI1,TC1D,ESW, & ESI,QVSI,INHGR,I3M0,BEST0,LAMIMIN,LAMIMAX,IPF,IPG,KDX, & ZETA2,ZETA4,KINV,AAI,BAI,IBA1,C1X2,VTA1,VTB1,IPH2,IPG2, & AVIA0,BEST,BDI,FDI,LTK,LQI,MDI IF (ICE_RHOI.EQ.0) THEN RHOI = RHOI0 VI1D = 0. ELSEIF (ICE_RHOI.EQ.1) THEN IF (VI1D.LT.ISMALL) THEN TC1D = TK1D-TK0C HID = MAX(MIN(NINT(ABS(TC1D)/0.25),120),0) ESW = MIN(0.99*P1D,POLYSVP(TK1D,0)) ESI = MIN(0.99*P1D,POLYSVP(TK1D,1)) IF (ESI.GT.ESW) ESI = ESW QVSI = 0.622*ESI/(P1D-ESI) INHGR = ITBLE(HID) RHOI = RHOI0*EXP(-3.*MAX((QV1D-QVSI)-5.E-5,0.)/INHGR) VI1D = QI1D/RHOI ENDIF RHOI = QI1D/VI1D IF (RHOI.LT.RHOIMIN) THEN RHOI = RHOIMIN ELSEIF (RHOI.GT.RHOIMAX) THEN RHOI = RHOIMAX ENDIF VI1D = QI1D/RHOI ELSEIF (ICE_RHOI.EQ.2) THEN RHOI = RHOI1 VI1D = 0. ENDIF IF (NI1D.LT.NSMALL) THEN LTK = LOG(TK1D) LQI = -1.*LOG(QI1D) MDI = EXP(-3.2653646+2.0539073*LTK-0.25*LQI)/1.E3 NI1D = 1.E9*QI1D*V2M3/RHOI/MDI**3. ENDIF IF (QI1D.GE.QSMAL1.AND.NI1D.LT.NSMAL1) THEN LTK = LOG(TK1D) LQI = -1.*LOG(QI1D) MDI = EXP(-3.2653646+2.0539073*LTK-0.25*LQI)/1.E3 NI1D = 1.E9*QI1D*V2M3/RHOI/MDI**3. ENDIF IF (I2M1D.GE.ASMALL.AND.I3M1D.GE.ISMALL) THEN KDX = MAX(0.,I2M1D**3./(I3M1D**2.*NI1D)) IF (KDX.GT.KCIMAX) THEN KDX = KCIMIN I2M1D = (KDX*NI1D*I3M1D**2.)**THRD ELSEIF (KDX.LT.KCIMIN) THEN KDX = KCIMIN I2M1D = (KDX*NI1D*I3M1D**2.)**THRD ENDIF AFAI = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAI = MIN(MAX(AFAI,AFAMIN),AFAMAX) ! LAMI = SQRT(NI1D*(AFAI+2.)*(AFAI+1.)/I2M1D) ! LAMI = (NI1D*(AFAI+1.)*(AFAI+2.)*(AFAI+3.)/I3M1D)**THRD LAMI = (AFAI+3.)*I2M1D/I3M1D LAMIMIN = (EXP(GAMLN(AFAI+4.)-GAMLN(AFAI+1.)))**THRD/DIMAX LAMIMAX = (EXP(GAMLN(AFAI+4.)-GAMLN(AFAI+1.)))**THRD/DIMIN IF (LAMI.LT.LAMIMIN) THEN LAMI = LAMIMIN NI1D = I3M1D*EXP(GAMLN(AFAI+1.)-GAMLN(AFAI+4.)+3.*LOG(LAMI)) ELSEIF (LAMI.GT.LAMIMAX) THEN LAMI = LAMIMAX NI1D = I3M1D*EXP(GAMLN(AFAI+1.)-GAMLN(AFAI+4.)+3.*LOG(LAMI)) ENDIF ELSEIF (I2M1D.LT.ASMALL.AND.I3M1D.GE.ISMALL) THEN IF (AFAI_3M.EQ.0) THEN AFAI = AFAI0 I2M1D = 0. ELSEIF (AFAI_3M.EQ.1) THEN AFAI = AFAI0 I2M1D = (KCIMIN*NI1D*I3M1D**2.)**THRD ELSEIF (AFAI_3M.EQ.2) THEN BDI = MIN(MAX((I3M1D/NI1D)**THRD*1.E3,DIMIN*1.E3),DIMAX*1.E3) FDI = 0.074015986+0.79866676*BDI-0.0094468892*LOG(NI1D)+ & 0.38235092*BDI**2.+0.00029811542*LOG(NI1D)**2.+ & 0.019052614*BDI*LOG(NI1D) KDX = MAX(KCIMIN,MIN(KCIMAX,(BDI/FDI)**3.)) AFAI = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAI = MIN(MAX(AFAI,AFAMIN),AFAMAX) ! AFAI = MAX(AFAMIN,12.*TANH(0.7*(BDI-1.7))+11.) I2M1D = 0. ENDIF LAMI = (EXP(GAMLN(AFAI+4.)-GAMLN(AFAI+1.))*NI1D/I3M1D)**THRD LAMIMIN = (EXP(GAMLN(AFAI+4.)-GAMLN(AFAI+1.)))**THRD/DIMAX LAMIMAX = (EXP(GAMLN(AFAI+4.)-GAMLN(AFAI+1.)))**THRD/DIMIN IF (LAMI.LT.LAMIMIN) THEN LAMI = LAMIMIN NI1D = I3M1D*EXP(GAMLN(AFAI+1.)-GAMLN(AFAI+4.)+3.*LOG(LAMI)) ELSEIF (LAMI.GT.LAMIMAX) THEN LAMI = LAMIMAX NI1D = I3M1D*EXP(GAMLN(AFAI+1.)-GAMLN(AFAI+4.)+3.*LOG(LAMI)) ENDIF ELSE IF (AFAI_3M.EQ.0) THEN AFAI = AFAI0 I2M1D = 0. ELSE BDI = (QI1D*V2M3/NI1D/RHOI)**THRD*1.E3 BDI = MAX(MIN(DIMAX*1.E3,BDI),DIMIN*1.E3) FDI = 7.4015986E-2+7.9866676E-1*BDI-9.4468892E-3*LOG( & NI1D)+3.8235092E-1*BDI**2.+2.9811542E-4*LOG(NI1D)** & 2.+1.9052614E-2*BDI*LOG(NI1D) KDX = MAX(KCIMIN,MIN(KCIMAX,(BDI/FDI)**3.)) AFAI = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAI = MIN(MAX(AFAI,AFAMIN),AFAMAX) ! AFAI = MAX(AFAMIN,12.*TANH(0.7*(BDI-1.7))+11.) IF (AFAI_3M.EQ.1) THEN KDX = (AFAI**2.+3.*AFAI+2.)/(AFAI**2.+6.*AFAI+9.) I2M1D = (KDX*NI1D*(QI1D*V2M3/RHOI)**2.)**THRD ELSEIF (AFAI_3M.EQ.2) THEN I2M1D = 0. ENDIF ENDIF GI1 = GAMLN(AFAI+1.) LAMI = (EXP(GAMLN(AFAI+4.)-GI1)*PI*RHOI*NI1D/QI1D/6.)**THRD LAMIMIN = (EXP(GAMLN(AFAI+4.)-GAMLN(AFAI+1.)))**THRD/DIMAX LAMIMAX = (EXP(GAMLN(AFAI+4.)-GAMLN(AFAI+1.)))**THRD/DIMIN IF (LAMI.LT.LAMIMIN) THEN LAMI = LAMIMIN NI1D = QI1D*V2M3/RHOI*EXP(GAMLN(AFAI+1.)-GAMLN(AFAI+4.)+ & 3.*LOG(LAMI)) ELSEIF (LAMI.GT.LAMIMAX) THEN LAMI = LAMIMAX NI1D = QI1D*V2M3/RHOI*EXP(GAMLN(AFAI+1.)-GAMLN(AFAI+4.)+ & 3.*LOG(LAMI)) ENDIF ENDIF MVDI = (EXP(GAMLN(AFAI+4.)-GAMLN(AFAI+1.)))**THRD/LAMI IF (I3M1D.GE.ISMALL.AND.FI1D.GE.ISMALL) THEN I3M0 = NI1D*DI0**3. IF (MVDI.GT.(DI0+1.E-7)) THEN IF (ICE_SHAPE.EQ.0) THEN ZETA = 0.; ADAGR = 1. FI1D = I3M1D ELSEIF (ICE_SHAPE.EQ.1) THEN ZETA = LOG(FI1D/I3M1D)/LOG(I3M1D/I3M0) IF (ZETA.GT.0.4) THEN ZETA = 0.4 FI1D = (I3M1D/I3M0)**ZETA*I3M1D ELSEIF (ZETA.LT.(-0.4)) THEN ZETA = -0.4 FI1D = (I3M1D/I3M0)**ZETA*I3M1D ENDIF ADAGR = (1.+2.*ZETA)/(1.-ZETA) ENDIF ELSE FI1D = I3M1D ZETA = 0.; ADAGR = 1. ENDIF ELSE I3M1D = QI1D*V2M3/RHOI IF (ICE_SHAPE.EQ.0) THEN FI1D = I3M1D ZETA = 0.; ADAGR = 1. ELSEIF (ICE_SHAPE.EQ.1) THEN TC1D = TK1D-TK0C HID = MAX(MIN(NINT(ABS(TC1D)/0.25),120),0) ADAGR = MAX(MIN(ITBLE(HID),2.),0.5)**THRD ZETA = (ADAGR-1.)/(ADAGR+2.) I3M0 = NI1D*DI0**3. FI1D = (I3M1D/I3M0)**ZETA*I3M1D ENDIF ENDIF IF ((ADAGR-1.).GE.SLIMIT) THEN AMI = PI*RHOI*DI0**(2.-2./ADAGR)/6. BMI = 2./ADAGR+1. ELSEIF ((1.-ADAGR).GE.SLIMIT) THEN AMI = PI*RHOI*DI0**(1.-ADAGR)/6. BMI = ADAGR+2. ELSEIF (ABS(ADAGR-1.).LT.SLIMIT) THEN AMI = PI*RHOI/6. BMI = BMI0 ENDIF IF (ICE_VTI.EQ.0) THEN AVI = AVI0 BVI = BVI0 ELSEIF (ICE_VTI.EQ.1) THEN KINV = (1.72E-5*(393./(TK1D+120.))*(TK1D/TK0C)**1.5)/RHO BEST0 = 2.*G*NI1D/(KINV**2.) IF ((ADAGR-1.).GE.SLIMIT) THEN AAI = PI/4./DI0**ZETA BAI = 3.*(ADAGR+1.)/(ADAGR+2.) IPH2 = 2.*3.*ADAGR/(ADAGR+2.) ZETA4 = 4.*(ADAGR-1.)/(ADAGR+2.) IBA1 = BMI+IPH2-BAI BEST = BEST0*AMI*EXP(GAMLN(IBA1+AFAI+1.)-GAMLN(AFAI+1.)- & IBA1*LOG(LAMI))/(AAI*DI0**ZETA4) C1X2 = VTC1*BEST**5.E-1 VTB1 = C1X2/(1.+C1X2)**5.E-1/((1.+C1X2)**5.E-1-1.)/2.- & VTA0*VTB0*BEST**VTB0/VTC2/(SQRT(1.+C1X2)-1.)**2. VTA1 = MAX((VTC2*((1.+C1X2)**5.E-1-1.)**2.-VTA0*BEST** & VTB0)/BEST**VTB1,0.) AVIA0 = VTA1*KINV**(1.-2.*VTB1) AVI = AVIA0*(2.*AMI*G/(AAI*DI0**ZETA4))**VTB1 BVI = VTB1*(BMI+IPH2-BAI)-1. ELSEIF ((1.-ADAGR).GE.SLIMIT) THEN ZETA2 = 2.*(ADAGR-1.)/(ADAGR+2.) AAI = PI*DI0**ZETA2/4. BAI = 2.*3./(ADAGR+2.) IPG2 = 2.*3./(ADAGR+2.) IBA1 = BMI+IPG2-BAI BEST = BEST0*AMI*DI0**ZETA2*EXP(GAMLN(IBA1+AFAI+1.)- & GAMLN(AFAI+1.)-IBA1*LOG(LAMI))/AAI C1X2 = VTC1*BEST**5.E-1 VTB1 = C1X2/(1.+C1X2)**5.E-1/((1.+C1X2)**5.E-1-1.)/2.- & VTA0*VTB0*BEST**VTB0/VTC2/(SQRT(1.+C1X2)-1.)**2. VTA1 = MAX((VTC2*((1.+C1X2)**5.E-1-1.)**2.-VTA0*BEST** & VTB0)/BEST**VTB1,0.) AVIA0 = VTA1*KINV**(1.-2.*VTB1) AVI = AVIA0*(2.*AMI*G*DI0**ZETA2/AAI)**VTB1 BVI = VTB1*(BMI+IPG2-BAI)-1. ELSEIF (ABS(ADAGR-1.).LT.SLIMIT) THEN AAI = PI/4. BAI = 2. IBA1 = BMI BEST = BEST0*AMI*EXP(GAMLN(IBA1+AFAI+1.)-GAMLN(AFAI+1.)- & IBA1*LOG(LAMI))/AAI C1X2 = VTC1*BEST**5.E-1 VTB1 = C1X2/(1.+C1X2)**5.E-1/((1.+C1X2)**5.E-1-1.)/2.- & VTA0*VTB0*BEST**VTB0/VTC2/(SQRT(1.+C1X2)-1.)**2. VTA1 = MAX((VTC2*((1.+C1X2)**5.E-1-1.)**2.-VTA0*BEST** & VTB0)/BEST**VTB1,0.) AVIA0 = VTA1*KINV**(1.-2.*VTB1) AVI = AVIA0*(2.*AMI*G/AAI)**VTB1 BVI = VTB1*BMI-1. ENDIF ENDIF END SUBROUTINE SOLVE_AFAI !====================================================================== ! !====================================================================== SUBROUTINE SOLVE_AFAS(TK1D,RHO,QS1D,QC1D,NS1D,VS1D,FS1D,S2M1D, & AFAS,LAMS,MVDS,RHOS,SASPR,AMS,AVS,BVS) !====================================================================== IMPLICIT NONE INTEGER :: TBIN,DBIN REAL :: TK1D,RHO,QS1D,QC1D,NS1D,VS1D,FS1D,S2M1D,AFAS,RHOS,SASPR, & LAMS,AMS,AVS,BVS,MVDS,LAMSMIN,LAMSMAX,KINV,C1X2,VTA1, & VTB1,BEST,S3M1D,KDX,TC1D,SBA1,AMS2,AAS2,AMS3,BMS3,AAS3, & BAS3,LAMD,BDS,FDS,LTK,LQS,MDS IF (ICE_RHOS.EQ.0) THEN RHOS = RHOS0 VS1D = 0. ELSEIF (ICE_RHOS.EQ.1) THEN RHOS = QS1D/(VS1D+ISMALL) IF (RHOS.GT.RHOIMAX) THEN RHOS = RHOIMAX ELSEIF (RHOS.LT.RHOIMIN) THEN RHOS = RHOIMIN ENDIF VS1D = QS1D/RHOS ELSEIF (ICE_RHOS.EQ.2) THEN LTK = LOG(TK1D) LQS = -1.*LOG(QS1D) IF (TK1D.LT.TK0C) THEN RHOS = 15740.702-6098.0087*LTK+503.33089*LQS+594.29913* & LTK**2.+1.9033961*LQS**2.-94.950429*LTK*LQS ELSE RHOS = EXP(-64808.666+23113.508*LTK-36.46632*LQS-2060.6024*& LTK**2.-0.005729458*LQS**2.+6.5057411*LTK*LQS) ENDIF RHOS = MIN(MAX(RHOS,RHOIMIN),RHOIMAX) VS1D = 0. ENDIF IF (NS1D.LT.NSMALL) THEN LTK = LOG(TK1D) LQS = -1.*LOG(QS1D) MDS = EXP(-123.23898+40.74706*LTK-3.0333477*LTK**2.- & 0.31219981*LQS+0.0012798222*LQS**2.)/1.E3 NS1D = 1.E9*QS1D*V2M3/RHOS/MDS**3. ENDIF AMS = PI*RHOS/6. S3M1D = QS1D/AMS IF (AGG_SHAPE.EQ.0) THEN SASPR = 1. FS1D = 0. ELSEIF (AGG_SHAPE.EQ.1) THEN SASPR = FS1D/S3M1D IF (SASPR.LT.SASMIN) THEN SASPR = SASMIN FS1D = SASPR*S3M1D ELSEIF (SASPR.GT.SASMAX) THEN SASPR = SASMAX FS1D = SASPR*S3M1D ENDIF ENDIF IF (S2M1D.GE.ASMALL.AND.S3M1D.GE.ISMALL) THEN KDX = MAX(0.,S2M1D**3./(S3M1D**2.*NS1D)) IF (KDX.GT.KCSMAX) THEN KDX = KCSMIN S2M1D = (KDX*NS1D*S3M1D**2.)**THRD ELSEIF (KDX.LT.KCSMIN) THEN KDX = KCSMIN S2M1D = (KDX*NS1D*S3M1D**2.)**THRD ENDIF AFAS = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAS = MIN(MAX(AFAS,AFAMIN),AFAMAX) ! LAMS = (NS1D*(AFAS+1.)*(AFAS+2.)*(AFAS+3.)/S3M1D)**THRD LAMS = (AFAS+3.)*S2M1D/S3M1D ELSE IF (AFAS_3M.EQ.0) THEN AFAS = AFAS0 S2M1D = 0. ELSEIF (AFAS_3M.EQ.1) THEN AFAS = AFAS0 S2M1D = (KCSMIN*NS1D*S3M1D**2.)**THRD ELSEIF (AFAS_3M.EQ.2) THEN BDS = MIN(MAX((S3M1D/NS1D)**THRD*1.E3,DSMIN*1.E3),DSMAX*1.E3) IF (TK1D.GE.TK0C) THEN FDS = -0.21911541+1.2739845*BDS+0.10141003*LOG(NS1D)+ & 0.30063818*BDS**2.-4.3857765E-3*LOG(NS1D)**2.- & 7.8801732E-2*BDS*LOG(NS1D) ELSE IF (QC1D.GE.1.E-8) THEN FDS = -1.1527014+2.9067645*BDS+0.25316062*LOG(NS1D)- & 0.17768557*BDS**2.-0.013117292*LOG(NS1D)**2.- & 0.17020429*BDS*LOG(NS1D) ELSE FDS = -0.2813929+1.7275463*BDS+0.045550156*LOG(NS1D)-& 0.16526226*BDS**2.-1.7699916E-3*LOG(NS1D)**2.- & 4.6441257E-2*BDS*LOG(NS1D) ENDIF ENDIF KDX = MAX(KCSMIN,MIN(KCSMAX,(BDS/FDS)**3.)) AFAS = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAS = MIN(MAX(AFAS,AFAMIN),AFAMAX) ! AFAS = MAX(AFAMIN,4.5*TANH(0.5*(BDS-5.))+5.5) S2M1D = 0. ENDIF LAMS = (EXP(GAMLN(AFAS+4.)-GAMLN(AFAS+1.))*NS1D/S3M1D)**THRD ENDIF LAMSMIN = (EXP(GAMLN(AFAS+4.)-GAMLN(AFAS+1.)))**THRD/DSMAX LAMSMAX = (EXP(GAMLN(AFAS+4.)-GAMLN(AFAS+1.)))**THRD/DSMIN IF (LAMS.LT.LAMSMIN) THEN LAMS = LAMSMIN NS1D = S3M1D*EXP(GAMLN(AFAS+1.)-GAMLN(AFAS+4.)+3.*LOG(LAMS)) ELSEIF (LAMS.GT.LAMSMAX) THEN LAMS = LAMSMAX NS1D = S3M1D*EXP(GAMLN(AFAS+1.)-GAMLN(AFAS+4.)+3.*LOG(LAMS)) ENDIF MVDS = (EXP(GAMLN(AFAS+4.)-GAMLN(AFAS+1.)))**THRD/LAMS IF (ICE_VTS.EQ.0) THEN AVS = AVS0 BVS = BVS0 ELSEIF (ICE_VTS.EQ.1.OR.ICE_VTS.EQ.2) THEN KINV = (1.72E-5*(393./(TK1D+120.))*(TK1D/TK0C)**1.5)/RHO IF (AGG_SHAPE.EQ.0) THEN AMS3 = AMS BMS3 = BMS AAS3 = PI/4. BAS3 = 2. LAMD = LAMS ELSEIF (AGG_SHAPE.EQ.1) THEN TC1D = TK1D-TK0C IF (ICE_VTS.EQ.1) THEN TBIN = MIN(MAX(INT(ABS(TC1D)/5.)-4,0),8) IF (MVDS.LT.1.5E-4) DBIN = 0 IF (MVDS.GE.1.5E-4.AND.MVDS.LT.3.E-4) DBIN = 1 IF (MVDS.GE.3.E-4.AND.MVDS.LT.1.E-3) DBIN = 2 IF (MVDS.GE.1.E-3.AND.MVDS.LT.3.E-3) DBIN = 3 IF (MVDS.GE.3.E-3) DBIN = 4 DBIN = MIN(MAX(DBIN,0),4) AMS3 = AMS1(TBIN*5+DBIN) BMS3 = BMS1(TBIN*5+DBIN) AAS3 = AAS1(TBIN*5+DBIN) BAS3 = BAS1(TBIN*5+DBIN) LAMD = MIN((AMS3*LAMS**3./AMS*EXP(GAMLN(BMS3+AFAS+1.)- & GAMLN(AFAS+4.)))**(1./BMS3),LAMS) ELSEIF (ICE_VTS.EQ.2) THEN BMS3 = 2.4+0.0085*MAX(TC1D,-65.) BAS3 = 2-0.19+0.0056*MAX(TC1D,-65.) AMS2 = 0.0102+0.00013*MAX(TC1D,-65.) AAS2 = 0.29+0.0035*MAX(TC1D,-65.) AMS3 = ((AMS2*(100.*MVDS)**BMS3)/1.E3)/(MVDS**BMS3) AAS3 = ((AAS2*(100.*MVDS)**BAS3)/1.E4)/(MVDS**BAS3) LAMD = LAMS ENDIF ENDIF SBA1 = BMS3+2.-BAS3 BEST = 2.*G*NS1D*AMS3*EXP(GAMLN(SBA1+AFAS+1.)-GAMLN(AFAS+1.)- & SBA1*LOG(LAMD))/(KINV**2.*AAS3) C1X2 = VTC1*BEST**5.E-1 VTB1 = C1X2/(1.+C1X2)**5.E-1/((1.+C1X2)**5.E-1-1.)/2. VTA1 = VTC2*((1+C1X2)**5.E-1-1.)**2./BEST**VTB1 AVS = VTA1*KINV**(1.-2.*VTB1)*(2.*AMS3*G/AAS3)**VTB1 BVS = VTB1*(BMS3+2.-BAS3)-1. ENDIF END SUBROUTINE SOLVE_AFAS !====================================================================== !====================================================================== SUBROUTINE SOLVE_AFAG(TK1D,RHO,QG1D,QC1D,NG1D,VG1D,G2M1D,LAMG, & AFAG,MVDG,RHOG,AMG,AVG,BVG) !====================================================================== IMPLICIT NONE REAL :: TK1D,RHO,QG1D,QC1D,NG1D,VG1D,G2M1D,G3M1D,RHOG,LAMG,AFAG, & AMG,AVG,BVG,BEST0,LAMGMIN,LAMGMAX,KDX,KINV,GG1,C1X2,VTA1,& VTB1,BEST,MVDG,GMLR,BDG,FDG,LTK,LQG,MDG IF (ICE_RHOG.EQ.0) THEN RHOG = RHOG1 VG1D = 0. ELSEIF (ICE_RHOG.EQ.1.OR.ICE_RHOG.EQ.2) THEN RHOG = QG1D/(VG1D+ISMALL) IF (TK1D.GE.TK0C) THEN IF (RHOG.GT.RHOH) THEN RHOG = RHOH ENDIF ELSE IF (RHOG.GT.RHOG0) THEN RHOG = RHOG0 ENDIF ENDIF IF (RHOG.LT.RHOIMIN) THEN RHOG = RHOIMIN ENDIF VG1D = QG1D/RHOG ENDIF IF (NG1D.LT.NSMALL) THEN LTK = LOG(TK1D) LQG = -1.*LOG(QG1D) MDG = EXP(-2205.8027+1225.8046*LTK-226.27995*LTK**2.+ & 13.929644*LTK**3.-0.25*LQG)/1.E3 NG1D = 1.E9*QG1D*V2M3/RHOG/MDG**3. ENDIF AMG = PI*RHOG/6. G3M1D = QG1D*V2M3/RHOG IF (G2M1D.GE.ASMALL.AND.G3M1D.GE.ISMALL) THEN KDX = MAX(0.,G2M1D**3./(G3M1D**2.*NG1D)) IF (KDX.GT.KCGMAX) THEN KDX = KCGMIN G2M1D = (KDX*NG1D*G3M1D**2.)**THRD ELSEIF (KDX.LT.KCGMIN) THEN KDX = KCGMIN G2M1D = (KDX*NG1D*G3M1D**2.)**THRD ENDIF AFAG = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAG = MIN(MAX(AFAG,AFAMIN),AFAMAX) ! LAMG = (NG1D*(AFAG+1.)*(AFAG+2.)*(AFAG+3.)/G3M1D)**THRD LAMG = (AFAG+3.)*G2M1D/G3M1D ELSE IF (AFAG_3M.EQ.0) THEN AFAG = AFAG0 G2M1D = 0. ELSEIF (AFAG_3M.EQ.1) THEN AFAG = AFAG0 G2M1D = (KCGMIN*NG1D*G3M1D**2.)**THRD ELSEIF (AFAG_3M.EQ.2) THEN BDG = MIN(MAX((G3M1D/NG1D)**THRD*1.E3,DGMIN*1.E3),DGMAX*1.E3) IF (TK1D.GE.TK0C) THEN FDG = 0.58006354+0.79661229*BDG-0.18394382*LOG(NG1D)+ & 0.067371044*BDG**2.+9.832945E-3*LOG(NG1D)**2.+ & 0.12433055*BDG*LOG(NG1D) ELSE IF (QC1D.GE.1.E-8) THEN FDG = 0.17363469+1.5044291*BDG-0.050639722*LOG(NG1D)+& 0.015101052*BDG**2.+2.5974719E-3*LOG(NG1D)**2.+& 0.01961464*BDG*LOG(NG1D) ELSE FDG = 0.59259317-0.89933515*BDG+2.0893032*BDG**2.- & 0.50305755*BDG**3.-2.2446793E-2*LOG(NG1D)+ & 2.7589047E-3*LOG(NG1D)**2. ENDIF ENDIF KDX = MAX(KCGMIN,MIN(KCGMAX,(BDG/FDG)**3.)) AFAG = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAG = MIN(MAX(AFAG,AFAMIN),AFAMAX) ! AFAG = MAX(AFAMIN,5.5*TANH(0.7*(BDG-4.5))+8.5) G2M1D = 0. ENDIF LAMG = (EXP(GAMLN(AFAG+4.)-GAMLN(AFAG+1.))*NG1D/G3M1D)**THRD ENDIF LAMGMIN = (EXP(GAMLN(AFAG+4.)-GAMLN(AFAG+1.)))**THRD/DGMAX LAMGMAX = (EXP(GAMLN(AFAG+4.)-GAMLN(AFAG+1.)))**THRD/DGMIN IF (LAMG.LT.LAMGMIN) THEN LAMG = LAMGMIN NG1D = G3M1D*EXP(GAMLN(AFAG+1.)-GAMLN(AFAG+4.)+3.*LOG(LAMG)) ELSEIF (LAMG.GT.LAMGMAX) THEN LAMG = LAMGMAX NG1D = G3M1D*EXP(GAMLN(AFAG+1.)-GAMLN(AFAG+4.)+3.*LOG(LAMG)) ENDIF MVDG = (EXP(GAMLN(AFAG+4.)-GAMLN(AFAG+1.)))**THRD/LAMG IF (ICE_VTG.EQ.0) THEN AVG = AVG0 BVG = BVG0 ELSEIF (ICE_VTG.EQ.1) THEN KINV = (1.72E-5*(393./(TK1D+120.))*(TK1D/TK0C)**1.5)/RHO BEST0 = 2.*G*NG1D/(KINV**2.) GG1 = GAMLN(AFAG+1.) BEST = BEST0*AMG*EXP(GAMLN(BMG+AFAG+1.)-GG1-BMG*LOG(LAMG))/AAW C1X2 = VTC1*BEST**5.E-1 VTB1 = C1X2/(1.+C1X2)**5.E-1/((1.+C1X2)**5.E-1-1.)/2. VTA1 = VTC2*((1+C1X2)**5.E-1-1.)**2./BEST**VTB1 AVG = VTA1*KINV**(1.-2.*VTB1)*(2.*AMG*G/AAW)**VTB1 BVG = VTB1*BMG-1. ENDIF END SUBROUTINE SOLVE_AFAG !====================================================================== ! !====================================================================== SUBROUTINE SOLVE_AFAH(TK1D,RHO,QH1D,NH1D,H2M1D,LAMH,AFAH,MVDH, & AVH,BVH) !====================================================================== IMPLICIT NONE REAL :: TK1D,RHO,QH1D,NH1D,H2M1D,H3M1D,LAMH,AFAH,AVH,BVH,BEST0, & GH1,KDX,MVDH,LAMHMIN,LAMHMAX,KINV,C1X2,VTA1,VTB1,BEST, & BDH,FDH,LTK,LQH IF (NH1D.LT.NSMALL) THEN LTK = LOG(TK1D) LQH = -1.*LOG(QH1D) NH1D = EXP(22.929406-4.2328364*LTK+0.30647567*LTK**2.- & 0.009233271*LTK**3.-0.25*LQH) ENDIF H3M1D = QH1D*V2M3/RHOH IF (H2M1D.GE.ASMALL.AND.H3M1D.GE.ISMALL) THEN KDX = MAX(0.,H2M1D**3./(H3M1D**2.*NH1D)) IF (KDX.GT.KCHMAX) THEN KDX = KCHMIN H2M1D = (KDX*NH1D*H3M1D**2.)**THRD ELSEIF (KDX.LT.KCHMIN) THEN KDX = KCHMIN H2M1D = (KDX*NH1D*H3M1D**2.)**THRD ENDIF AFAH = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAH = MIN(MAX(AFAH,AFAMIN),AFAMAX) ! LAMH = (NH1D*(AFAH+1.)*(AFAH+2.)*(AFAH+3.)/H3M1D)**THRD LAMH = (AFAH+3.)*H2M1D/H3M1D ELSE IF (AFAH_3M.EQ.0) THEN AFAH = AFAH0 H2M1D = 0. ELSEIF (AFAH_3M.EQ.1) THEN AFAH = AFAH0 H2M1D = (KCHMIN*NH1D*H3M1D**2.)**THRD ELSEIF (AFAH_3M.EQ.2) THEN BDH = MIN(MAX((H3M1D/NH1D)**THRD*1.E3,DHMIN*1.E3),DHMAX*1.E3) IF (TK1D.GE.TK0C) THEN FDH = 1.157754+0.37852874*BDH-0.11129737*LOG(NH1D)+ & 0.13929599*BDH**2.+8.1105237E-3*LOG(NH1D)**2.+ & 5.7432113E-2*BDH*LOG(NH1D) ELSE FDH = -0.48246793+2.0407077*BDH+2.2262969E-2*LOG(NH1D)- & 0.158389*BDH**2.-5.5545804E-3*LOG(NH1D)**2.+ & 2.9443577E-2*BDH*LOG(NH1D) ENDIF KDX = MAX(KCHMIN,MIN(KCHMAX,(BDH/FDH)**3.)) AFAH = (6.*KDX-3.+SQRT(8.*KDX+1.))/(2.-2.*KDX) AFAH = MIN(MAX(AFAH,AFAMIN),AFAMAX) ! AFAH = MAX(AFAMIN,3.7*TANH(0.3*(BDH-9.))+6.5) H2M1D = 0. ENDIF LAMH = (EXP(GAMLN(AFAH+4.)-GAMLN(AFAH+1.))*NH1D/H3M1D)**THRD ENDIF LAMHMIN = (EXP(GAMLN(AFAH+4.)-GAMLN(AFAH+1.)))**THRD/DHMAX LAMHMAX = (EXP(GAMLN(AFAH+4.)-GAMLN(AFAH+1.)))**THRD/DHMIN IF (LAMH.LT.LAMHMIN) THEN LAMH = LAMHMIN NH1D = H3M1D*EXP(GAMLN(AFAH+1.)-GAMLN(AFAH+4.)+3.*LOG(LAMH)) ELSEIF (LAMH.GT.LAMHMAX) THEN LAMH = LAMHMAX NH1D = H3M1D*EXP(GAMLN(AFAH+1.)-GAMLN(AFAH+4.)+3.*LOG(LAMH)) ENDIF MVDH = (EXP(GAMLN(AFAH+4.)-GAMLN(AFAH+1.)))**THRD/LAMH IF (ICE_VTH.EQ.0) THEN AVH = AVH0 BVH = BVH0 ELSEIF (ICE_VTH.EQ.1) THEN KINV = (1.72E-5*(393./(TK1D+120.))*(TK1D/TK0C)**1.5)/RHO BEST0 = 2.*G*NH1D/(KINV**2.) GH1 = GAMLN(AFAH+1.) BEST = BEST0*AMH*EXP(GAMLN(BMH+AFAH+1.)-GH1-BMH*LOG(LAMH))/AAW C1X2 = VTC1*BEST**5.E-1 VTB1 = C1X2/(1.+C1X2)**5.E-1/((1.+C1X2)**5.E-1-1.)/2. VTA1 = VTC2*((1+C1X2)**5.E-1-1.)**2./BEST**VTB1 AVH = VTA1*KINV**(1.-2.*VTB1)*(2.*AMH*G/AAW)**VTB1 BVH = VTB1*BMH-1. ENDIF END SUBROUTINE SOLVE_AFAH !====================================================================== ! !====================================================================== SUBROUTINE AERO_CONST(CCNTY) !====================================================================== IMPLICIT NONE INTEGER, INTENT(IN) :: CCNTY ! 1.MARINE(WHITBY,1978);2.CLEAN CONTINENTAL(WHITBY,1978); ! 3.AVERAGE CONTINENTAL(WHITBY,1978); ! 4.URBAN(WHITBY,1978); 5.ICE NUECEI INTEGER, PARAMETER :: ID_SENS = 1 ! sensetivity test ID INTEGER, PARAMETER :: N1 = NCCN, N2 = 5 ! AEROSOL MODE AND COMPONENT INTEGER :: IA,I,K,J,IM,IV,IB INTEGER, DIMENSION(NAER) :: IDAER REAL, DIMENSION(NAER) :: XMASS REAL, DIMENSION(N2) :: ASHX,CMASX,DNASX,AVANX ! AEROSOL SCALE HEIGHT. 3.57km-continental, 0.8km-marine; ! MOLAR WEIGHT OF AEROSOL (kg/mol); ! DENSITY OF SOLID AEROSOL (kg/m^3); ! VAN'T HOFF FACTOR REAL, DIMENSION(N1,N2) :: ZCCNX,CNMODX,CNSTDX ! AEROSOL TRI-MODAL SPECTRUM PARAMETERS REAL, DIMENSION(N2,3) :: SENX ! SENSITIVE TEST RATIO REAL :: XLOG3,XMTOT,XMAS,DUST_IN0,DUST_IN,DAIR,DZ0,ZH ! Height of grid point [m] DOUBLE PRECISION :: DMODE,D2STDV,DLNXX,DERF,DSQRT,DLOG DATA ZCCNX/3.4E8,6.E7,3.E6, 1.E9,8.E8,7.2E5, 6.4E9,2.3E9,3.2E6, &! Cheng et al. 2007, Table 1, number concentration 1.06E11,3.2E10,5.4E6, 4.02E-11,4.5E9,5.4E6/ DATA CNMODX/5.E-9,3.55E-8,3.1E-7,8.E-9,3.3E-8,4.6E-7,7.5E-9, &! Cheng et al. 2007, Table 1, mean size 3.8E-8,5.1E-7,7.E-9,2.7E-8,4.3E-7,4.5E-9,2.35E-8, & 5.5E-7/ DATA CNSTDX/0.47,0.69,0.99,0.47,0.74,0.79,0.53,0.69,0.77,0.59, &! Cheng et al. 2007, Table 1, standard deviation 0.77,0.79,0.53,0.75,0.74/ DATA ASHX/8.E2,3.57E3,3.57E3,2.E3,1.E4/ ! SCALING HEIGHT DATA CMASX/0.13214,0.13214,0.13214,0.13214,0.13214/ ! MOLAR WEIGHT OF AEROSOL (kg/mol) DATA DNASX/1769.,1769.,1769.,1769.,1769./ ! DENSITY OF SOLID AEROSOL (kg/m^3) DATA AVANX/3.,3.,3.,3.,3./ ! VAN'T HOFF FACTOR DATA SENX/1.,1.,1.,1.,1.,10.,10.,10.,10.,10.,0.1,0.1,0.1,0.1,0.1/ ! 1st, 2nd, and 3rd set for sensetivity test IF (CCNTY.EQ.1) THEN DNC0 = 46635.361; DNC1 = -25567.933 DNC2 = 4674.3534; DNC3 = -284.84377 ELSEIF (CCNTY.EQ.2) THEN DNC0 = -48575.888; DNC1 = 26634.802 DNC2 = -4862.961; DNC3 = 295.7886 ELSEIF (CCNTY.EQ.3) THEN DNC0 = -63646.587; DNC1 = 34778.66 DNC2 = -6326.6327; DNC3 = 383.28849 ELSEIF (CCNTY.EQ.4) THEN DNC0 = -6069.3697; DNC1 = 3667.889 DNC2 = -725.3412; DNC3 = 47.260735 ENDIF IBAER(1) = 1 DO I = 2,NAER IBAER(I) = IBAER(I-1)+NAERN(I-1) ENDDO TBLRC(1) = DLOG(1.D-9) ! 0.001 um TBLRC(NTBXA) = DLOG(1.D-3) ! 1000. um DO I = 2,NTBXA-1 TBLRC(I) = TBLRC(1)+DBLE(I-1)/DBLE(NTBXA-1)*(TBLRC(NTBXA)- & TBLRC(1)) ! ln(rc) ENDDO IDAER = (/CCNTY,5/) PRINT *, 'AEROSOL ID =',IDAER DO J = 1,NAER ASH(J) = ASHX(IDAER(J)) DNAS(J) = DNASX(IDAER(J)) CMAS(J) = CMASX(IDAER(J)) AVAN(J) = AVANX(IDAER(J)) BETA1(J) = CMW*DNAS(J)*AVAN(J)/(CMAS(J)*RHOW) XMTOT = 0. DO IM = 1,N1 ZCCN(IM,J) = ZCCNX(IM,IDAER(J))*1. ! change the factor to change the initial concentration CNMOD(IM,J) = CNMODX(IM,IDAER(J)) CNSTD(IM,J) = CNSTDX(IM,IDAER(J)) XLOG3 = CNMOD(IM,J)*CNMOD(IM,J)*CNMOD(IM,J) WMAS(IM,J) = ZCCN(IM,J)*XLOG3*EXP(4.5*CNSTD(IM,J)**2.) ! factor for mass, the 3rd momentum RFACT(IM,J) = 1./(C4PI3*DNAS(J)*XLOG3*EXP(4.5* & CNSTD(IM,J)**2.)) XMTOT = XMTOT+WMAS(IM,J) ENDDO DO IM = 1,N1 WMAS(IM,J) = WMAS(IM,J)/XMTOT ENDDO DO I = 1,NTBXA TBLXF(I,J) = DBLE(0.) ENDDO DO IM = 1,N1 DMODE = DLOG(DBLE(CNMOD(IM,J)))+DBLE(3.)*DBLE(CNSTD(IM,J))*DBLE(CNSTD(IM,J)) D2STDV = DSQRT(DBLE(2.))*DBLE(CNSTD(IM,J)) DO I = 1,NTBXA DLNXX = (TBLRC(I)-DMODE)/D2STDV TBLXF(I,J) = TBLXF(I,J)+DBLE(WMAS(IM,J))*DBLE(0.5)* & (1.D+0+DERF(DLNXX)) ENDDO ENDDO ! IM=1,N1 DO I = 1,NTBXA TBLXF(I,J) = DMAX1(1.D-20,TBLXF(I,J)) ENDDO CALL FIND_RC0(1.D-6,CNMOD(1,J),CNSTD(1,J),WMAS(1,J),RXMIN(J), & TBLXF(1,J),TBLRC) ! lower limit of dry aerosol fraction is set as 1.E-6 WRITE(*,*)'FOR DRY AEROSOL FRACTION BE > 1.E-6' WRITE(*,'(A30,I1,A8,F6.4,A3)')'MINIMUN R-CUTOFF FOR AEROSOL' & ,J,' IS SET ',RXMIN(J)*1.E6,' um' WRITE(*,*)'SENSITIVITY TEST ID:',ID_SENS !----------------- FACTOR FOR SENSETIVITY ------------------------------ SENS(J) = SENX(IDAER(J),ID_SENS) ENDDO END SUBROUTINE AERO_CONST !====================================================================== !====================================================================== SUBROUTINE INIT_AEROSOL(P,RHO,DZ,XLAND,QAERO,IDS,IDE,JDS,JDE,KDS,& KDE,IMS,IME,JMS,JME,KMS,KME,ITS,ITE,JTS,JTE,KTS,KTE) !====================================================================== IMPLICIT NONE INTEGER, INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME, & KMS,KME,ITS,ITE,JTS,JTE,KTS,KTE INTEGER :: ITF,JTF,IA,I,K,J,IM,K_PBL,IV,IB INTEGER, DIMENSION(NAER) :: IDAER REAL, INTENT(IN), DIMENSION(IMS:IME,JMS:JME) :: XLAND ! 1:land ; 2:ocean REAL, INTENT(IN), DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE) :: P,RHO,DZ REAL, DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE,NAERT) :: QAERO REAL, DIMENSION(ITS:ITE,JTS:JTE,NAER) :: RASH REAL, DIMENSION(NAER) :: XMASS REAL, PARAMETER :: P_PBL = 85000. ! set PBL pressure [Pa] REAL, PARAMETER :: BACKIN = 4.E5 ! initial ice nuclei number concentration (#/kg) REAL :: XLOG3,XMTOT,XMAS,DUST_IN0,DUST_IN,DAIR,DZ0,ZH ! Height of grid point [m] DOUBLE PRECISION :: DMODE,D2STDV,DLNXX,DERF,DSQRT,DLOG ITF = MIN(ITE,IDE-1) JTF = MIN(JTE,JDE-1) DO IA = 1,NAERT DO K = KTS,KTE DO J = JTS,JTF DO I = ITS,ITF QAERO(I,K,J,IA) = 0. ENDDO ENDDO ENDDO ENDDO DO J = 1,NAER XMASS(J) = 0. DO IM = 1,NCCN XMAS = ZCCN(IM,J)*CNMOD(IM,J)**3.*EXP(4.5*CNSTD(IM,J)**2.) XMASS(J) = XMASS(J)+XMAS ENDDO XMASS(J) = XMASS(J)*C4PI3*DNAS(J)*SENS(J) ENDDO IF (ID_DUST*ID_IN.NE.0) THEN ! IF DUST EXISTS DUST_IN0 = 0. ! DO IM = 1,NCCN ! DUST_IN0 = DUST_IN0+ZCCN(IM,ID_DUST) ! ENDDO ELSE ! NO DUST DUST_IN0 = BACKIN ENDIF DO J = JTS,JTF DO I = ITS,ITF K_PBL = 1 DO K = KTS,KTE IV = 1 DO IA = 1,NAER QAERO(I,K,J,IV) = XMASS(IA)/RHO(I,K,J) ! dry aerosol [kg/kg] RASH (I,J,IA) = 1./ASH(IA) ! increase scale height if necessary DUST_IN = DUST_IN0 ! dust # over land IF (IA.EQ.ID_IN) THEN ! IN IF (ID_DUST.NE.0) THEN ! DUST IS CONSIDERED AS IN ! IF (XLAND(I,J).EQ.2.) THEN ! over ocean ! DUST_IN = DUST_IN*0.01 ! dust # ! RASH (I,J,IA) = RASH (I,J,IA)*0.5 ! double scale height ! ENDIF ELSE DUST_IN = BACKIN ! set to background value RASH (I,J,IA) = RASH (I,J,IA)*0.5 ! double scale height ENDIF DAIR = RHO(I,K,J) DUST_IN = MAX(BACKIN,DUST_IN) QAERO(I,K,J,IV) = DUST_IN/DAIR ! 0.5 #/L [#/kg]*[kPa] ! tsai ELSE ! dust, NH42SO4, chemicals.... ! IF (XLAND(I,J).EQ.2.) THEN ! over ocean ! IF (IA.EQ.ID_NH42SO4) THEN ! ! QAERO(I,K,J,IV) = QAERO(I,K,J,IV)*0.05 ! reduce it ! RASH(I,J,IA) = RASH(I,J,IA)*0.5 ! double scale height ! ELSE ! dust or IN, not sulfate (the same) ! QAERO(I,K,J,IV) = QAERO(I,K,J,IV)*0.01 ! reduce it ! RASH(I,J,IA) = RASH (I,J,IA)*0.5 ! double scale height ! DUST_IN = DUST_IN0*0.01 ! dust # ! ENDIF ! ENDIF ENDIF IV = IV+NAERN(IA) ENDDO IF (P(I,K,J).LT.P_PBL) K_PBL = K ENDDO K_PBL = MIN(K_PBL,KTE-3+1) ! at least 3 layers are within PBL ZH = 0. DZ0 = 0. DO K = K_PBL-1,1,-1 DAIR = RHO(I,K,J) ZH = ZH+(DZ0+DZ(I,K,J))*0.5 DZ0 = DZ(I,K,J) IV = 1 DO IA = 1,NAER QAERO(I,K,J,IV)=QAERO(I,K,J,IV)*EXP(-RASH(I,J,IA)*ZH) ! dry aerosol [kg/kg]*[kPa] IV = IV+NAERN(IA) ENDDO ENDDO IV = 0 DO IA = 1,NAER DO IB = 2,MIN(2,NAERN(IA)) ! if NV>2 DO K = KTS,KTE QAERO(I,K,J,IV+IB) = QAERO(I,K,J,IV+1) ! set NV2: total = dry aerosol ENDDO ENDDO IV = IV+NAERN(IA) ENDDO ENDDO ! FOR I LOOPS ENDDO ! FOR J LOOPS PRINT *,'aerosols_init,I,J,NAER',(ITS+ITE)/2,(JTS+JTE)/2,NAER DO K = KTS,KTE-1 WRITE(*,'(I2,X,20(E12.6,X))') K,(QAERO((ITS+ITE)/2,K, & (JTS+JTE)/2,IV),IV=1,NAERT) ENDDO RETURN END SUBROUTINE INIT_AEROSOL !====================================================================== ! !======================================================================= SUBROUTINE MP_NTU(ITIMESTEP,TH,P,DZ,W,PII,DT_MP,SR,QV,QC,QR,QI, & QS,QG,QH,NC,NR,NI,NS,NG,NH,QDCN,QTCN,QCCN,QRCN,QNIN, & FI,FS,VI,VS,VG,AI,AS,AG,AH,I3M,RAINNC,RAINNCV,SNOWNC, & SNOWNCV,GRAPNC,GRAPNCV,HAILNC,HAILNCV,IDS,IDE,JDS,JDE,& KDS,KDE,IMS,IME,JMS,JME,KMS,KME,ITS,ITE,JTS,JTE,KTS, & KTE) !====================================================================== IMPLICIT NONE INTEGER, INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,IMS,IME,JMS,JME, & KMS,KME,ITS,ITE,JTS,JTE,KTS,KTE,ITIMESTEP REAL, INTENT(IN) :: DT_MP REAL, DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: RAINNC,RAINNCV,& SNOWNC,SNOWNCV,GRAPNC,GRAPNCV,HAILNC,HAILNCV,SR REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(IN) :: PII,DZ,W,P REAL, DIMENSION(IMS:IME,KMS:KME,JMS:JME), INTENT(INOUT) :: & TH,QV,QC,QR,QI,QS,QG,QH,NC,NR,NI,NS,NG,NH,VI,VS,VG,FI,FS, & AI,AS,AG,AH,I3M,QDCN,QTCN,QCCN,QRCN,QNIN INTEGER :: I,K,J,NK,ITF,JTF REAL :: DT,DTMN REAL, DIMENSION(KTS:KTE) :: TK1D,P1D,W1D,S1D,DZ1D,QV1D,QC1D,QR1D,& QI1D,QS1D,QG1D,QH1D,NC1D,NR1D,NI1D,NS1D,NG1D,NH1D,VI1D, & VS1D,VG1D,FI1D,FS1D,AI1D,AS1D,AG1D,AH1D,I3M1D,QDCN1D, & QTCN1D,QCCN1D,QRCN1D,QNIN1D REAL, DIMENSION(ITS:ITE,JTS:JTE) :: CLODNCV,ICENCV,VINCV,VSNCV, & VGNCV,FINCV,FSNCV,AINCV,ASNCV,AGNCV,AHNCV,I3MNCV REAL, DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE) :: TK DT = DT_MP DTMN = DT_MP/60. ITF = MIN(ITE,IDE-1) JTF = MIN(JTE,JDE-1) DO I = ITS,ITF DO J = JTS,JTF CLODNCV(I,J) = 0.; RAINNCV(I,J) = 0.; ICENCV(I,J) = 0. SNOWNCV(I,J) = 0.; GRAPNCV(I,J) = 0.; HAILNCV(I,J) = 0. VINCV(I,J) = 0.; VSNCV(I,J) = 0.; VGNCV(I,J) = 0. FINCV(I,J) = 0.; FSNCV(I,J) = 0.; AINCV(I,J) = 0. ASNCV(I,J) = 0.; AGNCV(I,J) = 0.; AHNCV(I,J) = 0. I3MNCV(I,J) = 0.; SR(I,J) = 0. !----------------- write data from 3D to 1D and assign upside down ----- DO K = KTS,KTE NK = KTE-K+1 S1D(K) = 0. TK1D(K) = TH(I,NK,J)*PII(I,NK,J) W1D(K) = 0.5*(W(I,NK,J)+W(I,NK+1,J)) ! W at half level QV1D(K) = QV(I,NK,J) DZ1D(K) = DZ(I,NK,J) P1D(K) = P(I,NK,J) QC1D(K) = QC(I,NK,J) QR1D(K) = QR(I,NK,J) QI1D(K) = QI(I,NK,J) QS1D(K) = QS(I,NK,J) QG1D(K) = QG(I,NK,J) QH1D(K) = QH(I,NK,J) NC1D(K) = NC(I,NK,J) NR1D(K) = NR(I,NK,J) NI1D(K) = NI(I,NK,J) NS1D(K) = NS(I,NK,J) NG1D(K) = NG(I,NK,J) NH1D(K) = NH(I,NK,J) VI1D(K) = VI(I,NK,J) VS1D(K) = VS(I,NK,J) VG1D(K) = VG(I,NK,J) FI1D(K) = FI(I,NK,J) FS1D(K) = FS(I,NK,J) AI1D(K) = AI(I,NK,J) AS1D(K) = AS(I,NK,J) AG1D(K) = AG(I,NK,J) AH1D(K) = AH(I,NK,J) I3M1D(K) = I3M(I,NK,J) QDCN1D(K) = MAX(RLIMIT,QDCN(I,NK,J)) QTCN1D(K) = MAX(RLIMIT,QTCN(I,NK,J)) QCCN1D(K) = MAX(RLIMIT,QCCN(I,NK,J)) QRCN1D(K) = MAX(RLIMIT,QRCN(I,NK,J)) QNIN1D(K) = MAX(RLIMIT,QNIN(I,NK,J)) ENDDO CALL NTU_MICRO(TK1D,QV1D,DZ1D,P1D,S1D,W1D,SR(I,J),ICENCV(I,J),& VINCV(I,J),VSNCV(I,J),VGNCV(I,J),FINCV(I,J),FSNCV(I,J), & AINCV(I,J),ASNCV(I,J),AGNCV(I,J),AHNCV(I,J),I3MNCV(I,J), & CLODNCV(I,J),RAINNC(I,J),RAINNCV(I,J),SNOWNC(I,J), & SNOWNCV(I,J),GRAPNC(I,J),GRAPNCV(I,J),HAILNC(I,J), & HAILNCV(I,J),DT,DTMN,QC1D,QR1D,QI1D,QS1D,QG1D,QH1D,NC1D, & NR1D,NI1D,NS1D,NG1D,NH1D,VI1D,VS1D,VG1D,FI1D,FS1D,AI1D, & AS1D,AG1D,AH1D,I3M1D,QDCN1D,QTCN1D,QCCN1D,QRCN1D,QNIN1D, & KTS,KTE) DO K = KTS,KTE NK = KTE-K+1 TH(I,NK,J) = TK1D(K)/PII(I,NK,J) QV(I,NK,J) = MAX(0.,QV1D(K)) QC(I,NK,J) = MAX(0.,QC1D(K)) QR(I,NK,J) = MAX(0.,QR1D(K)) QI(I,NK,J) = MAX(0.,QI1D(K)) QS(I,NK,J) = MAX(0.,QS1D(K)) QG(I,NK,J) = MAX(0.,QG1D(K)) QH(I,NK,J) = MAX(0.,QH1D(K)) NC(I,NK,J) = MAX(0.,NC1D(K)) NR(I,NK,J) = MAX(0.,NR1D(K)) NI(I,NK,J) = MAX(0.,NI1D(K)) NS(I,NK,J) = MAX(0.,NS1D(K)) NG(I,NK,J) = MAX(0.,NG1D(K)) NH(I,NK,J) = MAX(0.,NH1D(K)) VI(I,NK,J) = MAX(0.,VI1D(K)) VS(I,NK,J) = MAX(0.,VS1D(K)) VG(I,NK,J) = MAX(0.,VG1D(K)) FI(I,NK,J) = MAX(0.,FI1D(K)) FS(I,NK,J) = MAX(0.,FS1D(K)) AI(I,NK,J) = MAX(0.,AI1D(K)) AS(I,NK,J) = MAX(0.,AS1D(K)) AG(I,NK,J) = MAX(0.,AG1D(K)) AH(I,NK,J) = MAX(0.,AH1D(K)) I3M(I,NK,J) = MAX(0.,I3M1D(K)) QDCN(I,NK,J) = MAX(0.,QDCN1D(K)) QTCN(I,NK,J) = MAX(0.,QTCN1D(K)) QCCN(I,NK,J) = MAX(0.,QCCN1D(K)) QRCN(I,NK,J) = MAX(0.,QRCN1D(K)) QNIN(I,NK,J) = MAX(0.,QNIN1D(K)) ENDDO ENDDO ENDDO END SUBROUTINE MP_NTU !====================================================================== ! !====================================================================== SUBROUTINE NTU_MICRO(TK3D,QV3D,DZ3D,P3D,S3D,W3D,SR,ICENCV,VINCV, & VSNCV,VGNCV,FINCV,FSNCV,AINCV,ASNCV,AGNCV,AHNCV, & I3MNCV,CLODNCV,RAINNC,RAINNCV,SNOWNC,SNOWNCV,GRAPNC, & GRAPNCV,HAILNC,HAILNCV,DT,DTMN,QC3D,QR3D,QI3D,QS3D, & QG3D,QH3D,NC3D,NR3D,NI3D,NS3D,NG3D,NH3D,VI3D,VS3D, & VG3D,FI3D,FS3D,AI3D,AS3D,AG3D,AH3D,I3M3D,QDCN3D, & QTCN3D,QCCN3D,QRCN3D,QNIN3D,KTS,KTE) !====================================================================== IMPLICIT NONE INTEGER, INTENT(IN) :: KTS,KTE INTEGER :: K,A,IV,IV0,IV1,IV2,IV3,IV4,IM,IT,J INTEGER, DIMENSION(KTS:KTE) :: HID REAL, DIMENSION(KTS:KTE), INTENT(IN) :: W3D,P3D REAL, DIMENSION(KTS:KTE) :: TK3D,TC3D,QV3D,DZ3D,S3D,QC3D,QR3D, & QI3D,QS3D,QG3D,QH3D,QDCN3D,QTCN3D,QCCN3D,QRCN3D,QNIN3D, & NC3D,NR3D,NI3D,NS3D,NG3D,NH3D,VI3D,VS3D,VG3D,FI3D,FS3D, & AI3D,AS3D,AG3D,AH3D,I3M3D,TAIR,PRES,QVAP,ADAGR,AMI,BMI,AVI,& BVI,AMS,AVS,BVS,AMG,AVG,BVG,AVH,BVH,INHGR,DCLDMT,DCLDMC, & DCLDMR,DTAIR,DPDT,DQVAP,QACac,QACar,VT_QC,VT_QR,VT_QI, & VT_QS,VT_QG,VT_QH,VT_NC,VT_NR,VT_NI,VT_NS,VT_NG,VT_NH, & VT_VI,VT_VS,VT_VG,VT_FI,VT_FS,VT_AI,VT_AS,VT_AG,VT_AH, & VTI3M,TK0,P40,RHO,CPM,GQCTR,ESW,ESI,QVSW,QVSI,SSRW,SSRI, & XXLV,XXLS,XXLF,XDNC,XDNR,IPF,IPG,I3M0,ZETA,BEST,MVRC,MVRR, & SIGC,MNRC,KDX,BDR,BDI,BDS,BDG,BDH,R3M3D,S3M3D,G3M3D,H3M3D, & FDI,FDS,FDG,FDH,RHOI,RHOS,RHOG,IASPR,SASPR,MVDC,MVDR,MVDI, & MVDS,MVDG,MVDH,EFRC,EFRR,AFAC,AFAR,AFAI,AFAS,AFAG,AFAH, & LAMC,LAMR,LAMI,LAMS,LAMG,LAMH,LTK,LQC,LQR,LQI,LQS,LQG,LQH, & MDI,MDS,MDG REAL :: SR,ICENCV,CLODNCV,RAINNC,RAINNCV,SNOWNC,SNOWNCV,GRAPNC, & GRAPNCV,HAILNC,HAILNCV,DT,DTMN,SDTL,SDTS,DTL,DTS,CNDMAX, & TMP,XMASS,TEMPC,TEMPR,TEMPT,FINCV,FSNCV,VINCV,VSNCV, & VGNCV,AINCV,ASNCV,AGNCV,AHNCV,I3MNCV REAL, DIMENSION(KTS:KTE,3) :: VTMEAN LOGICAL, DIMENSION(KTS:KTE,NAER) :: DO_FR ! re-construct r-cutoff from XA & XC for once only REAL, DIMENSION(KTS:KTE,NAER) :: RX0 ! CUTOFF RADIUS (m) REAL, DIMENSION(KTS:KTE,NAERT) :: AERO ! AEROSOL MIXING RATIO REAL, DIMENSION(MAER) :: ABCD REAL, DIMENSION(NCCN,KTS:KTE,NAER) :: ZCCNS ! AEROSOL # OF EACH MODE (#/kg) !----- AT SUBSATURATION, REMOVE SMALL AMOUNTS OF HYDROMETEORS ---------- DO K = KTS,KTE AERO(K,1) = MAX(0.,QDCN3D(K)) AERO(K,2) = MAX(0.,QTCN3D(K)) AERO(K,3) = MAX(0.,QCCN3D(K)) AERO(K,4) = MAX(0.,QRCN3D(K)) AERO(K,5) = MAX(0.,QNIN3D(K)) TC3D(K) = TK3D(K)-TK0C ESW(K) = MIN(0.99*P3D(K),POLYSVP(TK3D(K),0)) ESI(K) = MIN(0.99*P3D(K),POLYSVP(TK3D(K),1)) IF (ESI(K).GT.ESW(K)) ESI(K) = ESW(K) SSRW(K) = QV3D(K)/(0.622*ESW(K)/(P3D(K)-ESW(K)))-1. SSRI(K) = QV3D(K)/(0.622*ESI(K)/(P3D(K)-ESI(K)))-1. XXLV(K) = 3.1484E6-2370.*TK3D(K) XXLS(K) = 3.15E6-2370.*TK3D(K)+0.3337E6 XXLF(K) = 2836310.8-(3.1484E6-2370.*TK3D(K)) CPM(K) = CP*(1.+0.887*QV3D(K)) HID(K) = MAX(MIN(NINT(ABS(TC3D(K))/0.25),120),0) !------------------------------------------------------------------------- IF (QC3D(K).GE.QSMALL.AND.NC3D(K).LT.NSMALL) THEN LTK(K) = LOG(TK3D(K)) LQC(K) = -1.*LOG(QC3D(K)) NC3D(K) = EXP(DNC0+DNC1*LTK(K)+DNC2*LTK(K)**2.+DNC3* & LTK(K)**3.-0.25*LQC(K)) IF (AFAC_3M.EQ.0) THEN AFAC(K) = AFAC0 ELSEIF (AFAC_3M.EQ.1) THEN MVRC(K) = (QC3D(K)/NC3D(K)/C4PI3W)**THRD MVRC(K) = MIN(MAX(MVRC(K),RCMIN),RCMAX) EFRC(K) = EXP(EFC1+EFC2*LOG(NC3D(K))+EFC3*LOG(MVRC(K))) KDX(K) = MAX(KCCMIN,MIN(KCCMAX,(MVRC(K)/EFRC(K))**3.)) ! KDX(K) = QC3D(K)*V2M3/(8.*EFRC(K)**3.*NC3D(K)*RHOW) KDX(K) = MIN(KCCMAX,MAX(KCCMIN,KDX(K))) AFAC(K) = (6.*KDX(K)-3.+SQRT(8.*KDX(K)+1.))/(2.-2.* & KDX(K)) AFAC(K) = MIN(MAX(AFAC(K),AFAMIN),AFAMAX) ELSEIF (AFAC_3M.EQ.2) THEN SIGC(K) = EXP(SIG1+SIG2*LOG(NC3D(K))+SIG3*LOG(QC3D(K))) MNRC(K) = EXP(MNR1+MNR2*LOG(NC3D(K))+MNR3*LOG(QC3D(K))) AFAC(K) = MIN(MAX(SIGC(K)/MNRC(K),AFAMIN),AFAMAX) ENDIF ENDIF IF (QR3D(K).GE.QSMALL.AND.NR3D(K).LT.NSMALL) THEN LTK(K) = LOG(TK3D(K)) LQR(K) = -1.*LOG(QR3D(K)) NR3D(K) = EXP(-5793.7852+3191.1171*LTK(K)-582.73279* & LTK(K)**2.+35.346854*LTK(K)**3.-0.25*LQR(K)) IF (AFAR_3M.EQ.0) THEN AFAR(K) = AFAR0 ELSEIF (AFAR_3M.EQ.1) THEN MVRR(K) = (QR3D(K)/NR3D(K)/C4PI3W)**THRD MVRR(K) = MIN(RRMAX,MAX(RRMIN,MVRR(K))) EFRR(K) = EXP(EFR1+EFR2*LOG(NR3D(K))+EFR3*LOG(MVRR(K))) KDX(K) = MAX(KCRMIN,MIN(KCRMAX,(MVRR(K)/EFRR(K))**3.)) ! KDX(K) = QR3D(K)*V2M3/(8.*EFRR(K)**3.*NR3D(K)*RHOW) KDX(K) = MIN(KCRMAX,MAX(KCRMIN,KDX(K))) AFAR(K) = (6.*KDX(K)-3.+SQRT(8.*KDX(K)+1.))/(2.-2.* & KDX(K)) AFAR(K) = MIN(MAX(AFAR(K),AFAMIN),AFAMAX) ELSEIF (AFAR_3M.EQ.2) THEN BDR(K) = (QR3D(K)*V2M3/RHOW/NR3D(K))**THRD BDR(K) = MIN(MAX(BDR(K),DRMIN),DRMAX) AFAR(K) = MAX(AFAMIN,19.*TANH(0.6*(1.E3*BDR(K)-1.8))+17.) ENDIF ENDIF IF (QI3D(K).GE.QSMALL.AND.NI3D(K).LT.NSMALL) THEN IF (ICE_RHOI.EQ.0) THEN RHOI(K) = RHOI0 VI3D(K) = 0. ELSEIF (ICE_RHOI.EQ.1) THEN QVSI(K) = 0.622*ESI(K)/(P3D(K)-ESI(K)) INHGR(K) = ITBLE(HID(K)) RHOI(K) = RHOI0*EXP(-3.*MAX((QV3D(K)-QVSI(K))-5.E-5, & 0.)/INHGR(K)) VI3D(K) = QI3D(K)/RHOI(K) ELSEIF (ICE_RHOI.EQ.2) THEN RHOI(K) = RHOI1 VI3D(K) = 0. ENDIF LTK(K) = LOG(TK3D(K)) LQI(K) = -1.*LOG(QI3D(K)) MDI(K) = EXP(-3.2653646+2.0539073*LTK(K)-0.25*LQI(K))/1.E3 NI3D(K) = 1.E9*QI3D(K)*V2M3/RHOI(K)/MDI(K)**3. I3M3D(K) = QI3D(K)*V2M3/RHOI(K) IF (ICE_SHAPE.EQ.0) THEN FI3D(K) = I3M3D(K) ELSEIF (ICE_SHAPE.EQ.1) THEN ADAGR(K) = (MAX(MIN(ITBLE(HID(K)),2.),0.5))**THRD ZETA(K) = (ADAGR(K)-1.)/(ADAGR(K)+2.) I3M0(K) = NI3D(K)*DI0**3. FI3D(K) = (I3M3D(K)/I3M0(K))**ZETA(K)*I3M3D(K) ENDIF IF (AFAI_3M.EQ.0) THEN AFAI(K) = AFAI0 AI3D(K) = 0. ELSE BDI(K) = (I3M3D(K)/NI3D(K))**THRD*1.E3 BDI(K) = MIN(MAX(BDI(K),DIMIN*1.E3),DIMAX*1.E3) FDI(K) = 7.4015986E-2+0.79866676*BDI(K)-9.4468892E-3* & LOG(NI3D(K))+0.38235092*BDI(K)**2.+ & 2.9811542E-4*LOG(NI3D(K))**2.+1.9052614E-2* & BDI(K)*LOG(NI3D(K)) KDX(K) = MAX(KCIMIN,MIN(KCIMAX,(BDI(K)/FDI(K))**3.)) AFAI(K) = (6.*KDX(K)-3.+SQRT(8.*KDX(K)+1.))/(2.-2.* & KDX(K)) AFAI(K) = MIN(MAX(AFAI(K),AFAMIN),AFAMAX) ! AFAI(K) = MAX(AFAMIN,12.*TANH(0.7*(BDI(K)-1.7))+11.) IF (AFAI_3M.EQ.1) THEN KDX(K) = (AFAI(K)**2.+3.*AFAI(K)+2.)/(AFAI(K)**2.+ & 6.*AFAI(K)+9.) AI3D(K) = (KDX(K)*NI3D(K)*I3M3D(K)**2.)**THRD ELSEIF (AFAI_3M.EQ.2) THEN AI3D(K) = 0. ENDIF ENDIF ENDIF IF (QS3D(K).GE.QSMALL.AND.NS3D(K).LT.NSMALL) THEN IF (ICE_RHOS.EQ.0) THEN RHOS(K) = RHOS0 VS3D(K) = 0. ELSE LTK(K) = LOG(TK3D(K)) LQS(K) = -1.*LOG(QS3D(K)) IF (TK3D(K).LT.TK0C) THEN RHOS(K) = 15740.702-6098.0087*LTK(K)+503.33089* & LQS(K)+594.29913*LTK(K)**2.+1.9033961* & LQS(K)**2.-94.950429*LTK(K)*LQS(K) ELSE RHOS(K) = EXP(-64808.666+23113.508*LTK(K)-36.46632* & LQS(K)-2060.6024*LTK(K)**2.-0.005729458* & LQS(K)**2.+6.5057411*LTK(K)*LQS(K)) ENDIF RHOS(K) = MIN(RHOIMAX,MAX(RHOIMIN,RHOS(K))) IF (ICE_RHOS.EQ.1) THEN VS3D(K) = QS3D(K)/RHOS(K) ELSEIF (ICE_RHOS.EQ.2) THEN VS3D(K) = 0. ENDIF ENDIF LTK(K) = LOG(TK3D(K)) LQS(K) = -1.*LOG(QS3D(K)) MDS(K) = EXP(-123.23898+40.74706*LTK(K)-3.0333477* & LTK(K)**2.-0.31219981*LQS(K)+0.0012798222* & LQS(K)**2.)/1.E3 NS3D(K) = 1.E9*QS3D(K)*V2M3/RHOS(K)/MDS(K)**3. S3M3D(K) = QS3D(K)*V2M3/RHOS(K) IF (AGG_SHAPE.EQ.0) THEN SASPR(K) = 1. FS3D(K) = 0. ELSEIF (AGG_SHAPE.EQ.1) THEN SASPR(K) = 0.7 FS3D(K) = SASPR(K)*S3M3D(K) ENDIF IF (AFAS_3M.EQ.0) THEN AFAS(K) = AFAS0 AS3D(K) = 0. ELSE BDS(K) = (S3M3D(K)/NS3D(K))**THRD*1.E3 BDS(K) = MIN(MAX(BDS(K),DSMIN*1.E3),DSMAX*1.E3) IF (TK3D(K).GE.TK0C) THEN FDS(K) = -0.21911541+1.2739845*BDS(K)+0.10141003* & LOG(NS3D(K))+0.30063818*BDS(K)**2.- & 4.3857765E-3*LOG(NS3D(K))**2.-0.078801732* & BDS(K)*LOG(NS3D(K)) ELSE IF (QC3D(K).GE.1.E-8) THEN FDS(K) = -1.1527014+2.9067645*BDS(K)+0.25316062* & LOG(NS3D(K))-0.17768557*BDS(K)**2.- & 0.013117292*LOG(NS3D(K))**2.-0.17020429* & BDS(K)*LOG(NS3D(K)) ELSE FDS(K) = -0.2813929+1.7275463*BDS(K)+0.045550156* & LOG(NS3D(K))-0.16526226*BDS(K)**2.- & 1.7699916E-3*LOG(NS3D(K))**2.- & 4.6441257E-2*BDS(K)*LOG(NS3D(K)) ENDIF ENDIF KDX(K) = MAX(KCSMIN,MIN(KCSMAX,(BDS(K)/FDS(K))**3.)) AFAS(K) = (6.*KDX(K)-3.+SQRT(8.*KDX(K)+1.))/(2.-2.* & KDX(K)) AFAS(K) = MIN(MAX(AFAS(K),AFAMIN),AFAMAX) ! AFAS(K) = MAX(AFAMIN,4.5*TANH(0.5*(BDS(K)-5.))+5.5) IF (AFAS_3M.EQ.1) THEN KDX(K) = (AFAS(K)**2.+3.*AFAS(K)+2.)/(AFAS(K)**2.+ & 6.*AFAS(K)+9.) AS3D(K) = (KDX(K)*NS3D(K)*S3M3D(K)**2.)**THRD ELSEIF (AFAS_3M.EQ.2) THEN AS3D = 0. ENDIF ENDIF ENDIF IF (QG3D(K).GE.QSMALL.AND.NG3D(K).LT.NSMALL) THEN IF (ICE_RHOG.EQ.0) THEN RHOG(K) = RHOG1 VG3D(K) = 0. ELSEIF (ICE_RHOG.EQ.1.OR.ICE_RHOG.EQ.2) THEN LTK(K) = LOG(TK3D(K)) LQG(K) = -1.*LOG(QG3D(K)) IF (TK3D(K).LT.TK0C) THEN RHOG(K) = EXP(-509.23219+187.98419*LTK(K)+1.5493489* & LQG(K)-17.147561*LTK(K)**2.-0.005829723* & LQG(K)**2.-0.26664109*LTK(K)*LQG(K)) ELSE RHOG(K) = EXP(-2894.6795+1050.5173*LTK(K)-11.536497* & LQG(K)-95.054373*LTK(K)**2.-0.001986964* & LQG(K)**2.+2.0502147*LTK(K)*LQG(K)) ENDIF VG3D(K) = QG3D(K)/RHOG(K) ENDIF LTK(K) = LOG(TK3D(K)) LQG(K) = -1.*LOG(QG3D(K)) MDG(K) = EXP(-2205.8027+1225.8046*LTK(K)-226.27995* & LTK(K)**2.+13.929644*LTK(K)**3.-0.25*LQG(K))/1.E3 NG3D(K) = 1.E9*QG3D(K)*V2M3/RHOG(K)/MDG(K)**3. G3M3D(K) = QG3D(K)*V2M3/RHOG(K) IF (AFAG_3M.EQ.0) THEN AFAG(K) = AFAG0 AG3D(K) = 0. ELSE BDG(K) = (G3M3D(K)/NG3D(K))**THRD*1.E3 BDG(K) = MIN(MAX(BDG(K),DGMIN*1.E3),DGMAX*1.E3) IF (TK3D(K).GE.TK0C) THEN FDG(K) = 0.58006354+0.79661229*BDG(K)-0.18394382* & LOG(NG3D(K))+6.7371044E-2*BDG(K)**2.+ & 9.832945E-3*LOG(NG3D(K))**2.+0.12433055* & BDG(K)*LOG(NG3D(K)) ELSE IF (QC3D(K).GE.1.E-8) THEN FDG(K) = 0.17363469+1.5044291*BDG(K)-0.050639722* & LOG(NG3D(K))+1.5101052E-2*BDG(K)**2.+ & 2.5974719E-3*LOG(NG3D(K))**2.+0.01961464*& BDG(K)*LOG(NG3D(K)) ELSE FDG(K) = 0.59259317-0.89933515*BDG(K)+2.0893032* & BDG(K)**2.-0.50305755*BDG(K)**3.- & 2.2446793E-2*LOG(NG3D(K))+2.7589047E-3* & LOG(NG3D(K))**2. ENDIF ENDIF KDX(K) = MAX(KCGMIN,MIN(KCGMAX,(BDG(K)/FDG(K))**3.)) AFAG(K) = (6.*KDX(K)-3.+SQRT(8.*KDX(K)+1.))/(2.-2.* & KDX(K)) AFAG(K) = MIN(MAX(AFAG(K),AFAMIN),AFAMAX) ! AFAG(K) = MAX(AFAMIN,5.5*TANH(0.7*(BDG(K)-4.5))+8.5) IF (AFAG_3M.EQ.1) THEN KDX(K) = (AFAG(K)**2.+3.*AFAG(K)+2.)/(AFAG(K)**2.+ & 6.*AFAG(K)+9.) AG3D(K) = (KDX(K)*NG3D(K)*G3M3D(K)**2.)**THRD ELSEIF (AFAG_3M.EQ.2) THEN AG3D = 0. ENDIF ENDIF ENDIF IF (QH3D(K).GE.QSMALL.AND.NH3D(K).LT.NSMALL) THEN LTK(K) = LOG(TK3D(K)) LQH(K) = -1.*LOG(QH3D(K)) NH3D(K) = EXP(22.929406-4.2328364*LTK(K)+0.30647567* & LTK(K)**2.-0.009233271*LTK(K)**3.-0.25*LQH(K)) H3M3D(K) = QH3D(K)*V2M3/RHOH IF (AFAH_3M.EQ.0) THEN AFAH(K) = AFAH0 AH3D(K) = 0. ELSE BDH(K) = (H3M3D(K)/NH3D(K))**THRD*1.E3 BDH(K) = MIN(MAX(BDH(K),DHMIN*1.E3),DHMAX*1.E3) IF (TK3D(K).GE.TK0C) THEN FDH(K) = 1.157754+0.37852874*BDH(K)-0.11129737* & LOG(NH3D(K))+0.13929599*BDH(K)**2.+ & 8.1105237E-3*LOG(NH3D(K))**2.+5.7432113E-2* & BDH(K)*LOG(NH3D(K)) ELSE FDH(K) = -0.48246793+2.0407077*BDH(K)+0.022262969* & LOG(NH3D(K))-0.158389*BDH(K)**2.- & 5.5545804E-3*LOG(NH3D(K))**2.+2.9443577E-2* & BDH(K)*LOG(NH3D(K)) ENDIF KDX(K) = MAX(KCHMIN,MIN(KCHMAX,(BDH(K)/FDH(K))**3.)) AFAH(K) = (6.*KDX(K)-3.+SQRT(8.*KDX(K)+1.))/(2.-2.* & KDX(K)) AFAH(K) = MIN(MAX(AFAH(K),AFAMIN),AFAMAX) ! AFAH(K) = MAX(AFAMIN,3.7*TANH(0.3*(BDH(K)-9.))+6.5) IF (AFAH_3M.EQ.1) THEN KDX(K) = (AFAH(K)**2.+3.*AFAH(K)+2.)/(AFAH(K)**2.+ & 6.*AFAH(K)+9.) AH3D(K) = (KDX(K)*NH3D(K)*H3M3D(K)**2.)**THRD ELSEIF (AFAH_3M.EQ.2) THEN AH3D(K) = 0. ENDIF ENDIF ENDIF !----------------------------------------------------------------------------- IF (SSRW(K).LT.-1.E-1) THEN IF (QR3D(K).LT.QSMAL1) THEN QV3D(K) = QV3D(K)+QR3D(K) TK3D(K) = TK3D(K)-QR3D(K)*XXLV(K)/CPM(K) QR3D(K) = 0.; NR3D(K) = 0. ENDIF IF (QC3D(K).LT.QSMAL1) THEN QV3D(K) = QV3D(K)+QC3D(K) TK3D(K) = TK3D(K)-QC3D(K)*XXLV(K)/CPM(K) QC3D(K) = 0.; NC3D(K) = 0. ENDIF ENDIF IF (SSRI(K).LT.-1.E-1) THEN IF (QI3D(K).LT.QSMAL1) THEN QV3D(K) = QV3D(K)+QI3D(K) TK3D(K) = TK3D(K)-QI3D(K)*XXLS(K)/CPM(K) QI3D(K) = 0.; NI3D(K) = 0.; I3M3D(K) = 0.; FI3D(K) = 0. VI3D(K) = 0.; AI3D(K) = 0. ENDIF IF (QS3D(K).LT.QSMAL1) THEN QV3D(K) = QV3D(K)+QS3D(K) TK3D(K) = TK3D(K)-QS3D(K)*XXLS(K)/CPM(K) QS3D(K) = 0.; NS3D(K) = 0.; VS3D(K) = 0.; AS3D(K) = 0. FS3D(K) = 0. ENDIF IF (QG3D(K).LT.QSMAL1) THEN QV3D(K) = QV3D(K)+QG3D(K) TK3D(K) = TK3D(K)-QG3D(K)*XXLS(K)/CPM(K) QG3D(K) = 0.; NG3D(K) = 0.; VG3D(K) = 0.; AG3D(K) = 0. ENDIF IF (QH3D(K).LT.QSMAL1) THEN QV3D(K) = QV3D(K)+QH3D(K) TK3D(K) = TK3D(K)-QH3D(K)*XXLS(K)/CPM(K) QH3D(K) = 0.; NH3D(K) = 0.; AH3D(K) = 0. ENDIF ENDIF !----- MELT VERY SMALL SNOW/GRAUPEL/HAIL MIXING RATIOS, ADD TO RAIN ---- IF (TK3D(K).GE.TK0C.AND.QS3D(K).LT.QLIMIT) THEN QR3D(K) = QR3D(K)+QS3D(K) NR3D(K) = NR3D(K)+NS3D(K) TK3D(K) = TK3D(K)-QS3D(K)*XXLF(K)/CPM(K) QS3D(K) = 0.; NS3D(K) = 0.; VS3D(K) = 0.; AS3D(K) = 0. FS3D(K) = 0. ENDIF IF (TK3D(K).GE.TK0C.AND.QG3D(K).LT.QLIMIT) THEN QR3D(K) = QR3D(K)+QG3D(K) NR3D(K) = NR3D(K)+NG3D(K) TK3D(K) = TK3D(K)-QG3D(K)*XXLF(K)/CPM(K) QG3D(K) = 0.; NG3D(K) = 0.; VG3D(K) = 0.; AG3D(K) = 0. ENDIF IF (TK3D(K).GE.TK0C.AND.QH3D(K).LT.QLIMIT) THEN QR3D(K) = QR3D(K)+QH3D(K) NR3D(K) = NR3D(K)+NH3D(K) TK3D(K) = TK3D(K)-QH3D(K)*XXLF(K)/CPM(K) QH3D(K) = 0.; NH3D(K) = 0.; AH3D(K) = 0. ENDIF ENDDO !------- INITIALIZE VARIABLE,ONLY T&QV ARE TREATED DIFFERENTLY --------- DO K = KTS,KTE TAIR(K) = TK3D(K) QVAP(K) = QV3D(K) DTAIR(K) = -W3D(K)*G/CP ! Lagrangian parcel dry adiabatic lapse rate DQVAP(K) = 0. PRES(K) = P3D(K) TK0(K) = TAIR(K)-DTAIR(K)*DT ! diabatically decend parcel P40(K) = PRES(K)*(TK0(K)/TAIR(K))**(CP/R) TAIR(K) = TK0(K) IV1 = 1 DO IV = 1,NAER IV2 = IV1+1 DO_FR(K,IV) = .TRUE. ! FIND CUT OFF RADIUS RX0(K,IV) = 99. ! initial R cutoff IF (NAERN(IV).GE.2) THEN ! have info of dry and total aerosol AERO(K,IV1) = MIN(AERO(K,IV1),AERO(K,IV2)) ! dry aerosol < total aerosol DO IM = 1,NCCN ! calculate CCN # [#/kg] from CCN mass ZCCNS(IM,K,IV) = AERO(K,IV2)*WMAS(IM,IV)*RFACT(IM,IV) ENDDO ELSE ! only info available is total aerosol DO IM = 1,NCCN ! calculate CCN # [#/kg] from CCN mass ZCCNS(IM,K,IV) = AERO(K,IV1)*WMAS(IM,IV)*RFACT(IM,IV) ENDDO ENDIF IV1 = IV1+NAERN(IV) ENDDO ENDDO ! INITILIZATION LOOPS SDTL = 0. ! initial integrated largedt timestep DO IT = 1,INT((DT-DT20S*0.01)/DT20S)+1 DTL = MIN(DT20S,DT-SDTL) ! initial largedt timestep DO K = KTS,KTE DTS = DTL ! initial smalldt timestep (could be changed later) SDTS = 0. ! initial integrated smalldt timestep TMP = W3D(K)*G/CP/TK0(K) ! 1-TAIR/TK0 XDNC(K) = 0.; XDNR(K) = 0.; GQCTR(K) = 0.; ADAGR(K) = 1. RHOI(K) = RHOI0; RHOS(K) = RHOS0; RHOG(K) = RHOG0 AFAC(K) = AFAC0; AFAR(K) = AFAR0; AFAI(K) = AFAI0 AFAS(K) = AFAS0; AFAG(K) = AFAG0; AFAH(K) = AFAH0 AVI(K) = AVI0; BVI(K) = BVI0; AVS(K) = AVS0; BVS(K) = BVS0 AVG(K) = AVG0; BVG(K) = BVG0; AVH(K) = AVH0; BVH(K) = BVH0 AMI(K) = AMI0; BMI(K) = BMI0; AMS(K) = AMS0; AMG(K) = AMG0 IASPR(K) = 1.; SASPR(K) = 1.; QACac(K) = 0.; QACar(K) = 0. 111 CONTINUE PRES(K) = P40(K)*(1.-TMP*(SDTL+SDTS))**(CP/R) DPDT(K) = -G*W3D(K)*PRES(K)/(R*TK0(K)*(1.-TMP*(SDTL+SDTS))) ESW(K) = MIN(0.99*PRES(K),POLYSVP(TAIR(K),0)) ESI(K) = MIN(0.99*PRES(K),POLYSVP(TAIR(K),1)) IF (ESI(K).GT.ESW(K)) ESI(K) = ESW(K) QVSW(K) = 0.622*ESW(K)/(PRES(K)-ESW(K)) QVSI(K) = 0.622*ESI(K)/(PRES(K)-ESI(K)) SSRW(K) = QVAP(K)/QVSW(K)-1. SSRI(K) = QVAP(K)/QVSI(K)-1. XXLV(K) = 3.1484E6-2370.*TAIR(K) CPM(K) = CP*(1.+0.887*QVAP(K)) RHO(K) = PRES(K)/(TAIR(K)*(1.+0.61*QVAP(K)))/R ! [KG/M**3] CNDMAX = (QVAP(K)-QVSW(K))/(1.+XXLV(K)**2.*QVAP(K)/ & (CPM(K)*RV*TAIR(K)**2.)) !--------------------------- CCN ACTIVATION. --------------------------- IF (SSRW(K).GT.RSMALL.AND.SSRW(K).GT.S3D(K)) THEN DCLDMT(K) = 0.; DCLDMC(K) = 0.; DCLDMR(K) = 0. IV1 = 1 DO IV = 1,NAER DO J = 1,NAERN(IV) ABCD(J) = AERO(K,IV1-1+J) ENDDO IF (AERO(K,IV1).GT.1.E-20.AND.NAERN(IV).GE.2.AND. & W3D(K).GT.1.E-5) THEN CALL ACTIVA(TAIR(K),W3D(K),NC3D(K),NR3D(K),ABCD, & QACac(K),QACar(K),RX0(K,IV),ZCCNS(1,K,IV), & SSRW(K),IV,DO_FR(K,IV)) DCLDMC(K) = DCLDMC(K)+QACac(K) DCLDMR(K) = DCLDMR(K)+QACar(K) DO J = 1,NAERN(IV) AERO(K,IV1-1+J) = ABCD(J) ENDDO ENDIF IV1 = IV1+NAERN(IV) ENDDO DCLDMT(K) = DCLDMC(K)+DCLDMR(K) IF (DCLDMT(K).GT.RLIMIT) THEN IF (DCLDMT(K).GT.CNDMAX) THEN DCLDMC(K) = DCLDMC(K)*(CNDMAX/DCLDMT(K)) DCLDMR(K) = DCLDMR(K)*(CNDMAX/DCLDMT(K)) DCLDMT(K) = DCLDMC(K)+DCLDMR(K) ENDIF QVAP(K) = MAX(0.,QVAP(K)-DCLDMT(K)) QC3D(K) = MAX(0.,QC3D(K)+DCLDMC(K)) QR3D(K) = MAX(0.,QR3D(K)+DCLDMR(K)) CPM(K) = CP*(1.+0.887*QVAP(K)) XXLV(K) = 3.1484E6-2370.*TAIR(K) TAIR(K) = TAIR(K)+XXLV(K)*DCLDMT(K)/CPM(K) ENDIF ESW(K) = MIN(0.99*PRES(K),POLYSVP(TAIR(K),0)) ESI(K) = MIN(0.99*PRES(K),POLYSVP(TAIR(K),1)) IF (ESI(K).GT.ESW(K)) ESI(K) = ESW(K) SSRW(K) = QVAP(K)/(0.622*ESW(K)/(PRES(K)-ESW(K)))-1. SSRI(K) = QVAP(K)/(0.622*ESI(K)/(PRES(K)-ESI(K)))-1. RHO(K) = PRES(K)/(TAIR(K)*(1.+0.61*QVAP(K)))/R ENDIF S3D(K) = MAX(S3D(K),SSRW(K)) !------------------------- ICE NUCLEATION ------------------------------ IF (TAIR(K).LT.(TK0C-2.).AND.SSRI(K).GT.1.E-4) THEN CALL ICENU(TAIR(K),PRES(K),DTS,RHO(K),QVAP(K),QI3D(K), & NI3D(K),VI3D(K),FI3D(K),AI3D(K),I3M3D(K),AERO(K,5)) ESW(K) = MIN(0.99*PRES(K),POLYSVP(TAIR(K),0)) ESI(K) = MIN(0.99*PRES(K),POLYSVP(TAIR(K),1)) IF (ESI(K).GT.ESW(K)) ESI(K) = ESW(K) SSRW(K) = QVAP(K)/(0.622*ESW(K)/(PRES(K)-ESW(K)))-1. SSRI(K) = QVAP(K)/(0.622*ESI(K)/(PRES(K)-ESI(K)))-1. RHO(K) = PRES(K)/(TAIR(K)*(1.+0.61*QVAP(K)))/R ENDIF !-----------MICROPHYSICAL TENDENCY CALCULATION IN SMALL TIMESTEP ------- IF (QC3D(K).GE.QSMALL.OR.QR3D(K).GE.QSMALL.OR.QI3D(K).GE. & QSMALL.OR.QS3D(K).GE.QSMALL.OR.QG3D(K).GE.QSMALL.OR. & QH3D(K).GE.QSMALL) THEN CALL SMALL_DT(DTL,DTS,SDTS,DTAIR(K),DQVAP(K),DPDT(K), & RHO(K),TAIR(K),PRES(K),QVAP(K),QC3D(K),QR3D(K), & QI3D(K),QS3D(K),QG3D(K),QH3D(K),NC3D(K),NR3D(K), & NI3D(K),NS3D(K),NG3D(K),NH3D(K),VI3D(K),VS3D(K), & VG3D(K),FI3D(K),FS3D(K),AI3D(K),AS3D(K),AG3D(K), & AH3D(K),I3M3D(K),SASPR(K),XDNC(K),XDNR(K),GQCTR(K)) ENDIF !----------------- UPDATE TIME & PROGNOSTIC VARIABLES ------------------ TAIR(K) = TAIR(K)+DTAIR(K)*DTS QVAP(K) = QVAP(K)+DQVAP(K)*DTS SDTS = SDTS+DTS IF (SDTS.GE.DTL) GOTO 333 GOTO 111 333 CONTINUE !----------------- RAIN DROP DEACTIVATION. ----------------------------- IF (XDNR(K).GT.1.) THEN IV0 = 0 DO IV = 1,NAER IF (NAERN(IV).EQ.4) THEN IV3 = IV0+3 IV4 = IV0+4 IF (NR3D(K).GT.1.) THEN XMASS = AERO(K,IV4)*MIN(1.,XDNR(K)/NR3D(K)) ELSE XMASS = AERO(K,IV4) ENDIF AERO(K,IV3) = AERO(K,IV3)+XMASS AERO(K,IV4) = AERO(K,IV4)-XMASS ENDIF IV0 = IV0+NAERN(IV) ENDDO ENDIF !----------------- CLOUD DROP DEACTIVATION. ---------------------------- IF (XDNC(K).GT.1.) THEN IV1 = 1 DO IV = 1,NAER IV2 = IV1+1 IV3 = IV2+1 IF (NAERN(IV).GE.3) THEN IF (NC3D(K).GT.1.) THEN XMASS = AERO(K,IV3)*MIN(1.,XDNC(K)/NC3D(K)) ELSE XMASS = AERO(K,IV3) ENDIF AERO(K,IV1) = AERO(K,IV1)+XMASS AERO(K,IV3) = AERO(K,IV3)-XMASS AERO(K,IV1) = MIN(AERO(K,IV1),AERO(K,IV2)) ELSEIF (NAERN(IV).EQ.2) THEN CALL DEACTIVA(XDNC(K),QC3D(K),NC3D(K),AERO(K,IV1),& AERO(K,IV2),AERO(K,IV3),RX0(K,IV), & ZCCNS(1,K,IV),IV,DO_FR(K,IV)) ELSE XMASS = AERO(K,IV1)*MIN(1.,XDNC(K)/ & (ZCCNS(2,K,IV)+ZCCNS(3,K,IV))) ! COARSE AND ACCUMULATION MODE AERO(K,IV1) = AERO(K,IV1)+XMASS ENDIF IV1 = IV1+NAERN(IV) ENDDO ENDIF PRES(K) = P40(K)*(1.-TMP*(SDTL+SDTS))**(CP/R) RHO(K) = PRES(K)/(TAIR(K)*(1.+0.61*QVAP(K)))/R ENDDO !-----------MICROPHYSICAL TENDENCY CALCULATION IN LARGE TIMESTEP ------- DO K = KTS,KTE IF (QC3D(K).GE.QSMALL.OR.QR3D(K).GE.QSMALL.OR.QI3D(K).GE. & QSMALL.OR.QS3D(K).GE.QSMALL.OR.QG3D(K).GE.QSMALL.OR. & QH3D(K).GE.QSMALL) THEN CALL LARGE_DT(DTL,TAIR(K),QVAP(K),PRES(K),RHO(K), & QC3D(K),QR3D(K),QI3D(K),QS3D(K),QG3D(K),QH3D(K), & NC3D(K),NR3D(K),NI3D(K),NS3D(K),NG3D(K),NH3D(K), & VI3D(K),VS3D(K),VG3D(K),FI3D(K),FS3D(K),AI3D(K), & AS3D(K),AG3D(K),AH3D(K),I3M3D(K),SASPR(K),GQCTR(K)) ENDIF IV0 = 0 DO IV = 1,NAER ! for all aerosol compoment IF (NAERN(IV).GE.4) THEN ! with aerosol in cloud/precip. recored IV3 = IV0+3 ! aerosols in clouds ll IV4 = IV0+4 ! aerosols in rains IF (GQCTR(K).GT.0.) THEN ! + from cloud to precip. XMASS = GQCTR(K)*AERO(K,IV3) ! + from cloud to precip. ELSE ! - from precip. to cloud XMASS = GQCTR(K)*AERO(K,IV4) ! - from precip. to cloud ENDIF AERO(K,IV3) = AERO(K,IV3)-XMASS ! update aerosol in cloud AERO(K,IV4) = AERO(K,IV4)+XMASS ! update aerosol in precipi. ENDIF IV0 = IV0+NAERN(IV) ENDDO ENDDO SDTL = SDTL+DTL ENDDO ! FOR IT LOOPS !----- AT SUBSATURATION, REMOVE SMALL AMOUNTS OF HYDROMETEORS ---------- DO K = KTS,KTE ESW(K) = MIN(0.99*PRES(K),POLYSVP(TAIR(K),0)) ESI(K) = MIN(0.99*PRES(K),POLYSVP(TAIR(K),1)) IF (ESI(K).GT.ESW(K)) ESI(K) = ESW(K) SSRW(K) = QVAP(K)/(0.622*ESW(K)/(PRES(K)-ESW(K)))-1. SSRI(K) = QVAP(K)/(0.622*ESI(K)/(PRES(K)-ESI(K)))-1. XXLV(K) = 3.1484E6-2370.*TAIR(K) XXLS(K) = 3.15E6-2.37E3*TAIR(K)+0.3337E6 CPM(K) = CP*(1.+0.887*QVAP(K)) IF (SSRW(K).LT.-1.E-1) THEN IF (QR3D(K).LT.QSMAL1) THEN QVAP(K) = QVAP(K)+QR3D(K) TAIR(K) = TAIR(K)-QR3D(K)*XXLV(K)/CPM(K) QR3D(K) = 0.; NR3D(K) = 0. ENDIF IF (QC3D(K).LT.QSMAL1) THEN QVAP(K) = QVAP(K)+QC3D(K) TAIR(K) = TAIR(K)-QC3D(K)*XXLV(K)/CPM(K) QC3D(K) = 0.; NC3D(K) = 0. ENDIF ENDIF IF (SSRI(K).LT.-1.E-1) THEN IF (QI3D(K).LT.QSMAL1) THEN QVAP(K) = QVAP(K)+QI3D(K) TAIR(K) = TAIR(K)-QI3D(K)*XXLS(K)/CPM(K) QI3D(K) = 0.; NI3D(K) = 0.; I3M3D(K) = 0.; FI3D(K) = 0. VI3D(K) = 0.; AI3D(K) = 0. ENDIF IF (QS3D(K).LT.QSMAL1) THEN QVAP(K) = QVAP(K)+QS3D(K) TAIR(K) = TAIR(K)-QS3D(K)*XXLS(K)/CPM(K) QS3D(K) = 0.; NS3D(K) = 0.; VS3D(K) = 0.; AS3D(K) = 0. FS3D(K) = 0. ENDIF IF (QG3D(K).LT.QSMAL1) THEN QVAP(K) = QVAP(K)+QG3D(K) TAIR(K) = TAIR(K)-QG3D(K)*XXLS(K)/CPM(K) QG3D(K) = 0.; NG3D(K) = 0.; VG3D(K) = 0.; AG3D(K) = 0. ENDIF IF (QH3D(K).LT.QSMAL1) THEN QVAP(K) = QVAP(K)+QH3D(K) TAIR(K) = TAIR(K)-QH3D(K)*XXLS(K)/CPM(K) QH3D(K) = 0.; NH3D(K) = 0.; AH3D(K) = 0. ENDIF ENDIF IF (TAIR(K).GT.TK0C.AND.QI3D(K).GT.0.) THEN QR3D(K) = QR3D(K)+QI3D(K) TAIR(K) = TAIR(K)-QI3D(K)*XXLF(K)/CPM(K) NR3D(K) = NR3D(K)+NI3D(K) QI3D(K) = 0.; NI3D(K) = 0.; I3M3D(K) = 0.; FI3D(K) = 0. VI3D(K) = 0.; AI3D(K) = 0. ENDIF IF (QC3D(K).LT.QSMALL.OR.NC3D(K).LT.NSMALL) THEN QVAP(K) = QVAP(K)+QC3D(K) TAIR(K) = TAIR(K)-QC3D(K)*XXLV(K)/CPM(K) QC3D(K) = 0.; NC3D(K) = 0. ENDIF IF (QR3D(K).LT.QSMALL.OR.NR3D(K).LT.NSMALL) THEN QVAP(K) = QVAP(K)+QR3D(K) TAIR(K) = TAIR(K)-QR3D(K)*XXLV(K)/CPM(K) QR3D(K) = 0.; NR3D(K) = 0. ENDIF IF (QI3D(K).LT.QSMALL.OR.NI3D(K).LT.NSMALL) THEN QVAP(K) = QVAP(K)+QI3D(K) TAIR(K) = TAIR(K)-QI3D(K)*XXLS(K)/CPM(K) QI3D(K) = 0.; NI3D(K) = 0.; I3M3D(K) = 0.; FI3D(K) = 0. VI3D(K) = 0.; AI3D(K) = 0. ENDIF IF (QS3D(K).LT.QSMALL.OR.NS3D(K).LT.NSMALL) THEN QVAP(K) = QVAP(K)+QS3D(K) TAIR(K) = TAIR(K)-QS3D(K)*XXLS(K)/CPM(K) QS3D(K) = 0.; NS3D(K) = 0.; VS3D(K) = 0.; AS3D(K) = 0. FS3D(K) = 0. ENDIF IF (QG3D(K).LT.QSMALL.OR.NG3D(K).LT.NSMALL) THEN QVAP(K) = QVAP(K)+QG3D(K) TAIR(K) = TAIR(K)-QG3D(K)*XXLS(K)/CPM(K) QG3D(K) = 0.; NG3D(K) = 0.; VG3D(K) = 0.; AG3D(K) = 0. ENDIF IF (QH3D(K).LT.QSMALL.OR.NH3D(K).LT.NSMALL) THEN QVAP(K) = QVAP(K)+QH3D(K) TAIR(K) = TAIR(K)-QH3D(K)*XXLS(K)/CPM(K) QH3D(K) = 0.; NH3D(K) = 0.; AH3D(K) = 0. ENDIF RHO(K) = PRES(K)/(TAIR(K)*(1.+0.61*QVAP(K)))/R VT_QC(K) = 0.; VT_NC(K) = 0.; VT_QR(K) = 0.; VT_NR(K) = 0. VT_QI(K) = 0.; VT_NI(K) = 0.; VT_QS(K) = 0.; VT_NS(K) = 0. VT_QG(K) = 0.; VT_NG(K) = 0.; VT_QH(K) = 0.; VT_NH(K) = 0. VT_VI(K) = 0.; VT_VS(K) = 0.; VT_VG(K) = 0.; VT_FI(K) = 0. VT_FS(K) = 0.; VT_AI(K) = 0.; VT_AS(K) = 0.; VT_AG(K) = 0. VT_AH(K) = 0.; VTI3M(K) = 0. IF (QC3D(K).GE.QSMALL.OR.QR3D(K).GE.QSMALL.OR.QI3D(K).GE. & QSMALL.OR.QS3D(K).GE.QSMALL.OR.QG3D(K).GE.QSMALL.OR. & QH3D(K).GE.QSMALL) THEN CALL SEDI_FALL(TAIR(K),PRES(K),QVAP(K),QC3D(K),QR3D(K), & QI3D(K),QS3D(K),QG3D(K),QH3D(K),NC3D(K),NR3D(K), & NI3D(K),NS3D(K),NG3D(K),NH3D(K),VI3D(K),VS3D(K), & VG3D(K),FI3D(K),FS3D(K),AI3D(K),AS3D(K),AG3D(K), & AH3D(K),I3M3D(K),VT_QC(K),VT_QR(K),VT_QI(K),VT_QS(K), & VT_QG(K),VT_QH(K),VT_NC(K),VT_NR(K),VT_NI(K),VT_NS(K),& VT_NG(K),VT_NH(K),VT_VI(K),VT_VS(K),VT_VG(K),VT_FI(K),& VT_FS(K),VT_AI(K),VT_AS(K),VT_AG(K),VT_AH(K),VTI3M(K),& SASPR(K),RHO(K)) ENDIF IF (QC3D(K).GE.QSMALL.AND.NC3D(K).GE.NSMALL) THEN MVDC(K) = (QC3D(K)*iAMW/NC3D(K))**THRD IF (MVDC(K).GT.DCR) THEN QR3D(K) = QR3D(K)+QC3D(K) NR3D(K) = NR3D(K)+NC3D(K) QC3D(K) = 0.; NC3D(K) = 0. ENDIF ENDIF IF (QR3D(K).GE.QSMALL.AND.NR3D(K).GE.NSMALL) THEN MVDR(K) = (QR3D(K)*iAMW/NR3D(K))**THRD IF (MVDR(K).LT.DCR) THEN QC3D(K) = QC3D(K)+QR3D(K) NC3D(K) = NC3D(K)+NR3D(K) QR3D(K) = 0.; NR3D(K) = 0. ENDIF ENDIF IF (QH3D(K).GE.QSMALL.AND.NH3D(K).GE.NSMALL) THEN MVDH(K) = (QH3D(K)*iAMH/NH3D(K))**THRD IF (MVDH(K).LT.DHMIN) THEN QG3D(K) = QG3D(K)+QH3D(K) NG3D(K) = NG3D(K)+NH3D(K) VG3D(K) = VG3D(K)+QH3D(K)/RHOG0 QH3D(K) = 0.; NH3D(K) = 0. IF (AH3D(K).GE.ASMALL.AND.AFAG_3M.EQ.1) THEN AG3D(K) = AG3D(K)+AH3D(K) AH3D(K) = 0. ENDIF ENDIF ENDIF !----------------- FOR TENDENCIES DUE TO PRECIPITATION/SEDIMENTATION --- TEMPC = QC3D(K)+QI3D(K) ! QC+QI IF (TEMPC.GT.0.) THEN VTMEAN(K,1) = (VT_QC(K)*QC3D(K)+VT_QI(K)*QI3D(K))/TEMPC ! SEDIMENTATION SPEED OF WET AEROSOL IN QC+QI ELSE VTMEAN(K,1) = 0. ENDIF TEMPR = QR3D(K)+QS3D(K)+QG3D(K)+QH3D(K) ! QR+QS+QG+QH IF (TEMPR.GT.0.) THEN VTMEAN(K,2) = (VT_QR(K)*QR3D(K)+VT_QS(K)*QS3D(K)+VT_QG(K)* & QG3D(K)+VT_QH(K)*QH3D(K))/TEMPR ! SEDIMENTATION SPEED OF WET AEROSOL IN QR+QS+QG+QH ELSE VTMEAN(K,2) = 0. ENDIF TEMPT = TEMPC+TEMPR ! QC+QR+QI+QS+QG IF (TEMPT.GT.0.) THEN VTMEAN(K,3) = (VTMEAN(K,1)*TEMPC+VTMEAN(K,2)*TEMPR)/TEMPT ! MEAN SEDIMENTATION SPEED OF WET AEROSOL ! IN QC+QR+QI+QS+QG+QH ELSE VTMEAN(K,3) = 0. ENDIF ENDDO !----------------- CALCULATED FALLING TERMS ---------------------------- CALL PTFLUX(QC3D(1),VT_QC(1),RHO(1),DZ3D(1),KTE,DT,DTMN,CLODNCV) CALL PTFLUX(QR3D(1),VT_QR(1),RHO(1),DZ3D(1),KTE,DT,DTMN,RAINNCV) CALL PTFLUX(QI3D(1),VT_QI(1),RHO(1),DZ3D(1),KTE,DT,DTMN,ICENCV) CALL PTFLUX(QS3D(1),VT_QS(1),RHO(1),DZ3D(1),KTE,DT,DTMN,SNOWNCV) CALL PTFLUX(QG3D(1),VT_QG(1),RHO(1),DZ3D(1),KTE,DT,DTMN,GRAPNCV) CALL PTFLUX(QH3D(1),VT_QH(1),RHO(1),DZ3D(1),KTE,DT,DTMN,HAILNCV) RAINNC = RAINNC+CLODNCV+RAINNCV+ICENCV+SNOWNCV+GRAPNCV+HAILNCV SNOWNC = SNOWNC+ICENCV+SNOWNCV GRAPNC = GRAPNC+GRAPNCV HAILNC = HAILNC+HAILNCV SR = (ICENCV+SNOWNCV+GRAPNCV+HAILNCV)/ & (CLODNCV+RAINNCV+ICENCV+SNOWNCV+GRAPNCV+HAILNCV+1.E-12) CALL FLFLUX(NC3D(1),VT_NC(1),RHO(1),DZ3D(1),KTE,DT) CALL FLFLUX(NR3D(1),VT_NR(1),RHO(1),DZ3D(1),KTE,DT) CALL FLFLUX(NI3D(1),VT_NI(1),RHO(1),DZ3D(1),KTE,DT) CALL FLFLUX(NS3D(1),VT_NS(1),RHO(1),DZ3D(1),KTE,DT) CALL FLFLUX(NG3D(1),VT_NG(1),RHO(1),DZ3D(1),KTE,DT) CALL FLFLUX(NH3D(1),VT_NH(1),RHO(1),DZ3D(1),KTE,DT) CALL PTFLUX(VI3D(1),VT_VI(1),RHO(1),DZ3D(1),KTE,DT,DTMN,VINCV) CALL PTFLUX(VS3D(1),VT_VS(1),RHO(1),DZ3D(1),KTE,DT,DTMN,VSNCV) CALL PTFLUX(VG3D(1),VT_VG(1),RHO(1),DZ3D(1),KTE,DT,DTMN,VGNCV) CALL PTFLUX(FI3D(1),VT_FI(1),RHO(1),DZ3D(1),KTE,DT,DTMN,FINCV) CALL PTFLUX(FS3D(1),VT_FS(1),RHO(1),DZ3D(1),KTE,DT,DTMN,FSNCV) CALL PTFLUX(AI3D(1),VT_AI(1),RHO(1),DZ3D(1),KTE,DT,DTMN,AINCV) CALL PTFLUX(AS3D(1),VT_AS(1),RHO(1),DZ3D(1),KTE,DT,DTMN,ASNCV) CALL PTFLUX(AG3D(1),VT_AG(1),RHO(1),DZ3D(1),KTE,DT,DTMN,AGNCV) CALL PTFLUX(AH3D(1),VT_AH(1),RHO(1),DZ3D(1),KTE,DT,DTMN,AHNCV) CALL PTFLUX(I3M3D(1),VTI3M(1),RHO(1),DZ3D(1),KTE,DT,DTMN,I3MNCV) IV0 = 0 DO IV = 1,NAER IF (NAERN(IV).GT.3) THEN CALL FLFLUX(AERO(1,IV0+3),VTMEAN(1,1),RHO(1),DZ3D(1),KTE,DT)! SEDI. OF AEROSOLS IN CLOUDS CALL FLFLUX(AERO(1,IV0+4),VTMEAN(1,2),RHO(1),DZ3D(1),KTE,DT)! SEDI. OF AEROSOLS IN PRECIP. ELSEIF (NAERN(IV).EQ.3) THEN CALL FLFLUX(AERO(1,IV0+3),VTMEAN(1,3),RHO(1),DZ3D(1),KTE,DT)! SEDI. OF AEROSOLS IN HYDROM. ENDIF IV0 = IV0+NAERN(IV) ENDDO !---------------------------------------------------------------------- DO K = KTS,KTE TK3D(K) = TAIR(K) QV3D(K) = QVAP(K) IF (QC3D(K).GE.QSMALL.AND.NC3D(K).GE.NSMALL) THEN CALL SOLVE_AFAC(TK3D(K),QC3D(K),NC3D(K),LAMC(K),MVDC(K), & AFAC(K)) ELSE QC3D(K) = 0.; NC3D(K) = 0.; MVDC(K) = 0.; AFAC(K) = 0. LAMC(K) = 0. ENDIF IF (QR3D(K).GE.QSMALL.AND.NR3D(K).GE.NSMALL) THEN CALL SOLVE_AFAR(TK3D(K),QR3D(K),NR3D(K),LAMR(K),MVDR(K), & AFAR(K)) ELSE QR3D(K) = 0.; NR3D(K) = 0.; MVDR(K) = 0.; LAMR(K) = 0. AFAR(K) = 0. ENDIF IF (QI3D(K).GE.QSMALL.AND.NI3D(K).GE.NSMALL) THEN CALL SOLVE_AFAI(TK3D(K),PRES(K),RHO(K),QV3D(K),QI3D(K), & NI3D(K),VI3D(K),FI3D(K),AI3D(K),I3M3D(K),ADAGR(K), & ZETA(K),LAMI(K),AFAI(K),MVDI(K),RHOI(K),AMI(K),BMI(K),& AVI(K),BVI(K),BEST(K)) IASPR(K) = FI3D(K)/I3M3D(K) ELSE QI3D(K) = 0.; NI3D(K) = 0.; MVDI(K) = 0.; I3M3D(K) = 0. FI3D(K) = 0.; VI3D(K) = 0.; AFAI(K) = 0.; IASPR(K) = 1. ADAGR(K) = 1.; RHOI(K) = 0.; AI3D(K) = 0.; LAMI(K) = 0. ENDIF IF (AS3D(K).LT.ASMALL) THEN AS3D(K) = 0. ENDIF IF (QS3D(K).GE.QSMALL.AND.NS3D(K).GE.NSMALL) THEN CALL SOLVE_AFAS(TK3D(K),RHO(K),QS3D(K),QC3D(K),NS3D(K), & VS3D(K),FS3D(K),AS3D(K),AFAS(K),LAMS(K),MVDS(K), & RHOS(K),SASPR(K),AMS(K),AVS(K),BVS(K)) ELSE QS3D(K) = 0.; NS3D(K) = 0.; VS3D(K) = 0.; AS3D(K) = 0. MVDS(K) = 0.; RHOS(K) = 0.; AFAS(K) = 0.; LAMS(K) = 0. FS3D(K) = 0.; SASPR(K) = 1. ENDIF IF (AG3D(K).LT.ASMALL) THEN AG3D(K) = 0. ENDIF IF (QG3D(K).GE.QSMALL.AND.NG3D(K).GE.NSMALL) THEN CALL SOLVE_AFAG(TK3D(K),RHO(K),QG3D(K),QC3D(K),NG3D(K), & VG3D(K),AG3D(K),LAMG(K),AFAG(K),MVDG(K),RHOG(K), & AMG(K),AVG(K),BVG(K)) ELSE QG3D(K) = 0.; NG3D(K) = 0.; MVDG(K) = 0.; VG3D(K) = 0. AG3D(K) = 0.; AFAG(K) = 0.; RHOG(K) = 0.; LAMG(K) = 0. ENDIF IF (AH3D(K).LT.ASMALL) THEN AH3D(K) = 0. ENDIF IF (QH3D(K).GE.QSMALL.AND.NH3D(K).GE.NSMALL) THEN CALL SOLVE_AFAH(TK3D(K),RHO(K),QH3D(K),NH3D(K),AH3D(K), & LAMH(K),AFAH(K),MVDH(K),AVH(K),BVH(K)) ELSE QH3D(K) = 0.; NH3D(K) = 0.; MVDH(K) = 0.; LAMH(K) = 0. AH3D(K) = 0.; AFAH(K) = 0. ENDIF IF ((QC3D(K)+QR3D(K)).LT.QSMALL) THEN AERO(K,1) = MAX(RLIMIT,AERO(K,1)+AERO(K,3)) AERO(K,1) = MIN(AERO(K,1),AERO(K,2)) AERO(K,3) = 0. ENDIF IF ((QR3D(K)+QS3D(K)+QG3D(K)+QH3D(K)).LT.QSMALL) THEN AERO(K,1) = MAX(RLIMIT,AERO(K,1)+AERO(K,4)) AERO(K,1) = MIN(AERO(K,1),AERO(K,2)) AERO(K,4) = 0. ENDIF QDCN3D(K) = AERO(K,1) QTCN3D(K) = AERO(K,2) QCCN3D(K) = AERO(K,3) QRCN3D(K) = AERO(K,4) QNIN3D(K) = AERO(K,5) !------------------------------------------------------------------------ ENDDO END SUBROUTINE NTU_MICRO !====================================================================== !====================================================================== SUBROUTINE ACTIVA(TK1D,W1D,NC1D,NR1D,ABCD,QACac,QACar,RX0,ZCCNS, & SSW,IAE,DO_FR) !====================================================================== IMPLICIT NONE INTEGER :: IAE ! (I) aerosol component REAL :: TK1D,NC1D,NR1D,QACac,QACar,NACac,NACar,RX0,SSW REAL, INTENT(IN) :: W1D DOUBLE PRECISION :: X1,X3,DMASS,DERF REAL, DIMENSION(NCCN) :: ZCCNS REAL, DIMENSION(MAER) :: ABCD LOGICAL :: DO_FR INTEGER :: IM REAL :: GAMAA,ALPHA,SSRW,ES0,RX9,RX1,XMASS,RACT,ZCN9,ZCN99, & ZCCN10,ZCCN15,ZCCN20,ZCCN25,RACT10,RACT15,RACT20,RACT25, & WR10,WR15,WR20,WR25,WR10TO15,WR15TO20,WR20TO25,WR25PLUS REAL, PARAMETER :: RS10 = 10.E-6, RS15 = 15.E-6, RS20 = 20.E-6, & RS25 = 25.E-6 QACac = 0.; QACar = 0.; NACac = 0.; NACar = 0. !------ decide the number of each mode and radius for CCN cut-off ------ IF (DO_FR) THEN CALL FIND_RC0(DBLE(ABCD(1)/ABCD(2)),CNMOD(1,IAE),CNSTD(1,IAE),& WMAS(1,IAE),RX0,TBLXF(1,IAE),TBLRC) DO_FR = .FALSE. ENDIF IF (RX0.EQ.1.E-9) RETURN ! RC=RC_MIN, no dry aerosol available RX1 = RX0 !----------------- Kohler Curve parameters ----------------------------- SSRW = MIN(SSW,0.03) ES0 = 0.0761-1.55E-4*(TK1D-TK0C) ! [J/M**2] ALPHA = 2.*ES0/(RV*TK1D*RHOW) ! [J/m^2]/[J/kg/K]/[K]/[kg/m^3] = [m] GAMAA = 4.*ALPHA**3./(27.*BETA1(IAE)) ! [m^3] RX9 = (GAMAA/(SSRW**2.))**THRD ! critical dry radius [m] RX9 = MAX(RX9,RXMIN(IAE)) ! set lower limit of RX IF (RX1.LE.RX9) RETURN !----------------- RAINDROP ACTIVATION --------------------------------- IF (RX1.GT.RS10) THEN ! force all aerosol > RS10 into rain DMASS = 0.D+0 DO IM = 1,NCCN X3 = DLNX(RX1,CNMOD(IM,IAE),CNSTD(IM,IAE),3) X1 = DLNX(RS10,CNMOD(IM,IAE),CNSTD(IM,IAE),3) DMASS = DMASS+DBLE(ZCCNS(IM))/DBLE(RFACT(IM,IAE))*(DERF(X3)-DERF(X1))/2.D+0 ! aerosol into rain ENDDO XMASS = MAX(0.,REAL(DMASS)) XMASS = MIN(ABCD(1),XMASS) IF (NAERN(IAE).GE.4) THEN ABCD(4) = ABCD(4)+XMASS ! update aerosol in rain ELSEIF (NAERN(IAE).EQ.3) THEN ABCD(3) = ABCD(3)+XMASS ! update aerosol in (cloud+rain) ENDIF ABCD(1) = ABCD(1)-XMASS ! update dry aerosol CALL RSWHITBY(W1D,RS10,ZCCN10,RACT10,BETA1(IAE),ALPHA,ZCCNS, & CNMOD(1,IAE),CNSTD(1,IAE)) CALL RSWHITBY(W1D,RS15,ZCCN15,RACT15,BETA1(IAE),ALPHA,ZCCNS, & CNMOD(1,IAE),CNSTD(1,IAE)) CALL RSWHITBY(W1D,RS20,ZCCN20,RACT20,BETA1(IAE),ALPHA,ZCCNS, & CNMOD(1,IAE),CNSTD(1,IAE)) CALL RSWHITBY(W1D,RS25,ZCCN25,RACT25,BETA1(IAE),ALPHA,ZCCNS, & CNMOD(1,IAE),CNSTD(1,IAE)) WR10 = C4PI3W*RACT10**3. ! [kg] WR15 = C4PI3W*RACT15**3. ! [kg] WR20 = C4PI3W*RACT20**3. ! [kg] WR25 = C4PI3W*RACT25**3. ! [kg] WR10TO15 = 0.5*(WR10+WR15)*(ZCCN10-ZCCN15) ! [kg/kg] WR15TO20 = 0.5*(WR15+WR20)*(ZCCN15-ZCCN20) ! [kg/kg] WR20TO25 = 0.5*(WR20+WR25)*(ZCCN20-ZCCN25) ! [kg/kg] WR25PLUS = WR25*ZCCN25 ! [kg/kg] NACar = ZCCN10 ! [1/kg] QACar = WR10TO15+WR15TO20+WR20TO25+WR25PLUS ! [kg/kg] RX1 = RS10 ! [m] !----------------- end sectionalize rain embryo ------------------------ !----------------- single category rain embryo ------------------------- NR1D = NR1D+NACar ! [ #/kg] ENDIF !----------------- CLOUD DROP ACTIVATION ------------------------------ IF (RX1.GT.RX9) THEN DMASS = 0.D+0 DO IM = 1,NCCN X3 = DLNX(RX1,CNMOD(IM,IAE),CNSTD(IM,IAE),3) X1 = DLNX(RX9,CNMOD(IM,IAE),CNSTD(IM,IAE),3) DMASS = DMASS+DBLE(ZCCNS(IM))/DBLE(RFACT(IM,IAE))*(DERF(X3)-DERF(X1))/2.D+0 ! aerosol into rain ENDDO XMASS = MAX(0.,REAL(DMASS)) XMASS = MIN(XMASS,ABCD(1)) IF (NAERN(IAE).GE.3) THEN ABCD(3) = ABCD(3)+XMASS ! update wet aerosol in (cloud+rain) ENDIF ABCD(1) = ABCD(1)-XMASS ! update dry aerosol CALL CCNWHITBY(W1D,RS10,RX9,ZCN99,RACT,BETA1(IAE),ALPHA,ZCCNS,& CNMOD(1,IAE),CNSTD(1,IAE)) CALL CCNWHITBY(W1D,RX1,RX9,ZCN9,RACT,BETA1(IAE),ALPHA,ZCCNS, & CNMOD(1,IAE),CNSTD(1,IAE)) NACac = MIN(ZCN9,MAX(ZCN99-NC1D,0.)) ! Total # should less than ZCN99 QACac = NACac*C4PI3W*RACT**3. ! [kg/kg] NC1D = NC1D+NACac ! [#/kg] RX1 = RX9 ENDIF RX0 = RX1 END SUBROUTINE ACTIVA !====================================================================== ! !----------------- FUNCTION FOR CALCULATE lnX -------------------------- FUNCTION DLNX(RX,XMODE,SIGMA,N) !====================================================================== INTEGER :: N REAL :: RX,XMODE,SIGMA DOUBLE PRECISION :: DRX,DXMODE,DSTDV,DLNX,DLOG,DSQRT DRX = DBLE(RX) DXMODE = DBLE(XMODE) DSTDV = DBLE(SIGMA) DLNX = (DLOG(DRX/DXMODE)-DSTDV**2.*DBLE(N))/(DSQRT(2.D+0)*DSTDV) END FUNCTION DLNX !----------------- FUNCTION FOR CALCULATE lnX -------------------------- FUNCTION DLNX2(DRX,XMODE,SIGMA,N) INTEGER :: N REAL :: XMODE,SIGMA DOUBLE PRECISION :: DRX,DXMODE,DSTDV,DLNX2,DLOG,DSQRT DXMODE = DBLE(XMODE) DSTDV = DBLE(SIGMA) DLNX2 = (DLOG(DRX/DXMODE)-DSTDV**2.*DBLE(N))/(DSQRT(2.D+0)*DSTDV) END FUNCTION DLNX2 !====================================================================== ! !====================================================================== SUBROUTINE RSWHITBY(W1D,RSX,Z,RACT,BETA1,ALPHA,ZCCN,CNMOD,CNSTD) !====================================================================== IMPLICIT NONE !----------------- calculate number of CCN with dry radii > R0 --------- INTEGER :: IM REAL :: RSX, &! smallest dry radius of rain embryo (input) RWS0, &! wet radius at S = 1 RACT, &! activation radius (output) Z, &! number of CCN that fits the criterion (output) BETA1,ALPHA,TEMP,TEMP1,TEMP2,TMPX,TEMP4,ERF REAL, INTENT(IN) :: W1D REAL, DIMENSION(NCCN) :: ZCCN,CNMOD,CNSTD !----- integrate the number concentration of aerosols with size > RACT-- Z = 0. DO IM = 1,NCCN TEMP1 = ZCCN(IM)/2. TEMP2 = SQRT2*CNSTD(IM) TEMP = (LOG(RSX/CNMOD(IM)))/TEMP2 Z = Z+TEMP1*(1.-ERF(TEMP)) END DO RWS0 = SQRT(BETA1*RSX**3./ALPHA) ! wet radius at S = 1 TMPX = LOG10(RSX*1.E6) TEMP4 = -0.61425115+(-0.66624878-0.17367658*TMPX)*TMPX RACT = 10.**TEMP4*RWS0*W1D**(-0.11782) END SUBROUTINE RSWHITBY !---------------------------------------------------------------------- SUBROUTINE CCNWHITBY(W1D,RSX,RS9,ZCN9,RACT,BETA1,ALPHA,ZCCN, & CNMOD,CNSTD) !====================================================================== IMPLICIT NONE !---- derive number of cloud drop activation according to Whitby's distributions INTEGER :: IM REAL :: RSX, &! previously activated smallest dry radius RS9, &! current activated smallest dry radius ZCN9, &! newly activated number concentration (#/mol) RACT, &! wet radius of the activated drop RWS0, &! wet radius at S = 1 BETA1,ALPHA,TMPX,TMP1,TMP2,TMP3,TMP4,TMP5,TMP6,ERF REAL, INTENT(IN) :: W1D REAL, DIMENSION(NCCN) :: ZCCN,CNMOD,CNSTD RWS0 = SQRT(BETA1*RS9**3./ALPHA) ! wet radius at S = 1 TMPX = LOG10(RSX*1.E6) TMP1 = -0.61425115+(-0.66624878-0.17367658*TMPX)*TMPX RACT = 10.**TMP1*RWS0*W1D**(-0.11782) ! activation size ZCN9 = 0. DO IM = 1,NCCN TMP2 = ZCCN(IM)/2. TMP3 = SQRT2*CNSTD(IM) TMP4 = (LOG(RSX/CNMOD(IM)))/TMP3 TMP5 = (LOG(RS9/CNMOD(IM)))/TMP3 TMP6 = TMP2*(ERF(TMP4)-ERF(TMP5)) ZCN9 = ZCN9+TMP6 END DO END SUBROUTINE CCNWHITBY !====================================================================== ! !====================================================================== SUBROUTINE DEACTIVA(NACcv,QC1D,NC1D,DCN,TCN,WCN,RC,ZCCNS,IAE, & DO_FR) !====================================================================== IMPLICIT NONE INTEGER :: IAE ! (I) aerosol component LOGICAL :: DO_FR REAL :: NC1D,NACcv,QC1D,DCN,TCN,WCN, &! (I/O) wet aerosol mixing ratio [kg/kg] RC ! (I/O) R-CUTOFF [um] REAL, DIMENSION(NCCN) :: ZCCNS ! (I) # of each mode in [#/kg] INTEGER :: I,IM DOUBLE PRECISION :: DNACcv, &! # of deactivated cloud droplet DMASS, &! mass of deactivated wet aerosol DSQRT2,X0,X1,X2,X3,S1,S2,DERF,DEXP,DSQRT !----------------- deactivation of cloud drops ------------------------- DSQRT2 = DSQRT(2.D+0) DNACcv = DBLE(MIN(NACcv,NC1D)) IF (QC1D.LT.RSMALL) THEN DNACcv = DBLE(NC1D) ! deactivate all cld drop for negligiable QC NC1D = RLIMIT ENDIF DMASS = 0.D+0 IF (DNACcv.GT.1.D+1) THEN !------ decide the number of each mode and radius for CCN cut-off ------ IF (DO_FR) THEN CALL FIND_RC0(DBLE(DCN/TCN),CNMOD(1,IAE),CNSTD(1,IAE), & WMAS(1,IAE),RC,TBLXF(1,IAE),TBLRC) DO_FR = .FALSE. ENDIF S1 = 0.D+0 DO IM = 1,NCCN X0 = DLNX(RC,CNMOD(IM,IAE),CNSTD(IM,IAE),0) S1 = S1+DBLE(ZCCNS(IM))*(1.D+0-DERF(X0))*5.D-1 ! # for all with r > RC ENDDO IF (DNACcv.GT.S1) THEN DCN = DCN+WCN WCN = 0. CALL FIND_RC0(DBLE(DCN/TCN),CNMOD(1,IAE),CNSTD(1,IAE), & WMAS(1,IAE),RC,TBLXF(1,IAE),TBLRC) RETURN ENDIF S1 = 0. DO IM = 1,NCCN X0 = DLNX(RC,CNMOD(IM,IAE),CNSTD(IM,IAE),0) S1 = S1+DBLE(ZCCNS(IM)/CNSTD(IM,IAE))/DSQRT2*DPDF(X0) ! # between 0 & RC ENDDO X1 = DBLE(RC) X2 = DBLE(RC)*DEXP(DNACcv/S1) DO I = 1,4 ! iteratively fixed dN 3 times S1 = 0.D+0 S2 = 0.D+0 DO IM = 1,NCCN X0 = DLNX2(X1,CNMOD(IM,IAE),CNSTD(IM,IAE),0) X3 = DLNX2(X2,CNMOD(IM,IAE),CNSTD(IM,IAE),0) S1 = S1+DBLE(ZCCNS(IM)/CNSTD(IM,IAE))/DSQRT2*DPDF(X3) S2 = S2+DBLE(ZCCNS(IM))*(DERF(X3)-DERF(X0))*5.D-1 ! # between X1 & X2 ENDDO DNACcv = DNACcv-S2 X1 = X2 X2 = X2*DEXP(DNACcv/S1) ENDDO X1 = DBLE(RC) X2 = MAX(X1,X2) DO IM = 1,NCCN X0 = DLNX2(X1,CNMOD(IM,IAE),CNSTD(IM,IAE),3) X3 = DLNX2(X2,CNMOD(IM,IAE),CNSTD(IM,IAE),3) DMASS = DMASS+DBLE(ZCCNS(IM))/DBLE(RFACT(IM,IAE))*(DERF(X3)-DERF(X0))*5.D-1 ! MASS between X0 & X3 ENDDO DCN = DCN+MIN(WCN,REAL(DMASS)) WCN = WCN-MIN(WCN,REAL(DMASS)) RC = REAL(X2) ENDIF END SUBROUTINE DEACTIVA !====================================================================== ! !====================================================================== SUBROUTINE ICENU(TK1D,P1D,DT,RHO,QV1D,QI1D,NI1D,VI1D,FI1D,AI1D, & I3M1D,QNIN) !====================================================================== IMPLICIT NONE INTEGER :: IDEPNU REAL :: TK1D,P1D,DT,RHO,QV1D,QI1D,NI1D,VI1D,FI1D,AI1D,I3M1D,QNIN,& QVSI,XXLS,ESI,QNDvi,NNDvi,VNDvi,FNDvi,ANDvi,INDvi,NVI0, & EPA,RGDEP,SFCTNV,ICED,GGDEP,SRI,ARDEP0,ARDEP,COSM1,GEOF1,& IJDEP0,IJDEP,DC1,DC2,DC3,TC1D,INR0,DANGLE,DACTE,QVDMAX, & NNUMAX,SSRI,CPM QNDvi = 0.;NNDvi = 0.;VNDvi = 0.;FNDvi = 0.;ANDvi = 0.;INDvi = 0. IDEPNU = 8 ! Ice deposition-nucleation equation option TC1D = TK1D-TK0C XXLS = 3.15E6-2.37E3*TK1D+0.3337E6 ESI = MIN(0.99*P1D,POLYSVP(TK1D,1)) QVSI = 0.622*ESI/(P1D-ESI) SRI = QV1D/QVSI SSRI = QV1D/QVSI-1. QVDMAX = (QV1D-QVSI)/(1.+XXLS**2.*QV1D/(CP*RV*TK1D**2.)) ! SUPERSATURATED WATER FIXED DUE TO TEMPERATURE CHANGE NNUMAX = 0.98*QVDMAX/MI0 ! INITIATION OF CLOUD ICE CONSTRAINED by 0.98*SUPICE ICED = 916.7-0.175*TC1D-5.E-4*TC1D**2. ! ICE DENSITY SRI = QV1D/QVSI ! SATURATION RATIO OVER ICE SFCTNV = ((76.1-0.155*TC1D)+(28.5+0.25*TC1D))*1.E-3 ! SURFACE TENSION OVER ICE/VAPOR IF (IDEPNU.EQ.0) THEN ! NO DEPOSITION NUCLEATION NNDvi = 0. ELSEIF (IDEPNU.EQ.1) THEN ! Fletcher (1962) NVI0 = 1.E-2*EXP(0.6*(TK0C-MAX(TK1D,2.46E2))) NNDvi = MIN(QNIN,NNUMAX,MAX(0.,NVI0/RHO-NI1D)) ELSEIF (IDEPNU.EQ.2) THEN ! Cooper NVI0 = 5.*EXP(0.304*(TK0C-MAX(TK1D,2.46E2))) NNDvi = MIN(QNIN,NNUMAX,MAX(0.,NVI0/RHO-NI1D)) ELSEIF (IDEPNU.EQ.3) THEN ! Huffman (1973) NVI0 = 1007.08*SSRI**4.5*1.E3 NNDvi = MIN(QNIN,NNUMAX,MAX(0.,NVI0/RHO-NI1D)) ELSEIF (IDEPNU.EQ.4) THEN ! Meyers et al (1992) NVI0 = 1.E3*EXP(1.296E1*SSRI-6.39E-1) NNDvi = MIN(QNIN,NNUMAX,MAX(0.,NVI0/RHO-NI1D)) ELSEIF (IDEPNU.EQ.5.AND.INSPEC.EQ.4) THEN ! Bacteria(P. syringae):Yankofsky 1981 NVI0 = QNIN*1.1776*EXP(-89318./((LOG(SRI))**2.*TK1D**3.)) NNDvi = MIN(QNIN,NNUMAX,MAX(0.,NVI0*DT)) ELSEIF (IDEPNU.EQ.6) THEN ! DUST (Field 2006)(Asia) NVI0 = QNIN*4.337E-1*EXP(-2.E6/((LOG(SRI))**2.*TK1D**3.)) NNDvi = MIN(QNIN,NNUMAX,MAX(0.,NVI0*DT)) ELSEIF (IDEPNU.EQ.7) THEN IF (INSPEC.EQ.1) THEN ! SOOT DC1 = -0.5411955; DC2 = 1.879918; DC3 = 1.607947 ! CHEN Table3 INR0 = 0.4E-7; DANGLE = 28.; DACTE = -20.E-20 ! HOOSE Table3 ELSEIF (INSPEC.EQ.2) THEN ! SAHARAN DUST DC1 = -0.3353619; DC2 = 1.990979; DC3 = 2.175539 ! CHEN Table3 INR0 = 1.75E-7; DANGLE = 5.06; DACTE = 3.35E-20 ! CHEN Table6 ELSEIF (INSPEC.EQ.3) THEN ! ASIAN DUST DC1 = -0.3598818; DC2 = 1.982032; DC3 = 2.025390 ! CHEN Table3 INR0 = 2.E-7; DANGLE = 8.1; DACTE = 1.82E-20 ! CHEN Table6 ENDIF EPA = ESI*SRI RGDEP = 2.*2.99E-26*SFCTNV/(ICED*BOLTZ*TK1D*LOG(SRI)) ! GERM RADIUS FOR DEPOSITION NUCLEATION GGDEP = C4PI3*SFCTNV*RGDEP**2. ! HOMOGENEOUS ENERGY ARDEP0 = SFCTNV**5.E-1*(BOLTZ*TK1D)**(-1.5)/1.E13 ARDEP = EPA**2.*ARDEP0/ICED COSM1 = LOG(1.-COS(DANGLE*PI/1.8E2)) GEOF1 = MIN(1.,EXP(DC1+DC2*COSM1+DC3*RGDEP/INR0)) IJDEP0 = EXP((-1.*DACTE-GEOF1*GGDEP)/(BOLTZ*TK1D)) IJDEP = ARDEP*INR0**2.*GEOF1**5.E-1*IJDEP0 NVI0 = QNIN*(1.-EXP(-IJDEP*DT)) NNDvi = MIN(QNIN,NNUMAX,MAX(0.,NVI0)) ELSEIF (IDEPNU.EQ.8) THEN NVI0 = MAX(0.,5.94E-5*((TK0C-TK1D)**3.33)*((1.E6/1.E6)** &! set naer05=1e6/kg (derived to cm-3 at STP below) (0.0264*(TK0C-TK1D)+3.3E-3))*1.E3) ! units of #/kg NNDvi = MIN(QNIN,NNUMAX,MAX(0.,NVI0/RHO-NI1D)) ENDIF QNDvi = NNDvi*MI0 VNDvi = QNDvi*iRHOI0 FNDvi = QNDvi*1.*iAMI0 ! ISOMETRIC ANDvi = (KCIMIN*NNDvi*INDvi**2.)**THRD INDvi = QNDvi*iAMI0 QNIN = MAX(0.,QNIN-NNDvi) QV1D = MAX(0.,QV1D-QNDvi) QI1D = MAX(0.,QI1D+QNDvi) NI1D = MAX(0.,NI1D+NNDvi) IF (ICE_RHOI.EQ.0.OR.ICE_RHOI.EQ.2) THEN VI1D = 0. ELSEIF (ICE_RHOI.EQ.1) THEN VI1D = MAX(0.,VI1D+VNDvi) ENDIF FI1D = MAX(0.,FI1D+FNDvi) I3M1D = MAX(0.,I3M1D+INDvi) IF (AFAI_3M.EQ.0.OR.AFAI_3M.EQ.2) THEN AI1D = 0. ELSEIF (AFAI_3M.EQ.1) THEN AI1D = MAX(0.,AI1D+ANDvi) ENDIF CPM = CP*(1.+0.887*QV1D) XXLS = 3.15E6-2370.*TK1D+0.3337E6 TK1D = TK1D+XXLS*QNDvi/CPM END SUBROUTINE ICENU !====================================================================== ! !====================================================================== SUBROUTINE SEDI_FALL(TK1D,P1D,QV1D,QC1D,QR1D,QI1D,QS1D,QG1D,QH1D,& NC1D,NR1D,NI1D,NS1D,NG1D,NH1D,VI1D,VS1D,VG1D,FI1D, & FS1D,AI1D,AS1D,AG1D,AH1D,I3M1D,VTQC,VTQR,VTQI,VTQS, & VTQG,VTQH,VTNC,VTNR,VTNI,VTNS,VTNG,VTNH,VTVI,VTVS, & VTVG,VTFI,VTFS,VTAI,VTAS,VTAG,VTAH,VTI3M,SASPR,RHO) !====================================================================== IMPLICIT NONE REAL :: TK1D,P1D,QV1D,QC1D,QR1D,QI1D,QS1D,QG1D,QH1D,NC1D,NR1D, & NI1D,NS1D,NG1D,NH1D,VI1D,VS1D,VG1D,FI1D,FS1D,AI1D,AS1D, & AG1D,AH1D,I3M1D,VTQC,VTQR,VTQI,VTQS,VTQG,VTQH,VTNC,VTNR, & VTNI,VTNS,VTNG,VTNH,VTVI,VTVS,VTVG,VTFI,VTFS,VTAI,VTAS, & VTAG,VTAH,VTI3M,QRHO,RHOAJ,RHO,MVRC,MVRR,GUC,GUR,LAMC, & LAMR,LAMI,LAMS,LAMG,LAMH,RHOI,RHOS,RHOG,AFAC,AFAR,AFAI, & AFAS,AFAG,AFAH,ADAGR,ZETA3,AMI,BMI,AMS,AMG,AVI,BVI,AVS, & BVS,AVG,BVG,AVH,BVH,ZETA,FSQC,FSQR,FSQI,FSQS,FSQG,FSQH, & FSNC,FSNR,FSNI,FSNS,FSNG,FSNH,FSVI,FSVS,FSVG,FSFI,FSAI, & FSAS,FSAG,FSAH,MVDC,MVDR,MVDI,MVDS,MVDG,MVDH,BEST,SASPR REAL, PARAMETER :: AVTC = 8.8462E+02, BVTC = 9.7593E+07 REAL, PARAMETER :: CVTC = -3.4249E+11, AVTR = 2.1454E+00 REAL, PARAMETER :: BVTR = -2.2812E-04, CVTR = 2.9676E-09 REAL, PARAMETER :: CQC1 = 2.0901E+01, CQC2 = 9.9111E-01 REAL, PARAMETER :: CQC3 = 4.4182E+00, CNC1 = 1.8276E+01 REAL, PARAMETER :: CNC2 = 1.0015E+00, CNC3 = 1.9838E+00 REAL, PARAMETER :: CQR1 = 1.5943E+01, CQR2 = 1.1898E+00 REAL, PARAMETER :: CQR3 = 4.0073E+00, CNR1 = 9.4791E+00 REAL, PARAMETER :: CNR2 = 9.7607E-01, CNR3 = 1.0858E+00 QRHO = SQRT(RHO) RHOAJ = (RHOSU/RHO)**0.54 IF (QC1D.GE.QSMALL) THEN CALL SOLVE_AFAC(TK1D,QC1D,NC1D,LAMC,MVDC,AFAC) IF (LIQ_VTC.EQ.0) THEN FSQC = EXP(GAMLN(BVC0+BMW+AFAC+1.)-GAMLN(BMW+AFAC+1.)- & BVC0*LOG(LAMC)) FSNC = EXP(GAMLN(BVC0+AFAC+1.)-GAMLN(AFAC+1.)-BVC0* & LOG(LAMC)) VTQC = RHOAJ*FSQC*AVC0 VTNC = RHOAJ*FSNC*AVC0 ELSEIF (LIQ_VTC.EQ.1) THEN MVRC = MIN(MAX((QC1D/NC1D/C4PI3W)**THRD,RCMIN),RCMAX) GUC = EXP(EXP(AFU+BFU*(LOG(MVRC))**3.+CFU*QRHO**3.)) VTQC = EXP(CQC1+CQC2*LOG(NC1D)+CQC3*LOG(MVRC))*GUC/QC1D VTNC = EXP(CNC1+CNC2*LOG(NC1D)+CNC3*LOG(MVRC))*GUC/NC1D ELSEIF (LIQ_VTC.EQ.2) THEN MVRC = MIN(MAX((QC1D/NC1D/C4PI3W)**THRD,RCMIN),RCMAX) GUC = EXP(EXP(AFU+BFU*(LOG(MVRC))**3.+CFU*QRHO**3.)) VTQC = MVRC*(AVTC+BVTC*MVRC+CVTC*MVRC**2.)*GUC VTNC = MVRC*(AVTC+BVTC*MVRC+CVTC*MVRC**2.)*GUC ENDIF VTQC = MIN(VTQC,VTCMAX) VTNC = MIN(VTNC,VTCMAX) ENDIF IF (QR1D.GE.QSMALL) THEN CALL SOLVE_AFAR(TK1D,QR1D,NR1D,LAMR,MVDR,AFAR) IF (LIQ_VTR.EQ.0) THEN FSQR = EXP(GAMLN(BVR0+BMW+AFAR+1.)-GAMLN(BMW+AFAR+1.)- & BVR0*LOG(LAMR)) FSNR = EXP(GAMLN(BVR0+AFAR+1.)-GAMLN(AFAR+1.)-BVR0* & LOG(LAMR)) VTQR = RHOAJ*FSQR*AVR0 VTNR = RHOAJ*FSNR*AVR0 ELSEIF (LIQ_VTR.EQ.1) THEN MVRR = MIN(MAX((QR1D/NR1D/C4PI3W)**THRD,RRMIN),RRMAX) GUR = EXP(EXP(AFU+BFU*(LOG(MVRR))**3.+CFU*QRHO**3.)) VTQR = EXP(CQR1+CQR2*LOG(NR1D)+CQR3*LOG(MVRR))*GUR/QR1D VTNR = EXP(CNR1+CNR2*LOG(NR1D)+CNR3*LOG(MVRR))*GUR/NR1D ELSEIF (LIQ_VTR.EQ.2) THEN MVRR = MIN(MAX((QR1D/NR1D/C4PI3W)**THRD,RRMIN),RRMAX) GUR = EXP(EXP(AFU+BFU*(LOG(MVRR))**3.+CFU*QRHO**3.)) VTQR = EXP(AVTR+BVTR/MVRR+CVTR/(MVRR**2.))*GUR VTNR = EXP(AVTR+BVTR/MVRR+CVTR/(MVRR**2.))*GUR ENDIF VTQR = MIN(VTQR,VTRMAX) VTNR = MIN(VTNR,VTRMAX) ENDIF IF (QI1D.GE.QSMALL) THEN CALL SOLVE_AFAI(TK1D,P1D,RHO,QV1D,QI1D,NI1D,VI1D,FI1D,AI1D, & I3M1D,ADAGR,ZETA,LAMI,AFAI,MVDI,RHOI,AMI,BMI,AVI,BVI, & BEST) FSQI = EXP(GAMLN(BVI+BMI+AFAI+1.)-GAMLN(BMI+AFAI+1.)-BVI* & LOG(LAMI)) FSNI = EXP(GAMLN(BVI+AFAI+1.)-GAMLN(AFAI+1.)-BVI*LOG(LAMI)) FSVI = EXP(GAMLN(BVI+AFAI+4.)-GAMLN(AFAI+4.)-BVI*LOG(LAMI)) VTQI = MIN(RHOAJ*FSQI*AVI,VTIMAX) VTNI = MIN(RHOAJ*FSNI*AVI,VTIMAX) VTVI = MIN(RHOAJ*FSVI*AVI,VTIMAX) IF (AI1D.GE.ASMALL) THEN FSAI = EXP(GAMLN(BVI+AFAI+3.)-GAMLN(AFAI+3.)-BVI*LOG(LAMI)) VTAI = MIN(RHOAJ*FSAI*AVI,VTIMAX) ENDIF IF (I3M1D.GE.ISMALL.AND.FI1D.GE.ISMALL) THEN ZETA3 = 3.*(ADAGR-1.)/(ADAGR+2.) FSFI = EXP(GAMLN(BVI+ZETA3+AFAI+4.)-GAMLN(ZETA3+AFAI+4.)- & BVI*LOG(LAMI)) VTFI = MIN(RHOAJ*FSFI*AVI,VTIMAX) VTI3M = MIN(RHOAJ*FSVI*AVI,VTIMAX) ENDIF ENDIF IF (QS1D.GE.QSMALL) THEN CALL SOLVE_AFAS(TK1D,RHO,QS1D,QC1D,NS1D,VS1D,FS1D,AS1D,AFAS, & LAMS,MVDS,RHOS,SASPR,AMS,AVS,BVS) FSQS = EXP(GAMLN(BVS+BMS+AFAS+1.)-GAMLN(BMS+AFAS+1.)-BVS* & LOG(LAMS)) FSNS = EXP(GAMLN(BVS+AFAS+1.)-GAMLN(AFAS+1.)-BVS*LOG(LAMS)) FSVS = EXP(GAMLN(BVS+AFAS+4.)-GAMLN(AFAS+4.)-BVS*LOG(LAMS)) VTQS = MIN(RHOAJ*FSQS*AVS,VTSMAX) VTNS = MIN(RHOAJ*FSNS*AVS,VTSMAX) VTVS = MIN(RHOAJ*FSVS*AVS,VTSMAX) VTFS = VTVS IF (AS1D.GE.ASMALL) THEN FSAS = EXP(GAMLN(BVS+AFAS+3.)-GAMLN(AFAS+3.)-BVS*LOG(LAMS)) VTAS = MIN(RHOAJ*FSAS*AVS,VTSMAX) ENDIF ENDIF IF (QG1D.GE.QSMALL) THEN CALL SOLVE_AFAG(TK1D,RHO,QG1D,QC1D,NG1D,VG1D,AG1D,LAMG,AFAG, & MVDG,RHOG,AMG,AVG,BVG) FSQG = EXP(GAMLN(BVG+BMG+AFAG+1.)-GAMLN(BMG+AFAG+1.)-BVG* & LOG(LAMG)) FSNG = EXP(GAMLN(BVG+AFAG+1.)-GAMLN(AFAG+1.)-BVG*LOG(LAMG)) FSVG = EXP(GAMLN(BVG+AFAG+4.)-GAMLN(AFAG+4.)-BVG*LOG(LAMG)) VTQG = MIN(RHOAJ*FSQG*AVG,VTGMAX) VTNG = MIN(RHOAJ*FSNG*AVG,VTGMAX) VTVG = MIN(RHOAJ*FSVG*AVG,VTGMAX) IF (AG1D.GE.ASMALL) THEN FSAG = EXP(GAMLN(BVG+AFAG+3.)-GAMLN(AFAG+3.)-BVG*LOG(LAMG)) VTAG = MIN(RHOAJ*FSAG*AVG,VTGMAX) ENDIF ENDIF IF (QH1D.GE.QSMALL) THEN CALL SOLVE_AFAH(TK1D,RHO,QH1D,NH1D,AH1D,LAMH,AFAH,MVDH,AVH,BVH) FSQH = EXP(GAMLN(BVH+BMH+AFAH+1.)-GAMLN(BMH+AFAH+1.)-BVH* & LOG(LAMH)) FSNH = EXP(GAMLN(BVH+AFAH+1.)-GAMLN(AFAH+1.)-BVH*LOG(LAMH)) VTQH = MIN(RHOAJ*FSQH*AVH,VTHMAX) VTNH = MIN(RHOAJ*FSNH*AVH,VTHMAX) IF (AH1D.GE.ASMALL) THEN FSAH = EXP(GAMLN(BVH+AFAH+3.)-GAMLN(AFAH+3.)-BVH*LOG(LAMH)) VTAH = MIN(RHOAJ*FSAH*AVH,VTHMAX) ENDIF ENDIF END SUBROUTINE SEDI_FALL !====================================================================== ! !====================================================================== SUBROUTINE PTFLUX(Q1D,VT1D,RHO,DZ,NK,DT,DTMN,PRT1D) !====================================================================== IMPLICIT NONE INTEGER :: K,NK,NS,NSTEP INTEGER, PARAMETER :: MAXSTP = 1000 REAL :: PRT1D, &! input/output ; accumulative prcipitation (cm) DT, &! time step for integration(=2*DTMN*60 for leapfrog scheme) DTMN ! time step (min) REAL, DIMENSION(NK) :: Q1D,VT1D,RHO,DZ,DQDT !------------ FOR TENDENCIES DUE TO PRECIPITATION/SEDIMENTATION ------- NSTEP = 1 DO K = 1,NK NSTEP = MAX(NSTEP,INT(VT1D(K)*DT/DZ(K)+1.)) IF (NSTEP.GT.80) THEN PRINT *,'IN PREFLUX',K,NSTEP,DT,VT1D(K),DZ(K),Q1D(K),RHO(K) ENDIF ENDDO IF (NSTEP.GT.MAXSTP) THEN PRINT *,'NSTEP FOR PRECIP. IS: ',NSTEP,VT1D,DT,DZ STOP ENDIF DO NS = 1,NSTEP DQDT(1) = -VT1D(1)*Q1D(1)/DZ(1) DO K = 2,NK,1 DQDT(K) = (VT1D(K-1)*RHO(K-1)*Q1D(K-1)- & VT1D(K)*RHO(K)*Q1D(K))/(DZ(K)*RHO(K)) ENDDO PRT1D = PRT1D+VT1D(NK)*RHO(NK)*Q1D(NK)*DTMN*60./REAL(NSTEP) ! accumulate precipitation [kg/m^2]; 1 kg/m^2 = 1 mm DO K = 1,NK Q1D(K) = Q1D(K)+DQDT(K)*DT/REAL(NSTEP) ENDDO ENDDO END SUBROUTINE PTFLUX !====================================================================== ! !====================================================================== SUBROUTINE FLFLUX(Q1D,VT1D,RHO,DZ,NK,DT) !====================================================================== IMPLICIT NONE INTEGER :: NSTEP,K,NK,NS INTEGER, PARAMETER :: MAXSTP = 1000 REAL :: DT REAL, DIMENSION(NK) :: Q1D,VT1D,RHO,DZ,DQDT !-- INPUT DATA FOR STEPSIZE RK4 CALCULATION FOR TENDENCIES DUE TO PRECIPITATION/SEDIMENTATION NSTEP = 1 DO K = 1,NK NSTEP = MAX(NSTEP,INT(VT1D(K)*DT/DZ(K)+1.)) ENDDO DO NS = 1,NSTEP DQDT(1) = -VT1D(1)*Q1D(1)/DZ(1) DO K = 2,NK,1 DQDT(K) = (VT1D(K-1)*RHO(K-1)*Q1D(K-1)- & VT1D(K)*RHO(K)*Q1D(K))/(DZ(K)*RHO(K)) ENDDO DO K = 1,NK Q1D(K) = Q1D(K)+DQDT(K)*DT/REAL(NSTEP) ENDDO ENDDO END SUBROUTINE FLFLUX !====================================================================== ! !====================================================================== SUBROUTINE SMALL_DT(DT,DTS,SDTS,DTKDT,DQVDT,DPDT,RHO,TK1D,P1D, & QV1D,QC1D,QR1D,QI1D,QS1D,QG1D,QH1D,NC1D,NR1D,NI1D, & NS1D,NG1D,NH1D,VI1D,VS1D,VG1D,FI1D,FS1D,AI1D,AS1D, & AG1D,AH1D,I3M1D,SASPR,XDNC,XDNR,GQCTR) !====================================================================== IMPLICIT NONE INTEGER :: I,HID REAL :: DT,DTS,SDTS, &! TIME STEP(s)& SUB-TIME STEP(s)& DTS MUST