Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

CMiPhy.f90 @ 4400

Last change on this file since 4400 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR

Integration of MAR physics

File size: 175.1 KB

Rev	Line
[2089]	1	subroutine CMiPhy
	2
	3	!------------------------------------------------------------------------------+
	4	! Mon 17-Jun-2013 MAR \|
	5	! MAR CMiPhy \|
	6	! subroutine CMiPhy contains the MAR Cloud Microphysical Scheme \|
	7	! \|
	8	! version 3.p.4.1 created by H. Gallee, Thu 21-Mar-2013 \|
	9	! Last Modification by H. Gallee, Mon 17-Jun-2013 \|
	10	! \|
	11	!------------------------------------------------------------------------------+
	12	! \|
	13	! INPUT / OUTPUT: qv__DY(kcolp,mzp) : air specific humidity [kg/kg] \|
	14	! ^^^^^^^^^^^^^^ qw__CM(kcolp,mzp) : cloud drops [kg/kg] \|
	15	! qi__CM(kcolp,mzp) : ice crystals Concentration [kg/kg] \|
	16	! qs__CM(kcolp,mzp) : Snow Particl. Concentration [kg/kg] \|
	17	! qr__CM(kcolp,mzp) : Rain Drops Concentration [kg/kg] \|
	18	! \|
	19	! (to be added: qg__CM(kcolp,mzp) : Graupels Concentration [kg/kg])\|
	20	! \|
	21	! CCNwCM(kcolp,mzp) : cloud droplets number [Nb/m3] \|
	22	! CCNiCM(kcolp,mzp) : ice crystals number [Nb/m3] \|
	23	! \|
	24	! CFraCM(kcolp,mzp) : cloud fraction [-] \|
	25	! \|
	26	! RainCM(kcolp ) : rain Precipitation [m w.e.] \|
	27	! SnowCM(kcolp ) : snow Precipitation [m w.e.] \|
	28	! Ice_CM(kcolp ) : ice Precipitation [m w.e.] \|
	29	! \|
	30	! qid_CM(kcolp,mzp) : Ice Water Formation [kg/kg] \|
	31	! qwd_CM(kcolp,mzp) : Liquid Water Formation [kg/kg] \|
	32	! \|
	33	! INPUT : wa__DY(kcolp,mzp) : Vertical Wind Speed [m/s] \|
	34	! ^^^^^ roa_DY(kcolp,mzp) : Air Density [T/m3] \|
	35	! qvsiCM(kcolp,mzp) : Satur.specific humidity (ICE) [kg/kg] \|
	36	! qvswCM(kcolp,mzp) : Satur.specific humidity (LIQ.) [kg/kg] \|
	37	! \|
	38	! OUTPUT: \|
	39	! ^^^^^^ \|
	40	! \|
	41	! \|
	42	! REFER. : 1) Ntezimana, unpubl.thes.LLN, 115 pp, 1993 \|
	43	! ^^^^^ 2) Lin et al. JCAM 22, 1065--1092, 1983 \|
	44	! (very similar, except that graupels are represented) \|
	45	! 3) Emde and Kahlig, Annal.Geophys. 7, 405-- 414, 1989 \|
	46	! 4) Levkov et al., Contr.Atm.Phys. 65, 35-- 57, 1992 \|
	47	! 5) Meyers et al., JAM 31, 708-- 731, 1992 \|
	48	! (Primary Ice-Nucleation Parameterization) \|
	49	! 6) Delobbe and Gallee, BLM 89, 75-- 107 1998 \|
	50	! (Partial Condensation Scheme) \|
	51	! \|
	52	! # OPTIONS: #hy Additional IF/THEN/ELSE added precluding vectorization \|
	53	! # ^^^^^^^ #qg Graupel Conservation Equation (to verify & include) \|
	54	! # #cn Intercept Parameter / snow gamma distribution = fct(Ta) \|
	55	! # #kk Limitation of SCu fraction \|
	56	! # #VW Cloud Drop. Sediment.(Duynkerke .. 1995, JAS 52, p.2763 \|
	57	! # #pp Emde & Kahlig Ice Crystal Deposition (not included) \|
	58	! # #wi QSat modified by qw/qi Ratio (not included) \|
	59	! \|
	60	! # DEBUG: #WH Additional Output (Each Process is detailled) \|
	61	! # ^^^^^ #WQ FULL Output (Each Process is detailled) \|
	62	! # #wH Additional Output \|
	63	! # #wh Additional Output (Include write CMiPhy_Debug.h) \|
	64	! # #EW Additional Output (Energy and Water Conservation) \|
	65	! # #ew Additional Output (Energy and Water Conservation) \|
	66	! \|
	67	! REMARK : the sign '~' indicates that reference must be verified \|
	68	! ^^^^^^^^ \|
	69	! CAUTION: Partial Condensation Scheme NOT validated \|
	70	! ^^^^^^^ for SCu -- Cu Transition \|
	71	! erf fonction is erroneous on HP \|
	72	! \|
	73	!------------------------------------------------------------------------------+
	74
	75
	76
	77	! Global Variables
	78	! ================
	79
	80	use Mod_Real
	81	use Mod_PHY____dat
	82	use Mod_PHY____grd
	83	use Mod_PHY_CM_dat
	84	use Mod_PHY_CM_grd
	85	use Mod_PHY_CM_kkl
	86	use Mod_PHY_DY_kkl
	87	use Mod_PHY_AT_kkl
	88
	89
	90
	91	! Local Variables
	92	! ================
	93
	94	use Mod_CMiPhy_loc
	95
	96
	97	IMPLICIT NONE
	98
	99
	100	logical :: Heter_Freezng = .FALSE. ! .TRUE. => Levkov et al. (1992) Heterogeneous Freezing of Cloud Droplets
	101	logical :: Homog_Sublima = .FALSE. ! .TRUE. => Levkov et al. (1992) Homogeneous Sublimation of Cloud Ice Particles
	102	logical :: Meyers = .TRUE. ! .TRUE. => Meyers et al. (1992) Ice Nucleation I
	103	logical :: AUTO_w_Sundqv = .TRUE. ! .TRUE. => Sundqvist (1988) Autoconversion Cloud Droplets --> Rain
	104	logical :: AUTO_w_LiouOu = .FALSE. ! .TRUE. => Liou and Ou (1989) Autoconversion Cloud Droplets --> Rain (Tropical SCu ONLY)
	105	logical :: AUTO_w_LinAll = .FALSE. ! .TRUE. => Lin et al. (1983) Autoconversion Cloud Droplets --> Rain
	106	logical :: AUTO_i_Levkov = .TRUE. ! .TRUE. => Levkov et al. (1992) Autoconversion Cloud Ice --> Snow
	107	logical :: AUTO_i_LevkXX = .TRUE. ! .TRUE. => Levkov et al. (1992) Autoconversion Cloud Ice --> Snow (Deposition.Growth)
	108	logical :: AUTO_i_EmdeKa = .FALSE. ! .TRUE. => Emde & Kahlig (1989) Autoconversion Cloud Ice --> Snow
	109	logical :: AUTO_i_Sundqv = .FALSE. ! .TRUE. => Sundqvist Autoconversion Cloud Ice --> Snow
	110	! .FALSE. => Emde & Kahlig (1989) Bergeron-Findeisen Process
	111	logical :: fracSC = .FALSE. ! .TRUE. => Delobbe SCu Fractional Cloudiness may be set up if Frac__Clouds = .TRUE.
	112	logical :: fraCEP = .FALSE. ! .TRUE. => ECMWF SCu Fractional Cloudiness
	113	logical :: HalMos = .TRUE. ! .TRUE. => Levkov et al. (1992) Ice Nucleation II (Hallet-Mossop Process)
	114	logical :: graupel_shape = .TRUE. ! .TRUE. => Snow Particles Shape: Graupellike Snow Flakes of Hexagonal Type
	115	logical :: planes__shape = .FALSE. ! .TRUE. => Snow Particles Shape: Unrimed Side Planes
	116	logical :: aggrega_shape = .FALSE. ! .TRUE. => Snow Particles Shape: Aggregates of unrimed radiating assemblages
	117
	118	logical :: NO_Vec = .TRUE. ! .TRUE. => Preference of IF/THEN/ELSE to sign and max/min Functions
	119
	120	real(kind=real8) :: rad_ww ! Cloud Droplets Radius [...]
	121	real(kind=real8) :: qvs_wi ! Saturation Specific Humididy over Liquid Water [kg/kg]
	122	! Ref.: Emde & Kahlig 1989, Ann.Geophys. 7, p.407 (5)
	123	real(kind=real8) :: BNUCVI ! Nucleation I: Deposition & Condensation-Freez. Nucl. [kg/kg]
	124	real(kind=real8) :: SSat_I ! Nucleation I: Sursaturation % ICE [%]
	125
	126	real(kind=real8) :: Flag_T_NuId ! Flag: T < T_NuId => Nucleation I (Dep./Cond.) may occur [-]
	127	real(kind=real8) :: CCNiId ! [Nb/m3]
	128
	129	real(kind=real8) :: Flag___NuIc ! Flag: 1 => Nucleation I (Cont.Frz.) may occur [-]
	130	real(kind=real8) :: Flag_T_NuIc ! Flag: T < T_NuIc => Nucleation I (Cont.Frz.) may occur [-]
	131	real(kind=real8) :: qw__OK ! Flag: 1 => Non-zero cloud droplets Concentration
	132	real(kind=real8) :: qi__OK ! Flag: 1 => Non-zero cloud ice Concentration
	133	real(kind=real8) :: CCNiOK ! Flag: 1 => Non-zero cloud ice Particles
	134	real(kind=real8) :: CCNiIc ! [Nb/m3]
	135
	136	real(kind=real8) :: Flag_TmaxHM ! Flag: T < TmaxHM => Nucleation II (Hall-Mossop) may occur [-]
	137	real(kind=real8) :: Flag_TminHM ! Flag: T > TminHM => Nucleation II (Hall-Mossop) may occur [-]
	138	real(kind=real8) :: Flag_wa__HM ! Flag: w > wa__HM => Nucleation II (Hall-Mossop) may occur [-]
	139	real(kind=real8) :: SplinJ ! Hallet-Mossop Theory (Levkov et al., 1992,
	140	real(kind=real8) :: SplinP ! Contr.Atm.Phy. 65, p.40)
	141
	142	real(kind=real8) :: Flag_Ta_Neg ! Flag: T < 0.0dgC [-]
	143	real(kind=real8) :: Flag_TqwFrz ! Flag: T < Temper. => Instantaneous Freezing may occur [-]
	144
	145	real(kind=real8) :: RHumid ! Relative Humidity [-]
	146
	147	real(kind=real8) :: qi_Nu1,qi_Nu2 ! Nucleation ( 0 > Ta > -35 dgC ), Generations 1 and 2 [kg/kg]
	148	real(kind=real8) :: qi_Nuc ! Nucleation ( 0 > Ta > -35 dgC ), Generation (effective) [kg/kg]
	149	real(kind=real8) :: NuIdOK ! Nucleation I: Contact -Freez. Nucl. [-]
	150	real(kind=real8) :: BSPRWI ! Nucleation I: Contact -Freez. Nucl. [kg/kg]
	151	real(kind=real8) :: BHAMWI ! Nucleation II: Hallett-Mossop Ice-Multiplic. Process [kg/kg]
	152	real(kind=real8) :: BNUFWI ! Heterogeneous Freezing of Cloud Droplets [kg/kg]
	153	real(kind=real8) :: BDEPVI ! Ice Crystals Sublim. [kg/kg]
	154	real(kind=real8) :: Flag_SURSat ! Flag = 1/0 if (SUR/sub)Saturation FLAG [-]
	155	real(kind=real8) :: Flag_SUBSat ! Flag = 1/0 if (SUB/sur)Saturation FLAG [-]
	156	real(kind=real8) :: Flag_Sublim ! Flag = 1/0 for Sublimation Occurence or not [-]
	157	real(kind=real8) :: RH_Ice ! Relative Humidity vs Ice [-]
	158	real(kind=real8) :: RH_Liq ! Relative Humidity vs Water [-]
	159	real(kind=real8) :: DenDi1,DenDi2 ! dqi/dt\|sublimation: terms of the denominator [...]
	160	real(kind=real8) :: dqsiqv ! SUBsaturation qsi - qv [kg/kg]
	161	real(kind=real8) :: dqvSUB ! SURsaturation qv - qsiEFF [kg/kg]
	162	real(kind=real8) :: dqvDUM ! Water Vapor Variation, dummy variable [kg/kg]
	163	real(kind=real8) :: dqiDUM ! Ice Crystals Variation, dummy variable [kg/kg]
	164
	165	real(kind=real8) :: qCloud ! qw + qi [kg/kg]
	166	real(kind=real8) :: coefC2 ! Coefficient of Ek and Mahrt parameter C2_EkM [../..]
	167	real(kind=real8) :: pa_hPa,es_hPa ! Pressure, Pressure of Vapor at Saturat., over Water [hPa]
	168	real(kind=real8) :: Qsat_L ! Specific Concentration of Vapor at Saturat., over Water [kg/kg]
	169	! (even for temperatures smaller than freezing pt)
	170	real(kind=real8) :: t_qvqw ! qv+qw mixing ratio used in Partial Condensation Scheme [kg/kg]
	171	real(kind=real8) :: d_qvqw ! qv+qw mixing ratio variation [kg/kg]
	172	real(kind=real8) :: Kdqvqw ! qv+qw vertical turbulent Flux [kg/kg m/s]
	173	real(kind=real8) :: ww_TKE ! vertical Velocity Variance [m2/s2]
	174	real(kind=real8) :: RH_TKE ! Relative Humidity Variance [../..]
	175	real(kind=real8) :: qt_TKE ! qv+qw Variance [../..]
	176	real(kind=real8) :: TLiqid ! Liquid Temperature [K]
	177	real(kind=real8) :: CFr_t1,CFr_t2 ! [-]
	178	real(kind=real8) :: CFrCoe ! [-]
	179	real(kind=real8) :: CFraOK ! [-]
	180	real(kind=real8) :: qwCFra ! CloudAveraged Liquid Water Mixing Ratio [kg/kg]
	181	real(kind=real8) :: qwMesh ! Mesh Averaged Liquid Water Mixing Ratio [kg/kg]
	182	real(kind=real8) :: dwMesh ! Mesh Averaged Liquid Water Mixing Ratio Variation [kg/kg]
	183	real(kind=real8) :: signdw ! Sign of Liquid Water Mixing Ratio Variation (-1/1) [-]
	184	real(kind=real8) :: Flag_dqwPos ! Flag = 1/0 if Liquid Water Mixing Ratio Variation >/< 0 [-]
	185	real(kind=real8) :: updatw ! [-]
	186	real(kind=real8) :: SCuLim ! Fraction Limit [-]
	187	real(kind=real8) :: ARGerf ! Argument of erf function (used in partial condensation Scheme)
	188	real(kind=real8) :: OUTerf ! OUTPUT of erf function
	189	real(kind=real8) :: erf ! erf function
	190	real(kind=real8) :: dwTUR4,dwTURi !
	191	real(kind=real8) :: dwTUR3,dwTUR2 !
	192	real(kind=real8) :: dwTUR8,dwTURc !
	193	real(kind=real8) :: dwTUR5,dwTUR1 !
	194
	195	real(kind=real8) :: Di_Pri ! Pristine Ice Diameter
	196	real(kind=real8) :: c1saut,cnsaut ! Ice Crystals AUToconv. Parameters
	197	real(kind=real8) :: dtsaut ! Ice Crystals AUToconv. Time Scale
	198	real(kind=real8) :: ps_AUT ! Ice Crystals AUToconv. (BDEPIS, BAGRIS,...) Rate [kg/kg/s]
	199	real(kind=real8) :: qs_AUT ! Ice Crystals AUToconv. (BDEPIS, BAGRIS,...) [kg/kg]
	200
	201	real(kind=real8) :: Flag_Ta_Pos ! Flag = 1/0 if Ta >/< 273.15 K FLAG [-]
	202	real(kind=real8) :: Flag_qiMELT ! Flag = 1/0 for qi Melting or not FLAG [-]
	203	real(kind=real8) :: qxMelt ! Potential qi to Melt [kg/kg]
	204	real(kind=real8) :: qiMELT ! Effective qi to Melt [kg/kg]
	205	real(kind=real8) :: CiMelt ! Effective CCNi removed by qi Melting [nb/m3]
	206
	207	real(kind=real8) :: Flag_qsMELT ! Flag = 1/0 for qs Melting or not FLAG [-]
	208	real(kind=real8) :: xCoefM !
	209	real(kind=real8) :: AcoefM,BcoefM !
	210	real(kind=real8) :: dTMELT ! Temperature Offset before Snow Particles Melting [K]
	211	real(kind=real8) :: qsMELT ! Melt of Snow Particles [kg/kg]
	212
	213	real(kind=real8) :: qs_ACW ! Accretion of Cloud Drop. by Snow Particl. [kg/kg]
	214	real(kind=real8) :: Flag_qs_ACW ! Accretion of Cloud Drop. by Snow Particl. FLAG [-]
	215	real(kind=real8) :: Flag_qs_ACI ! Accretion of Cloud Ice by Snow Particl. FLAG [-]
	216	real(kind=real8) :: effACI ! Accretion of Cloud Ice by Snow Particl. Efficiency [-]
	217	real(kind=real8) :: ps_ACI ! Accretion of Cloud Ice by Snow Particl. Rate [kg/kg/s]
	218	real(kind=real8) :: qs_ACI ! Accretion of Cloud Ice by Snow Particl. [kg/kg]
	219	real(kind=real8) :: CNsACI ! Accretion of Cloud Ice by Snow Particl. CCNi decrease [nb/m3]
	220	real(kind=real8) :: Flag_qs_ACR ! Accretion of Snow by Rain FLAG [-]
	221	real(kind=real8) :: coeACR ! Accretion of Snow by Rain Sedimentation Coeffic. [...]
	222	real(kind=real8) :: qs_ACR ! Accretion of Snow by Rain [kg/kg]
	223	real(kind=real8) :: qs_ACR_R ! Accretion of Snow by Rain => Rain [kg/kg]
	224	real(kind=real8) :: qs_ACR_S ! => Graupels [kg/kg]
	225	real(kind=real8) :: Flag_qr_ACS ! Accretion of Rain by Snow Particl. => Snow FLAG [-]
	226	real(kind=real8) :: coeACS ! Accretion of Rain by Snow Sedimentation Coeffic. [...]
	227	real(kind=real8) :: pr_ACS ! Accretion of Rain by Snow Particl. => Snow Rate [kg/kg/s]
	228	real(kind=real8) :: qr_ACS ! Accretion of Rain by Snow Particl. => Snow [kg/kg]
	229	real(kind=real8) :: qr_ACS_S ! Accretion of Rain by Snow Particl. => Snow [kg/kg]
	230	real(kind=real8) :: ps_SUB ! Deposition/Sublim. on/of Snow Particl. Rate [kg/kg/s]
	231	real(kind=real8) :: qs_SUB ! Deposition/Sublim. on/of Snow Particl. [kg/kg]
	232	real(kind=real8) :: ls_NUM ! Sublimation of Snow: NUMerator of Snow Distribut.Coef. [...]
	233	real(kind=real8) :: ls_DEN ! Sublimation of Snow: DENominator of Snow Distribut.Coef. [...]
	234
	235	real(kind=real8) :: Flag_Freeze ! Freezing of Rain FLAG [-]
	236	real(kind=real8) :: ps_FRZ ! Freezing of Rain Rate [kg/kg/s]
	237	real(kind=real8) :: qs_FRZ ! Freezing of Rain [kg/kg]
	238
	239	real(kind=real8) :: rwMEAN ! Droplets Autoconversion: Mean Radius
	240
	241	real(kind=real8) :: th_AUT ! Cloud Droplets AUToconv. Threshold [-]
	242	real(kind=real8) :: pr_AUT ! Cloud Droplets AUToconv. Rate [kg/kg/s]
	243	real(kind=real8) :: qr_AUT ! Cloud Droplets AUToconv. [kg/kg]
	244	real(kind=real8) :: Flag_qr_ACW ! Accretion of Cloud Drop. by Rain FLAG [-]
	245	real(kind=real8) :: pr_ACW ! Accretion of Cloud Drop. by Rain Rate [kg/kg/s]
	246	real(kind=real8) :: qr_ACW ! Accretion of Cloud Drop. by Rain [kg/kg]
	247	real(kind=real8) :: Flag_qr_ACI ! Accretion of Cloud Drop. by Rain => Rain FLAG [-]
	248	real(kind=real8) :: pr_ACI ! Accretion of Cloud Ice by Rain => Rain Rate [kg/kg/s]
	249	real(kind=real8) :: qr_ACI ! Accretion of Cloud Ice by Rain => Rain [kg/kg]
	250	real(kind=real8) :: CNrACI ! Accretion of Cloud Ice by Rain CCNi decrease [nb/m3]
	251	real(kind=real8) :: pi_ACR ! Accretion of Cloud Ice by Rain => Snow Rate [kg/kg/s]
	252	real(kind=real8) :: qi_ACR ! Accretion of Cloud Ice by Rain => Snow [kg/kg]
	253	real(kind=real8) :: Flag_DryAir ! 1 => RH_Liq < 1 FLAG [-]
	254	real(kind=real8) :: Flag_qr_EVP ! Evaporation of Rain FLAG [-]
	255	real(kind=real8) :: lr_NUM ! Evaporation of Rain: NUMerator of rain Distribut.Coef. [...]
	256	real(kind=real8) :: lr_DEN ! Evaporation of Rain: DENominator of rain Distribut.Coef. [...]
	257	real(kind=real8) :: pr_EVP ! Evaporation of Rain Rate [kg/kg/s]
	258	real(kind=real8) :: qr_EVP ! Evaporation of Rain [kg/kg]
	259
	260	real(kind=real8) :: effACS ! Accretion of Snow by Graupels Efficiency [-]
	261
	262	real(kind=real8) :: a_rodz ! Air Mass
	263	real(kind=real8) :: qwFlux ! qw Sedimentation Flux Coefficient
	264	real(kind=real8) :: qwrodz ! Droplets Mass
	265	real(kind=real8) :: wRatio ! Droplets Mass Ratio [-]
	266	real(kind=real8) :: qiFlux ! qi Sedimentation Flux Coefficient
	267	real(kind=real8) :: qirodz ! Crystals Mass
	268	real(kind=real8) :: iRatio ! Crystals Mass Ratio [-]
	269	real(kind=real8) :: qsFlux ! qs Sedimentation Flux Coefficient
	270	real(kind=real8) :: qsrodz ! Snow Mass
	271	real(kind=real8) :: sRatio ! Snow Mass Ratio [-]
	272	real(kind=real8) :: qrFlux ! qr Sedimentation Flux Coefficient
	273	real(kind=real8) :: qrrodz ! Rain Mass
	274	real(kind=real8) :: rRatio ! Rain Mass Ratio [-]
	275
	276	real(kind=real8) :: Vw_MAX ! MAX Sedimentation Velocity of Droplets
	277	real(kind=real8) :: Flag_Fall_i ! Sedimentation of Ice FLAG [-]
	278	real(kind=real8) :: Vi_MAX ! MAX Sedimentation Velocity of Ice Particles
	279	real(kind=real8) :: Vs_MAX,VsMMAX ! MAX Sedimentation Velocity of Snow Particles
	280	real(kind=real8) :: Vs__OK ! Sedimentation Velocity of Snow Particles
	281	real(kind=real8) :: Vr_MAX,VrMMAX ! MAX Sedimentation Velocity of Rain Drops
	282	real(kind=real8) :: Vr__OK ! Sedimentation Velocity of Rain Drops
	283	real(kind=real8) :: dtPrec ! Sedimentation Time Step
	284	real(kind=real8) :: d_Snow ! Sedimented Snow Part. over 1 time step [m]
	285	real(kind=real8) :: d_Rain ! Sedimented Rain Drops over 1 time step [m]
	286
	287	real(kind=real8) :: SatiOK !
	288	real(kind=real8) :: FlagNu ! 1 =>Nucleation for -35 dgC < Ta < 0 dgc
	289
	290	real(kind=real8) :: Qw0_OK ! HYDROMETEORS INPUT STATUS
	291	real(kind=real8) :: Qi0_OK !
	292	real(kind=real8) :: Qi0qOK !
	293	real(kind=real8) :: Ci0_OK !
	294	real(kind=real8) :: Ci0cOK !
	295	real(kind=real8) :: Qs0_OK !
	296	real(kind=real8) :: Qr0_OK !
	297
	298	real(kind=real8) :: qiBerg ! Bergeron-Findeisen Process: avaiklable qi
	299	real(kind=real8) :: qwBerg ! Bergeron-Findeisen Process: avaiklable qw
	300	real(kind=real8) :: qxBerg ! Bergeron-Findeisen Process: potential qi Generation
	301	real(kind=real8) :: a1Berg,a2Berg ! Bergeron-Findeisen Process: parameters
	302	real(kind=real8) :: a0Berg,afBerg ! Bergeron-Findeisen Process: integration constants
	303
	304	real(kind=real8) :: WaterB ! Vertically Integrated Water Budget
	305
	306	real(kind=real8) :: argEXP !
	307
	308	integer :: k ,ikl !
	309	integer :: i_Berg !
	310	integer :: nItMAX,itFall !
	311	integer :: ikl_io,io__Pt !
	312
	313
	314	! Debug Variables
	315	! ---------------
	316
	317	! #wH character(len=70) :: debugH !
	318	! #wH character(len=10) :: proc_1,proc_2 !
	319	! #wH character(len=10) :: proc_3,proc_4 !
	320	! #wH real(kind=real8) :: procv1,procv2 !
	321	! #wH real(kind=real8) :: procv3,procv4 !
	322	! #wH integer :: kv ,nl !
	323
	324
	325
	326
	327	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	328	! !
	329	! ALLOCATION !
	330	! ========== !
	331
	332	IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN !
	333
	334	allocate ( qw_io0( mzp) ) ! Droplets Concentration entering CMiPhy [kg/kg]
	335	allocate ( qi_io0( mzp) ) ! Ice Part. Concentration entering CMiPhy [kg/kg]
	336	allocate ( qs___0(kcolp,mzp) ) ! Snow Part. Concentration entering CMiPhy [kg/kg]
	337	! #qg allocate ( qg___0(kcolp,mzp) ) ! Graupels Concentration entering CMiPhy [kg/kg]
	338	allocate ( qr___0(kcolp,mzp) ) ! Rain Drops Concentration entering CMiPhy [kg/kg]
	339	allocate ( Ta_dgC(kcolp,mzp) ) ! Air Temperature [dgC]
	340	allocate ( sqrrro(kcolp,mzp) ) ! sqrt(roa(mzp)/roa(k)) [-]
	341	allocate ( qsiEFF(kcolp,mzp) ) ! EFFective Saturation Specific Humidity over Ice [kg/kg]
	342	allocate ( Fletch(kcolp,mzp) ) ! Monodisperse Nb of hexagonal Plates, Fletcher (1962) [-]
	343	allocate ( lamdaS(kcolp,mzp) ) ! Marshall-Palmer distribution parameter for Snow Particl.
	344	! #qg allocate ( lamdaG(kcolp,mzp) ) ! Marshall-Palmer distribution parameter for Graupels
	345	allocate ( lamdaR(kcolp,mzp) ) ! Marshall-Palmer distribution parameter for Rain Drops
	346	allocate ( ps_ACR(kcolp,mzp) ) ! Accretion of Snow by Rain Rate [kg/kg/s]
	347	allocate ( ps_ACW(kcolp,mzp) ) ! Accretion of Cloud Drop. by Snow Particl. Rate [kg/kg/s]
	348	allocate ( FallVw(kcolp,mzp) ) ! Sedimentation Velocity of Droplets
	349	allocate ( FallVi(kcolp,mzp) ) ! Sedimentation Velocity of Ice Particles
	350	allocate ( FallVs(kcolp,mzp) ) ! Sedimentation Velocity of Snow Particles
	351	! #qg allocate ( FallVg(kcolp,mzp) ) ! Sedimentation Velocity of Snow Particles
	352	allocate ( FallVr(kcolp,mzp) ) ! Sedimentation Velocity of Rain Drops
	353	allocate ( qwLoss( mzp) ) ! Mass Loss related to Sedimentation of Rain Droplets
	354	allocate ( qiLoss( mzp) ) ! Mass Loss related to Sedimentation of Ice Crystals
	355	allocate ( qsLoss( mzp) ) ! Mass Loss related to Sedimentation of Snow Particles
	356	allocate ( qrLoss( mzp) ) ! Mass Loss related to Sedimentation of Rain Drops
	357	! #wH allocate ( debugV( mzp,16) ) ! Debug Variable (of 16 microphysical processes)
	358
	359	! #WH allocate ( wihm1(mzp) ) ! Cloud Droplets Freezing
	360	! #WH allocate ( wihm2(mzp) ) ! Ice Crystals Homogeneous Sublimation
	361	! #WH allocate ( wicnd(mzp) ) ! Ice Crystals Nucleation (Emde & Kahlig)
	362	! #WH allocate ( widep(mzp) ) ! Ice Crystals Growth Bergeron Process (Emde & Kahlig)
	363	! #WH allocate ( wisub(mzp) ) ! Ice Crystals Sublimation (Levkov)
	364	! #WH allocate ( wimlt(mzp) ) ! Ice Crystals Melting
	365	! #WH allocate ( wwevp(mzp) ) ! Water Vapor Condensation / Evaporation (Fractional Cloudiness)
	366	! #WH allocate ( wraut(mzp) ) ! Cloud Droplets AUTO-Conversion
	367	! #WH allocate ( wsaut(mzp) ) ! Ice Crystals AUTO-Conversion
	368	! #WH allocate ( wracw(mzp) ) ! Accretion of Cloud Droplets by Rain, Ta > 0, --> Rain
	369	! #WH allocate ( wsacw(mzp) ) ! Accretion of Cloud Droplets by Rain, Ta < 0, --> Snow
	370	! #WH allocate ( wsaci(mzp) ) ! Accretion of Ice Crystals by Snow --> Snow
	371	! #WH allocate ( wraci(mzp) ) ! Accretion of Ice Crystals by Rain --> Snow
	372	! #WH allocate ( wiacr(mzp) ) ! Accretion of Rain by Ice Crystals --> Snow
	373	! #WH allocate ( wsacr(mzp) ) ! Accretion of Rain by Snow --> Snow
	374	! #WH allocate ( wracs(mzp) ) ! Accretion of Snow by Rain --> Snow, Rain
	375	! #WH allocate ( wrevp(mzp) ) ! Rain Drops Evaporation
	376	! #WH allocate ( wssub(mzp) ) ! Snow Particles Sublimation
	377	! #WH allocate ( wsmlt(mzp) ) ! Snow Particles Melting
	378	! #WH allocate ( wsfre(mzp) ) ! Rain Drops Freezing
	379
	380	END IF !
	381	! !
	382	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	383
	384
	385
	386
	387	! +++++++++++++++++++
	388	DO ikl=1,kcolp
	389	! +++++++++++++++++++
	390
	391
	392
	393
	394	! Debug
	395	! ~~~~~
	396	! #wH DO k=mz1_CM,mzp
	397	! #wH debugH( 1:35) = 'CMiPhy: Debugged Variables: Initial'
	398	! #wH debugH(36:70) = ' '
	399	! #wH proc_1 = 'R.Hum W[%]'
	400	! #wH procv1 = 0.1qv__DY(ikl,k)/(RHcrit qvswCM(ikl,k))
	401	! #wH proc_2 = 'R.Hum I[%]'
	402	! #wH procv2 = 0.1qv__DY(ikl,k)/(RHcrit qvsiCM(ikl,k))
	403	! #wH proc_3 = ' '
	404	! #wH procv3 = 0.
	405	! #wH proc_4 = ' '
	406	! #wH procv4 = 0.
	407
	408	! #wh include 'CMiPhy_Debug.h'
	409
	410	! #wH DO kv=1,16
	411	! #wH debugV(k,kv) = 0.
	412	! #wH ENDDO
	413	! #wH END DO
	414
	415
	416
	417
	418	! Vertical Integrated Energy and Water Content
	419	! ============================================
	420
	421	! #EW enr0EW(ikl) = 0.0
	422	! #EW wat0EW(ikl) = 0.0
	423
	424	! #EW DO k=1,mzp
	425	! #EW enr0EW(ikl) = enr0EW( ikl) &
	426	! #EW& +(Ta__CM(ikl,k) &
	427	! #EW& -(qw__CM(ikl,k)+qr__CM(ikl,k))*Lv_Cpd &
	428	! #EW& -(qi__CM(ikl,k)+qs__CM(ikl,k))*Ls_Cpd) &
	429	! #EW& * dsigmi(k)
	430	! #EW wat0EW(ikl) = wat0EW( ikl) &
	431	! #EW& +(qv__DY(ikl,k) &
	432	! #EW& + qw__CM(ikl,k)+qr__CM(ikl,k) &
	433	! #EW& + qi__CM(ikl,k)+qs__CM(ikl,k) ) &
	434	! #EW& * dsigmi(k)
	435	! #EW END DO
	436
	437	! #EW mphyEW(ikl) =' '
	438	! .. mphy2D --> '12345678901234567890'
	439
	440	! #ew enr0EW(ikl) = enr0EW(ikl) * psa_DY(ikl) * Grav_I
	441	! #EW wat0EW(ikl) = wat0EW(ikl) * psa_DY(ikl) * Grav_I
	442	! .. wat0EW [m] contains an implicit factor 1.d3 [kPa-->Pa] /ro_Wat
	443
	444
	445	! #WH VsMMAX = 0.0
	446	! #WH VrMMAX = 0.0
	447
	448
	449
	450
	451	! Set lower limit on Hydrometeor Concentration
	452	! ============================================
	453
	454	! #hy IF (NO_Vec) THEN
	455
	456	! #hy DO k=mz1_CM,mzp
	457
	458	! #hy IF (qw__CM(ikl,k).lt.epsn) THEN
	459	! #hy qv__DY(ikl,k) = qv__DY(ikl,k)+qw__CM(ikl,k)
	460	! #hy Ta__CM(ikl,k) = Ta__CM(ikl,k)-qw__CM(ikl,k)*Lv_Cpd
	461	! #hy qwd_CM(ikl,k)= qwd_CM(ikl,k)-qw__CM(ikl,k)
	462	! #hy qw__CM(ikl,k) = 0.0
	463	! #hy END IF
	464
	465	! #hy IF (qr__CM(ikl,k).lt.epsn) THEN
	466	! #hy qv__DY(ikl,k) = qv__DY(ikl,k)+qr__CM(ikl,k)
	467	! #hy Ta__CM(ikl,k) = Ta__CM(ikl,k)-qr__CM(ikl,k)*Lv_Cpd
	468	! #hy qwd_CM(ikl,k) = qwd_CM(ikl,k)-qr__CM(ikl,k)
	469	! #hy qr__CM(ikl,k) = 0.0
	470	! #hy END IF
	471
	472	! #hy IF (qi__CM(ikl,k).lt.epsn.or.CCNiCM(ikl,k).lt.un_1) THEN
	473	! #hy qv__DY(ikl,k) = qv__DY(ikl,k)+qi__CM(ikl,k)
	474	! #hy Ta__CM(ikl,k) = Ta__CM(ikl,k)-qi__CM(ikl,k)*Ls_Cpd
	475	! #hy qid_CM(ikl,k) = qid_CM(ikl,k)-qi__CM(ikl,k)
	476	! #hy qi__CM(ikl,k) = 0.0
	477	! #hy CCNiCM(ikl,k) = 0.0
	478	! #hy END IF
	479
	480	! #hy IF (qs__CM(ikl,k).lt.epsn) THEN
	481	! #hy qv__DY(ikl,k) = qv__DY(ikl,k)+qs__CM(ikl,k)
	482	! #hy Ta__CM(ikl,k) = Ta__CM(ikl,k)-qs__CM(ikl,k)*Ls_Cpd
	483	! #hy qid_CM(ikl,k) = qid_CM(ikl,k)-qs__CM(ikl,k)
	484	! #hy qs__CM(ikl,k) = 0.0
	485	! #hy END IF
	486	! #hy END DO
	487
	488	! #hy ELSE
	489
	490	DO k=mz1_CM,mzp
	491
	492	Qw0_OK = max(zer0,sign(un_1,epsn-qw__CM(ikl,k)))*qw__CM(ikl,k)
	493	qw__CM(ikl,k) = qw__CM(ikl,k) -Qw0_OK
	494	qwd_CM(ikl,k) = qwd_CM(ikl,k) -Qw0_OK
	495	qv__DY(ikl,k) = qv__DY(ikl,k) +Qw0_OK
	496	Ta__CM(ikl,k) = Ta__CM(ikl,k) -Qw0_OK*Lv_Cpd
	497
	498	Qr0_OK = max(zer0,sign(un_1,epsn-qr__CM(ikl,k)))*qr__CM(ikl,k)
	499	qr__CM(ikl,k) = qr__CM(ikl,k) -Qr0_OK
	500	qwd_CM(ikl,k) = qwd_CM(ikl,k) -Qr0_OK
	501	qv__DY(ikl,k) = qv__DY(ikl,k) +Qr0_OK
	502	Ta__CM(ikl,k) = Ta__CM(ikl,k) -Qr0_OK*Lv_Cpd
	503
	504	Qi0qOK = max(zer0,sign(un_1,epsn-qi__CM(ikl,k)))
	505	Ci0cOK = max(zer0,sign(un_1,un_1-CCNiCM(ikl,k)))
	506
	507	Ci0_OK = max(Ci0cOK,Qi0qOK)
	508	Qi0_OK = Ci0_OK*qi__CM(ikl,k)
	509
	510	CCNiCM(ikl,k) = Ci0_OK*CCNiCM(ikl,k)
	511	qi__CM(ikl,k) = qi__CM(ikl,k) - Qi0_OK
	512	qid_CM(ikl,k) = qid_CM(ikl,k) - Qi0_OK
	513	qv__DY(ikl,k) = qv__DY(ikl,k) + Qi0_OK
	514	Ta__CM(ikl,k) = Ta__CM(ikl,k) - Qi0_OK*Ls_Cpd
	515
	516	Qs0_OK = max(zer0,sign(un_1,epsn-qs__CM(ikl,k)))*qs__CM(ikl,k)
	517	qs__CM(ikl,k) = qs__CM(ikl,k) -Qs0_OK
	518	qid_CM(ikl,k) = qid_CM(ikl,k) -Qs0_OK
	519	qv__DY(ikl,k) = qv__DY(ikl,k) +Qs0_OK
	520	Ta__CM(ikl,k) = Ta__CM(ikl,k) -Qs0_OK*Ls_Cpd
	521
	522	END DO
	523
	524	! #hy END IF
	525
	526
	527
	528
	529	! Initial Concentrations
	530	! ======================
	531
	532	DO k=1,mzp
	533	Ta_dgC(ikl,k) = Ta__CM(ikl,k) -Tf_Sno
	534	Fletch(ikl,k) = 1.e-2exp(-0.6Ta_dgC(ikl,k)) ! Ice Crystals Number (Fletcher, 1962)
	535
	536	qr___0(ikl,k) = qr__CM(ikl,k)
	537	qs___0(ikl,k) = qs__CM(ikl,k)
	538	! #qg qg___0(ikl,k) = qg__CM(ikl,k)
	539
	540	! #WH IF (ikl.eq.ikl0CM(1)) THEN
	541	! #WH qw_io0(k) = qw__CM(ikl,k)
	542	! #WH qi_io0(k) = qi__CM(ikl,k)
	543	! #WH END IF
	544
	545	END DO
	546
	547
	548
	549
	550	! Saturation Specific Humidity
	551	! ============================
	552
	553	DO k=1,mzp
	554
	555	qsiEFF(ikl,k) = RHcrit * qvsiCM(ikl,k ) ! Saturation Specific Humidity over Ice
	556
	557	sqrrro(ikl,k) = sqrt((psa_DY(ikl )+pt__DY) &!
	558	& /(roa_DY(ikl,k )*R_DAir &!
	559	& * Ta__CM(ikl,mzp))) !
	560
	561
	562
	563
	564	! Hydrometeors Fall Velocities
	565	! ==============================
	566
	567	! Cloud Droplets Fall Velocity (Calcul de la Vitesse Terminale Moyenne) FALL VELOCITY
	568	! ----------------------------
	569
	570	! #VW IF (qw__CM(ikl,k).ge.epsn) THEN
	571
	572	! #VW CCNwCM(ikl,k) = 1.2d+8 ! ASTEX case (Duynkerke et al. 1995, JAS 52, p.2763)
	573
	574	! #VW qwCFra = qw__CM(ikl,k) / max(CFrMIN ,CFraCM(ikl,k))
	575	! #VW dwTUR4 = 4.5 qwTURB qwTURB
	576	! #VW dwTUR1 = 12.5 qwTURB qwTURB
	577	! #VW dwTURi = qwCFra * roa_DY(ikl,k) &
	578	! #VW& * 6.0d+0/(piNmbrCCNwCM(ikl,k)exp(dwTUR4))
	579	! #VW dwTUR5 = exp(R_5by3*log(dwTURi))
	580
	581	! #VW FallVw(ikl,k) = 1.19d8* piNmbrCCNwCM(ikl,k) dwTUR5 &
	582	! #VW& * exp(dwTUR1)/(24.0 roa_DY(ikl,k) qwCFra)
	583	! #VW ELSE
	584	! #VW FallVw(ikl,k) = 0.00
	585	! #VW END IF
	586
	587
	588	! Rain Fall Velocity FALL VELOCITY
	589	! ----------------------------
	590
	591	lamdaR(ikl,k) = exp(0.25log((piNmbrn0___r) &! Marshall-Palmer Distribution Parameter for Rain
	592	& / (roa_DY(ikl,k)*max( epsn ,qr__CM(ikl,k))))) ! Ref.: Emde and Kahlig 1989, Ann.Geoph. 7, p.407 (3)
	593	! Note that a simplification occurs
	594	! between the 1000. factor of rho, and rho_water=1000.
	595
	596	! #hy IF (qr__CM(ikl,k).gt.epsn) THEN
	597
	598	Vr__OK = max(zer0,sign(un_1, qr__CM(ikl,k) - epsn)) ! Vr__OK = 1.0 if qr__CM(ikl,k) > epsn
	599	! = 0.0 otherwise
	600
	601	FallVr(ikl,k) = Vr__OK392. sqrrro(ikl,k) &! Terminal Fall Velocity for Rain
	602	& / exp(0.8 *log(lamdaR(ikl,k))) ! 392 = a Gamma[4+b] / 6
	603	! where a = 842. b = 0.8
	604
	605	! #hy ELSE
	606	! #hy FallVr(ikl,k)= 0.0
	607	! #hy END IF
	608
	609
	610	! Snow Fall Velocity (c and d parameters: see Locatelli and Hobbs, 1974, JGR: table 1 p.2188) FALL VELOCITY
	611	! ------------------
	612
	613	! #cn n0___s = min(2.e8 &
	614	! #cn& ,2.e6exp(-.12min(0.,Ta_dgC(ikl,k))))
	615
	616	lamdaS(ikl,k) = exp(0.25log((0.50piNmbr*n0___s) &! Marshall-Palmer distribution parameter for Snow
	617	& / (roa_DY(ikl,k)*max(epsn,qs__CM(ikl,k)))))! Ref.: Emde and Kahlig 1989, Ann.Geoph. 7, p.407 (3)
	618	! Levkov et al. 1992, Cont.Atm.Phys. 65(1) p.37 (5) (rho_snow)
	619	! Note that a partial simplification occurs
	620	! between the 1000. factor of rho, and rho_snow=500.
	621
	622	! #hy IF (qs__CM(ikl,k).gt.epsn) THEN
	623
	624	Vs__OK = max(zer0,sign(un_1, qs__CM(ikl,k) - epsn)) ! Vs__OK = 1.0 if qs__CM(ikl,k) > epsn
	625	! = 0.0 otherwise
	626
	627	IF (graupel_shape) THEN
	628	FallVs(ikl,k) = Vs__OK2.19 sqrrro(ikl,k) &! Terminal Fall Velocity for Graupellike Snow Flakes of Hexagonal Type
	629	& / exp(0.25*log(lamdaS(ikl,k))) ! 2.19 = c Gamma[4+d] / 6
	630	! where c = 4.836, d = 0.25
	631	! = 0.86 1000.*0.25
	632	ELSE IF (planes__shape) THEN
	633	FallVs(ikl,k) = Vs__OK2976.sqrrro(ikl,k) &! Terminal Fall Velocity for Unrimed Side Planes
	634	& / exp(0.99*log(lamdaS(ikl,k))) ! 2976.= c Gamma[4+d] / 6
	635	! where c = 755.9, d = 0.99
	636	! = 0.81 1000.*0.99
	637
	638	ELSE IF (aggrega_shape) THEN
	639	FallVs(ikl,k) = Vs__OK20.06sqrrro(ikl,k) &! Terminal Fall Velocity for Aggregates of unrimed radiating assemblages
	640	& / exp(0.41*log(lamdaS(ikl,k))) ! 2976.= c Gamma[4+d] / 6
	641	! where c = 755.9, d = 0.41
	642	! = 0.69 1000.*0.41
	643	ELSE
	644	STOP 'Snow Particles Shape is not defined'
	645	END IF
	646
	647	! #hy ELSE
	648	! #hy FallVs(ikl,k) = 0.0 ! FALL VELOCITY
	649	! #hy END IF
	650
	651
	652	! Graupel Fall Velocity
	653	! ---------------------
	654
	655	! #qg IF (qg__CM(ikl,k).ge.epsn) THEN
	656	! #qg lamdaG(ikl,k) =exp(0.250log((piNmbrn0___g) &! Marshall-Palmer distribution parameter for Graupel
	657	! #qg& /(roa_DY(ikl,k)*max(epsn ,qg__CM(ikl,k))))) ! Note that a simplification occurs
	658	! between the 1000. factor of rho, and rho_ice=1000.
	659	! #qg FallVg(ikl,k) = 25.1 *sqrrro(ikl,k) &! 25.1 = c Gamma[4+d] / 6
	660	! #qg& / exp(0.57log(lamdaG(ikl,k))) ! where c = 4.836 = 1.10 1000.**0.57 and d = 0.57
	661	! ! Hexagonal Graupel, Locatelli and Hobbs, 1974, JGR: table 1 p.2188:
	662
	663	! #qg ELSE
	664	! #qg FallVg(ikl,k) = 0.0
	665	! #qg lamdaG(ikl,k) = 0.0
	666	! #qg END IF
	667
	668	END DO
	669
	670
	671	!=============================================================================== CLOUD ICE PARTICLES
	672	! ++++++++++++++++++++
	673	! Microphysical Processes affecting non Precipitating Cloud Particles
	674	! ===================================================================
	675
	676	DO k=mz1_CM,mzp
	677
	678
	679	! Homogeneous Nucleation by Cloud Dropplets Solidification ! BFREWI
	680	! Reference: Emde and Kahlig 1989, Ann.Geoph. 7, p.407 (11) ! Levkov (24) p.40
	681	! ---------------------------------------------------------
	682
	683	! #wH Flag_Ta_Neg= 0.
	684	! #wH qw__OK = 0.
	685
	686	! #hy IF (Ta_dgC(ikl,k).lt.TqwFrz)THEN
	687
	688	Flag_TqwFrz=max(zer0,-sign(un_1,Ta_dgC(ikl,k) - TqwFrz)) ! Flag_TqwFrz = 1.0 if Ta_dgC(ikl,k) < TqwFrz
	689	! = 0.0 otherwise
	690
	691	! #EW IF(Flag_TqwFrz.gt.eps6) THEN
	692	! #EW mauxEW = mphyEW(ikl)
	693	! #EW mauxEW(01:01) = 'i'
	694	! #EW mphyEW(ikl) = mauxEW
	695	! #EW END IF
	696
	697	qw__OK = qw__CM(ikl,k) * Flag_TqwFrz
	698	qi__CM(ikl,k) = qi__CM(ikl,k) + qw__OK
	699	CCNiCM(ikl,k) = CCNiCM(ikl,k) + roa_DY(ikl,k) *qw__OK/qw_VOL
	700	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Lc_Cpd *qw__OK
	701
	702	! #WQ write(6,*) 'Qihm1',qw__CM(ikl,k), &
	703	! #WQ& ' Qi' ,qi__CM(ikl,k), &
	704	! #WQ& ' CcnI' ,CCNiCM(ikl,k),it_EXP,ikl,k
	705	! #WH IF (ikl.eq.ikl0CM(1)) wihm1(k) = qw__OK
	706
	707	qw__CM(ikl,k) = qw__CM(ikl,k) - qw__OK
	708
	709	! #hy END IF
	710
	711
	712	! Heterogeneous Freezing of Cloud Droplets ! BNUFWI CLOUD ICE PARTICLES
	713	! Reference: Levkov et al., 1992 (21) p.40 ! Levkov (21) p.40 ++++++++++++++++++++
	714	! ----------------------------------------
	715
	716	IF (Heter_Freezng) THEN
	717	! #hy IF (Ta_dgC(ikl,k).lt.0.0) THEN
	718
	719	Flag_Ta_Neg=max(zer0,-sign(un_1,Ta_dgC(ikl,k) - 0.0)) ! Flag_Ta_Neg = 1.0 if Ta_dgC(ikl,k) < 0.00dgC
	720	! = 0.0 otherwise
	721
	722	argEXP = min(max(ea_MIN , -Ta_dgC(ikl,k)) ,ea_MAX)
	723	BNUFWI = Flag_Ta_Neg*(exp(argEXP) -1. ) &!
	724	& * 100.* qw__CM(ikl,k) *qw_VOL
	725	BNUFWI = min( BNUFWI , qw__CM(ikl,k) )
	726
	727	qi__CM(ikl,k) = qi__CM(ikl,k) + BNUFWI
	728	CCNiCM(ikl,k) = CCNiCM(ikl,k) + roa_DY(ikl,k)*BNUFWI/qw_VOL
	729	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Lc_Cpd *BNUFWI
	730	qw__CM(ikl,k) = qw__CM(ikl,k) - BNUFWI
	731
	732
	733	! Debug
	734	! ~~~~~
	735	! #wH debugH( 1:35) = 'Homo+Hetero Nucleation by Droplets '
	736	! #wH debugH(36:70) = 'Solidification (BFREWI+BNUFWI) '
	737	! #wH proc_1 = 'BFREWI '
	738	! #wH procv1 = Flag_TqwFrz
	739	! #wH proc_2 = 'BNUFWI '
	740	! #wH procv2 = BNUFWI
	741	! #wH proc_3 = ' '
	742	! #wH procv3 = 0.
	743	! #wH proc_4 = ' '
	744	! #wH procv4 = 0.
	745	! #wh include 'CMiPhy_Debug.h'
	746	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	747	! #wH& debugV(k,01) = qw__OK+BNUFWI
	748
	749	! #hy END IF
	750	END IF
	751
	752
	753	!=============================================================================== CLOUD ICE PARTICLES
	754	! ++++++++++++++++++++
	755	! Homogeneous Sublimation ! XXXXXX
	756	! Reference: Emde and Kahlig 1989, Ann.Geoph. 7, p.407 (12) ! Levkov
	757	! ---------------------------------------------------------
	758
	759	! #EW IF (Flag_TqwFrz.gt.eps6) THEN
	760	! #EW mauxEW = mphyEW(ikl)
	761	! #EW mauxEW(02:02) = 'I'
	762	! #EW mphyEW(ikl) = mauxEW
	763	! #EW END IF
	764
	765	IF (Homog_Sublima) THEN
	766	dqvSUB = (qv__DY(ikl,k)-qsiEFF(ikl,k)) &!
	767	& /(1.00 +1.733e7qsiEFF(ikl,k) &! 1.733e7=LsLs*0.622/Cpa/Ra with Ls = 2833600 J/kg
	768	& /(Ta__CM(ikl,k)*Ta__CM(ikl,k))) !
	769
	770	dqvSUB = Flag_TqwFrz*max(zer0,dqvSUB)
	771	dqvDUM = dqvSUB
	772
	773	qi__CM(ikl,k) = qi__CM(ikl,k) + dqvDUM
	774	qid_CM(ikl,k) = qid_CM(ikl,k) + dqvDUM
	775	! CCNiCM(ikl,k) : NO VARIATION
	776	qv__DY(ikl,k) = qv__DY(ikl,k) - dqvDUM
	777	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Ls_Cpd * dqvDUM
	778
	779	! Full Debug
	780	! ~~~~~~~~~~
	781	! #WQ write(6,*) 'Qihm2',dqvDUM, &
	782	! #WQ& ' Qi' ,qi__CM(ikl,k), &
	783	! #WQ& ' CcnI' ,CCNiCM(ikl,k),it_EXP,ikl,k
	784	! #WH if (ikl.eq.ikl0CM(1)) wihm2(k) = dqvDUM
	785
	786	! Debug
	787	! ~~~~~
	788	! #wH debugH( 1:35) = 'Emde and Kahlig: Homogeneous Sublim'
	789	! #wH debugH(36:70) = 'ation '
	790	! #wH proc_1 = 'dQv g/kg'
	791	! #wH procv1 = dqvDUM
	792	! #wH proc_2 = ' '
	793	! #wH procv2 = 0.
	794	! #wH proc_3 = ' '
	795	! #wH procv3 = 0.
	796	! #wH proc_4 = 'CCNI/1.e15'
	797	! #wH procv4 = CCNiCM(ikl,k)*1.e-18
	798	! #wh include 'CMiPhy_Debug.h'
	799	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	800	! #wH& debugV(k,01) = dqvDUM + debugV(k,01)
	801
	802	END IF
	803	END DO
	804
	805
	806
	807
	808	!=============================================================================== CLOUD ICE PARTICLES, MEYERS
	809	! ++++++++++++++++++++
	810	! Nucleation I: Deposition & Condensation-Freezing Nucleat.
	811	! Source : Water Vapor ! BNUCVI
	812	! Reference: Meyers et al., 1992, JAM 31, (2.4) p.712 ! Levkov (20) p.40
	813	! -----------------------------------------------------------
	814
	815	IF (Meyers) THEN
	816	DO k=mz1_CM,mzp
	817
	818	! #wH NuIdOK = 0.
	819	! #wH CCNiId = 0.
	820	! #wH BNUCVI = 0.
	821	! #wH BSPRWI = 0.
	822	! #wH BHAMWI = 0.
	823
	824	! #hy IF (Ta_dgC(ikl,k).lt.T_NuId) THEN
	825	Flag_T_NuId=max(zer0,-sign(un_1,Ta_dgC(ikl,k) - T_NuId))! Flag_T_NuId = 1.0 if Ta_dgC(ikl,k) < T_NuId
	826	! ! = 0.0 otherwise
	827
	828	dqvDUM = qv__DY(ikl,k) - qsiEFF(ikl,k) ! Sursaturat.
	829
	830	! #hy IF (dqvDUM.gt.0.) THEN
	831	SatiOK = max(zer0,sign(un_1,dqvDUM)) ! SatiOK = 1.0 if qv__DY(ikl,k) > qsiEFF
	832	! ! = 0.0 otherwise
	833	dqvDUM = max(zer0, dqvDUM)
	834
	835	NuIdOK = Flag_T_NuId * SatiOK
	836
	837	SSat_I = 1.e2*dqvDUM / qsiEFF(ikl,k) ! Sursaturat.%I
	838	SSat_I = min(SSat_I,SSImax) !
	839	CCNiId = 1.0e3 * exp(a_NuId+b_NuId*SSat_I) ! Meyers et al. 1992 JAM, 2.4
	840	CCNiId = max(CCNiId-CCNiCM(ikl,k),zer0) &!
	841	& * NuIdOK !
	842	CCNiCM(ikl,k) = CCNiId+CCNiCM(ikl,k) !
	843	dqiDUM = 1.e-15* CCNiId/roa_DY(ikl,k) ! 1.e-15 = 0.001 * Initial Ice Crystal Mass
	844	dqiDUM = min(dqiDUM , dqvDUM)
	845	qi__CM(ikl,k) = qi__CM(ikl,k) + dqiDUM
	846	qid_CM(ikl,k) = qid_CM(ikl,k) + dqiDUM
	847	qv__DY(ikl,k) = qv__DY(ikl,k) - dqiDUM
	848	Ta__CM(ikl,k) = Ta__CM(ikl,k) + dqiDUM*Ls_Cpd
	849	BNUCVI = dqiDUM
	850
	851	! #hy END IF
	852	! #hy END IF
	853
	854
	855	! Nucleation I: Contact -Freezing Nucleat. CLOUD ICE PARTICLES, MEYERS
	856	! Source : Cloud Dropplets ! BSPRWI ++++++++++++++++++++
	857	! Reference: Meyers et al., 1992, JAM 31, (2.6) p.713 ! Levkov (20) p.40
	858	! -----------------------------------------------------------
	859
	860	! #wH CCNiIc = 0.
	861	! #wH dqiDUM = 0.
	862
	863	! #hy IF (qw__CM(ikl,k).gt.0.) THEN
	864	qw__OK = max(zer0,sign(un_1,qw__CM(ikl,k))) ! qw__OK = 1.0 if qw__CM(ikl,k) > 0.
	865	! = 0.0 otherwise
	866
	867	! #hy IF (Ta_dgC(ikl,k).lt.T_NuIc) THEN
	868	Flag_T_NuIc = max(zer0,-sign(un_1,Ta_dgC(ikl,k) - T_NuIc))
	869	! Flag_T_NuIc = 1.0 if Ta_dgC(ikl,k) < T_NuIc
	870	! = 0.0 otherwise
	871
	872	Flag___NuIc = Flag_T_NuIc * qw__OK
	873
	874	CCNiIc = 1.e3 * Flag___NuIc &! Contact-Freez
	875	& * exp(a_NuIc &! Potent.Nuclei
	876	& -b_NuIc &! Meyers et al.
	877	& *Ta_dgC(ikl,k)) ! 1992 JAM, 2.6
	878	rad_ww = (1.e3 * roa_DY(ikl,k) &! Drop. Radius
	879	& * qw__CM(ikl,k) &!
	880	& * .2e-11 ) ** 0.33 !
	881	! .2e-11 = 1. / (1.2e+8 * 1.e3 * 4.19)
	882	! CCNwCM (ASTEX) ro_w 4 pi /3
	883	CCNiIc = 603.2e+3 * CCNiIc * rad_ww &! Levkov et al.
	884	& * roa_DY(ikl,k) ! 1992 CAM,(23)
	885	! 603.2e3 = 4.0e-7 * 4 pi * 1.2e+8 * 1.e3
	886	! DFar CCNwCM fact(rolv)
	887	CCNiCM(ikl,k) = CCNiCM(ikl,k) &!
	888	& + CCNiIc !
	889	dqiDUM = 1.e-15 * CCNiIc/roa_DY(ikl,k)
	890	! 1.e-15 = 1.0e-3 * Ice Crystal Mass
	891	dqiDUM = min( qw__CM(ikl,k) , dqiDUM)
	892	qi__CM(ikl,k) = qi__CM(ikl,k) + dqiDUM
	893	qw__CM(ikl,k) = qw__CM(ikl,k) - dqiDUM
	894	Ta__CM(ikl,k) = Ta__CM(ikl,k) + dqiDUM*Lc_Cpd
	895	BSPRWI = dqiDUM
	896
	897	! #hy END IF
	898	! #hy END IF
	899
	900
	901	! Nucleation II: Hallett-Mossop Ice-Multiplication Process ! BSPRWI CLOUD ICE PARTICLES, MEYERS
	902	! Reference: Levkov et al., 1992, Contr.Atm.Ph.65,(25) p.40 ! Levkov (25) p.40 ++++++++++++++++++++
	903	! -----------------------------------------------------------
	904
	905	IF (HalMos) THEN
	906	! #hy IF (Ta_dgC(ikl,k).lt.TmaxHM.AND. &
	907	! #hy& Ta_dgC(ikl,k).gt.TminHM.AND. &
	908	! #hy& wa__DY(ikl,k).gt.wa__HM ) THEN
	909	Flag_TmaxHM = max(zer0,-sign(un_1,Ta_dgC(ikl,k) - TmaxHM))
	910	! Flag_TmaxHM = 1.0 if Ta_dgC(ikl,k) < TmaxHM
	911	! = 0.0 otherwise
	912
	913	Flag_TminHM = max(zer0, sign(un_1,Ta_dgC(ikl,k) - TminHM))
	914	! Flag_TminHM = 1.0 if Ta_dgC(ikl,k) > TminHM
	915	! = 0.0 otherwise
	916
	917	Flag_wa__HM = max(zer0, sign(un_1,wa__DY(ikl,k) - wa__HM))
	918	! Flag_wa__HM = 1.0 if wa__DY(ikl,k) > wa__HM
	919	! = 0.0 otherwise
	920
	921	! #cn n0___s = min(2.e8,2.e6exp(-.12min(Ta_dgC(ikl,k),0.)))
	922
	923	SplinJ = 1.358e12 *qw__CM(ikl,k) &!
	924	& n0___s /(lamdaS(ikl,k)*.33)
	925	! 1.358e12=pi Gamma(3.5) g ro_s /(3 Cd *4.19e-12)
	926	! [=3.14 3.3233625 9.810.1 /(3 0.6 *4.19e-12)]
	927	SplinP = 0.003 * (1. - 0.05 SplinJ) Flag_TmaxHM &!
	928	& * Flag_TminHM * Flag_wa__HM !
	929	SplinP = max(zer0, SplinP)
	930
	931	dqiDUM = 1.e-15 * SplinP / roa_DY(ikl,k) ! 1.e-15 = 1.0e-3 * Ice Crystal Mass
	932	SplinP = (min(1.0,qs__CM(ikl,k)/max(dqiDUM,epsn))) *SplinP
	933	CCNiCM(ikl,k) = CCNiCM(ikl,k) +SplinP
	934	dqiDUM = min(qs__CM(ikl,k), dqiDUM)
	935	qi__CM(ikl,k) = qi__CM(ikl,k) + dqiDUM
	936	qid_CM(ikl,k) = qid_CM(ikl,k) + dqiDUM
	937	qs__CM(ikl,k) = qs__CM(ikl,k) - dqiDUM
	938	BHAMWI = dqiDUM
	939	! #hy END IF
	940	END IF
	941
	942
	943	! Debug
	944	! ~~~~~
	945	! #wH debugH( 1:35) = 'Meyers: Nucl. I, Depot & Cond-Freez'
	946	! #wH debugH(36:70) = 'Nucl. / Freez / Nucl. II / Bergeron'
	947	! #wH proc_1 = 'dQi1 Meyer'
	948	! #wH procv1 = BNUCVI
	949	! #wH proc_2 = 'dQi2 Meyer'
	950	! #wH procv2 = BSPRWI
	951	! #wH proc_3 = 'dQi Ha-Mos'
	952	! #wH procv3 = BHAMWI
	953	! #wH proc_4 = ' '
	954	! #wH procv4 = 0.
	955	! #wh include 'CMiPhy_Debug.h'
	956	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	957	! #wH& debugV(k,02) = BNUCVI + BSPRWI + BHAMWI
	958
	959	END DO
	960
	961
	962
	963
	964	!===============================================================================
	965
	966	ELSE
	967
	968	!=============================================================================== CLOUD ICE PARTICLES, EMDE & KAHLIG
	969	! ++++++++++++++++++++
	970	! Ice Crystals Nucleation Process between 0.C and -35.C
	971	! (each crystal has a mass equal or less than 10d-12 kg)
	972	! Reference: Emde and Kahlig 1989, Ann.Geoph. 7, p.408 (13)
	973	! ---------------------------------------------------------
	974
	975	DO k=mz1_CM,mzp
	976
	977	! #wH qi_Nu1 = 0.
	978	! #wH qi_Nu2 = 0.
	979	! #wH qi_Nuc = 0.
	980
	981	! #hy IF (Ta_dgC(ikl,k).gt.TqwFrz) THEN
	982	Flag_TqwFrz = max(zer0, sign(un_1,Ta_dgC(ikl,k) - TqwFrz))
	983	! Flag_TqwFrz = 1.0 if Ta_dgC(ikl,k) > TqwFrz
	984	! = 0.0 otherwise
	985
	986	! #hy IF (Ta_dgC(ikl,k).lt.0.e0 ) THEN
	987	Flag_Ta_Neg = max(zer0,-sign(un_1,Ta_dgC(ikl,k) ))
	988	! Flag_Ta_Neg = 1.0 if Ta_dgC(ikl,k) < 0.e0
	989	! = 0.0 otherwise
	990
	991	! #hy IF (qv__DY(ikl,k).gt.qsiEFF(ikl,k)) THEN
	992	SatiOK = max(zer0, sign(un_1,qv__DY(ikl,k) - qsiEFF(ikl,k)))
	993	! SatiOK = 1.0 if qv__DY(ikl,k) > qsiEFF(ikl,k)
	994	! = 0.0 otherwise
	995
	996	FlagNu = Flag_TqwFrz * Flag_Ta_Neg * SatiOK
	997
	998	! #EW IF(FlagNu.gt.eps6) THEN
	999	! #EW mauxEW = mphyEW(ikl)
	1000	! #EW mauxEW(03:03) = 'I'
	1001	! #EW mphyEW(ikl) = mauxEW
	1002	! #EW END IF
	1003
	1004	qi_Nu1 = FlagNu * 1.d-15 * Fletch(ikl,k) /roa_DY(ikl,k)
	1005	! qi_Nu1 : amount of nucleated ice crystals (first condition)
	1006
	1007	qi_Nu1 = qi_Nu1*max(zer0,sign(un_1,qi_Nu1-qi__CM(ikl,k)))
	1008
	1009	qi_Nu2 = ( qv__DY(ikl,k)-qsiEFF(ikl,k)) &
	1010	& /(1.0d0+1.733d7*qsiEFF(ikl,k) &
	1011	& /(Ta__CM(ikl,k)*Ta__CM(ikl,k)))
	1012	qi_Nu2 = FlagNu * max(zer0 ,qi_Nu2) ! amount of nucleated ice crystals (second condition)
	1013
	1014	qi_Nuc = min(qi_Nu1,qi_Nu2)
	1015
	1016	qi__CM(ikl,k) = qi__CM(ikl,k) + qi_Nuc
	1017	qid_CM(ikl,k) = qid_CM(ikl,k) + qi_Nuc
	1018	CCNiCM(ikl,k) = CCNiCM(ikl,k) + roa_DY(ikl,k) *qi_Nuc &
	1019	& *1.e15
	1020	qv__DY(ikl,k) = qv__DY(ikl,k) - qi_Nuc
	1021	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Ls_Cpd *qi_Nuc
	1022
	1023	! Full Debug
	1024	! ~~~~~~~~~~
	1025	! #WQ write(6,*) 'QiCnd',qi_Nuc, &
	1026	! #WQ& ' Qi' ,qi__CM(ikl,k), &
	1027	! #WQ& ' CcnI' ,CCNiCM(ikl,k),it_EXP,ikl,k
	1028	! #WH IF (ikl.eq.ikl0CM(1)) wicnd(k) = qi_Nuc
	1029
	1030	! Debug
	1031	! ~~~~~
	1032	! #wH debugH( 1:35) = 'Emde and Kahlig: Ice Crystals Nucle'
	1033	! #wH debugH(36:70) = 'ation Process between 0.C and -35.C'
	1034	! #wH proc_1 = 'Qicnd1 '
	1035	! #wH procv1 = qi_Nu1
	1036	! #wH proc_2 = 'Qicnd2 '
	1037	! #wH procv2 = qi_Nu2
	1038	! #wH proc_3 = 'Qicnd g/kg'
	1039	! #wH procv3 = qi_Nuc
	1040	! #wH proc_4 = ' '
	1041	! #wH procv4 = 0.
	1042	! #wh include 'CMiPhy_Debug.h'
	1043	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1))&
	1044	! #wH& debugV(k,02) = qi_Nuc
	1045
	1046
	1047	! #hy END IF
	1048	! #hy END IF
	1049	! #hy END IF
	1050
	1051
	1052	END DO
	1053
	1054	END IF
	1055
	1056
	1057
	1058
	1059	!============================================================================== CLOUD PARTICLES, MIXED PHASE
	1060	! +++++++++++++++
	1061	! Bergeron Process (water vapor diffusion-deposition on ice crystals)
	1062	! Reference: Koenig 1971, J.A.S. 28, p.235
	1063	! Emde and Kahlig 1989, Ann.Geoph. 7, p.408 (14)
	1064	! ---------------------------------------------------------
	1065
	1066	IF (.NOT.AUTO_i_LevkXX) THEN
	1067
	1068	DO k=mz1_CM,mzp
	1069
	1070	! #wH qi0_OK = 0.
	1071	! #wH qxBerg = 0.
	1072	! #wH qwBerg = 0.
	1073
	1074	! #hy IF (qi__CM(ikl,k).gt.epsn
	1075	! #hy& .AND. Ta_dgC(ikl,k).lt.0.e0) THEN
	1076
	1077	qi0_OK = max(zer0, sign(un_1,qi__CM(ikl,k) - epsn))
	1078	! qi0_OK = 1.0 if qi__CM(ikl,k) > epsn
	1079	! = 0.0 otherwise
	1080
	1081	Flag_Ta_Neg = max(zer0,-sign(un_1,Ta_dgC(ikl,k) ))
	1082	! Flag_Ta_Neg = 1.0 if Ta_dgC(ikl,k) < 0.e0
	1083	! = 0.0 otherwise
	1084
	1085	qi0_OK = Flag_Ta_Neg * qi0_OK
	1086
	1087	! #EW IF(qi0_OK.gt.eps6) THEN
	1088	! #EW mauxEW = mphyEW(ikl)
	1089	! #EW mauxEW(04:04) = 'i'
	1090	! #EW mphyEW(ikl) = mauxEW
	1091	! #EW END IF
	1092
	1093	i_Berg = abs(Ta_dgC(ikl,k)-un_1)
	1094	i_Berg = min(i_Berg,31)
	1095	i_Berg = max(i_Berg, 1)
	1096	a1Berg = aa1(i_Berg)
	1097	a2Berg = aa2(i_Berg)
	1098
	1099	a0Berg = 1.d+3roa_DY(ikl,k)qi__CM(ikl,k) / Fletch(ikl,k)
	1100	afBerg =(a1Berg (1.0-a2Berg) dt__CM &! analytical integration of (14) p.408
	1101	& +a0Berg(1.0-a2Berg))(1.0/(1.0-a2Berg)) ! Emde and Kahlig 1989, Ann.Geoph. 7
	1102	qxBerg =(1.d-3*Fletch(ikl,k)/roa_DY(ikl,k)) &!
	1103	& *(afBerg-a0Berg) !
	1104	qxBerg = max(zer0,qxBerg)
	1105
	1106	qwBerg = max(zer0,qw__CM(ikl,k)) ! qwBerg : to avoid the use of qwd_CM < 0.
	1107
	1108	qxBerg = qi0_OK*min(qwBerg,qxBerg)
	1109	qi__CM(ikl,k)= qi__CM(ikl,k) +qxBerg
	1110	! CCNiCM(ikl,k):NO VARIATION
	1111
	1112	qw__CM(ikl,k)= qw__CM(ikl,k) -qxBerg
	1113	Ta__CM(ikl,k)= Ta__CM(ikl,k)+Lc_Cpd *qxBerg
	1114
	1115	! Full Debug
	1116	! ~~~~~~~~~~
	1117	! #WQ write(6,*) 'QiDep',qxBerg,
	1118	! #WQ& ' Qi' ,qi__CM(ikl,k),
	1119	! #WQ& ' CcnI' ,CCNiCM(ikl,k),it_EXP,ikl,k
	1120	! #WH IF (ikl.eq.ikl0CM(1)) widep(k)= qxBerg
	1121
	1122	! Debug
	1123	! ~~~~~
	1124	! #wH debugH( 1:35) = 'Bergeron Process (water vapor diffu'
	1125	! #wH debugH(36:70) = 'sion-deposition on ice crystals) '
	1126	! #wH proc_1 = 'qi0_OK ICE'
	1127	! #wH procv1 = qi0_OK
	1128	! #wH proc_2 = 'Qicnd g/kg'
	1129	! #wH procv2 = qwBerg
	1130	! #wH proc_3 = 'Qidep g/kg'
	1131	! #wH procv3 = qxBerg
	1132	! #wH proc_4 = ' '
	1133	! #wH procv4 = 0.
	1134	! #wh include 'CMiPhy_Debug.h'
	1135	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1136	! #wH& debugV(k,02) = qxBerg + debugV(k,02)
	1137
	1138
	1139	! #hy END IF
	1140
	1141	END DO
	1142
	1143	END IF
	1144
	1145
	1146
	1147
	1148	!=============================================================================== CLOUD ICE PARTICLES
	1149
	1150	! Ice Crystals Sublimation ! BDEPVI
	1151	! Reference: Emde and Kahlig, 1989 p.408 (15) ! Levkov (27) p.40
	1152	! -------------------------------------------
	1153
	1154	DO k=mz1_CM,mzp
	1155
	1156	! #wH BDEPVI = 0.
	1157
	1158	! #hy IF (qsiEFF(ikl,k).gt.qv__DY(ikl,k)) THEN
	1159
	1160	dqsiqv = qsiEFF(ikl,k) - qv__DY(ikl,k)
	1161	! #pp Flag_SUBSat = max(zer0,sign(un_1,dqsiqv))
	1162	! Flag_SUBSat = 1.0 if qsiEFF(ikl,k) > qv__DY(ikl,k)
	1163	! = 0.0 otherwise
	1164
	1165	! #hy IF (qi__CM(ikl,k).gt.epsn) THEN
	1166
	1167	qi0_OK = max(zer0,sign(un_1, qi__CM(ikl,k) - epsn))
	1168	! qi0_OK = 1.0 if qi__CM(ikl,k) > epsn
	1169	! = 0.0 otherwise
	1170
	1171	Flag_Sublim = qi0_OK &
	1172	! #pp& * Flag_SUBSat &
	1173	& + 0.0
	1174
	1175	! #EW IF(Flag_Sublim.gt.eps6) THEN
	1176	! #EW mauxEW = mphyEW(ikl)
	1177	! #EW mauxEW(05:05) = 'V'
	1178	! #EW mphyEW(ikl) = mauxEW
	1179	! #EW END IF
	1180
	1181	RH_Ice = qv__DY(ikl,k) / qsiEFF(ikl,k)
	1182	DenDi1 = 6.959d+11 / (Ta__CM(ikl,k)*Ta__CM(ikl,k))
	1183	! 6.959e+11
	1184	! = [Ls=2833600J/kg] * Ls / [kT=0.025W/m/K] / [Rv=461.J/kg/K]
	1185	! kT: Air thermal Conductivity
	1186	DenDi2 = 1.0d0 / (1.875d-2roa_DY(ikl,k)qsiEFF(ikl,k))
	1187	! 1.875d-5: Water Vapor Diffusivity in Air
	1188	BDEPVI = dt__CM (1.-RH_Ice)4.0Di_Hex Fletch(ikl,k) &!
	1189	& /(DenDi1+DenDi2)
	1190	BDEPVI = max(BDEPVI, -qv__DY(ikl,k)) ! H2O deposit.limit = H2O content
	1191	BDEPVI = min(BDEPVI, qi__CM(ikl,k)) ! qi sublim. limit = qi content
	1192	BDEPVI = min(BDEPVI, dqsiqv ) &! qi sublim. limit = Saturation
	1193	& * Flag_Sublim
	1194
	1195	qi__CM(ikl,k) = qi__CM(ikl,k) - BDEPVI
	1196	qid_CM(ikl,k) = qid_CM(ikl,k) - BDEPVI
	1197	qv__DY(ikl,k) = qv__DY(ikl,k) + BDEPVI
	1198	Ta__CM(ikl,k) = Ta__CM(ikl,k) - Ls_Cpd *BDEPVI
	1199
	1200	! Full Debug
	1201	! ~~~~~~~~~~
	1202	! #WQ write(6,*) 'QiSub', BDEPVI, &
	1203	! #WQ& ' Qi' , qi__CM(ikl,k), &
	1204	! #WQ& ' CcnI' ,CCNiCM(ikl,k),it_EXP,ikl,k
	1205	! #WH IF (ikl.eq.ikl0CM(1)) wisub(k) = BDEPVI
	1206
	1207	! #hy END IF
	1208	! #hy END IF
	1209
	1210	! Debug
	1211	! ~~~~~
	1212	! #wH debugH( 1:35) = 'Emde and Kahlig: Ice Crystals Subli'
	1213	! #wH debugH(36:70) = 'mation '
	1214	! #wH proc_1 = 'Qisub g/kg'
	1215	! #wH procv1 = BDEPVI
	1216	! #wH proc_2 = 'R.Hum I[%]'
	1217	! #wH procv2 = 0.1 * RH_Ice
	1218	! #wH proc_3 = ' '
	1219	! #wH procv3 = 0.
	1220	! #wH proc_4 = ' '
	1221	! #wH procv4 = 0.
	1222	! #wh include 'CMiPhy_Debug.h'
	1223	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1224	! #wH& debugV(k,03) = -BDEPVI
	1225
	1226	END DO
	1227
	1228	DO k=mz1_CM,mzp
	1229	IF (qi__CM(ikl,k).le.0.e0) THEN
	1230	qi__CM(ikl,k) = 0.e0
	1231	CCNiCM(ikl,k) = 0.e0
	1232	END IF
	1233	END DO
	1234
	1235
	1236
	1237
	1238	!===============================================================================
	1239
	1240	! Ice Crystals Instantaneous Melting
	1241	! ----------------------------------
	1242
	1243	DO k=mz1_CM,mzp
	1244
	1245	! #wH qiMELT = 0.
	1246	! #wH CiMelt = 0.
	1247
	1248	! #hy IF (Ta_dgC(ikl,k).gt.0.e0)THEN
	1249
	1250	Flag_Ta_Pos = max(zer0, sign(un_1,Ta_dgC(ikl,k) ))
	1251	! Flag_Ta_Pos = 1.0 if Ta_dgC(ikl,k) > 0.e0
	1252	! = 0.0 otherwise
	1253
	1254	! #hy IF (qi__CM(ikl,k).gt.epsn)THEN
	1255	qi0_OK = max(zer0, sign(un_1,qi__CM(ikl,k) - epsn))
	1256	! qi0_OK = 1.0 if qi__CM(ikl,k) > epsn
	1257	! = 0.0 otherwise
	1258
	1259	Flag_qiMELT = Flag_Ta_Pos * qi0_OK
	1260
	1261	! #EW IF(Flag_qiMELT .gt.eps6) THEN
	1262	! #EW mauxEW = mphyEW(ikl)
	1263	! #EW mauxEW(06:06) = 'w'
	1264	! #EW mphyEW(ikl) = mauxEW
	1265	! #EW END IF
	1266
	1267	qxMelt = Ta_dgC(ikl,k) / Lc_Cpd
	1268	qiMELT = min(qi__CM(ikl,k) , qxMelt)*Flag_qiMELT
	1269	CiMelt = CCNiCM(ikl,k) * qiMELT &
	1270	& /max(qi__CM(ikl,k) , epsn)
	1271	qi__CM(ikl,k) = qi__CM(ikl,k) - qiMELT
	1272	CCNiCM(ikl,k) = CCNiCM(ikl,k) - CiMelt
	1273	qw__CM(ikl,k) = qw__CM(ikl,k) + qiMELT
	1274	Ta__CM(ikl,k) = Ta__CM(ikl,k) - Lc_Cpd *qiMELT
	1275
	1276	! Full Debug
	1277	! ~~~~~~~~~~
	1278	! #WQ write(6,*) 'QiMlt',qiMELT, &
	1279	! #WQ& ' Qi' ,qi__CM(ikl,k), &
	1280	! #WQ& ' CcnI' ,CCNiCM(ikl,k),it_EXP,ikl,k
	1281	! #WH IF (ikl.eq.ikl0CM(1)) wimlt(k) = qiMELT
	1282
	1283	! #hy END IF
	1284	! #hy END IF
	1285
	1286	! Debug
	1287	! ~~~~~
	1288	! #wH debugH( 1:35) = 'Emde and Kahlig: Ice Crystals Insta'
	1289	! #wH debugH(36:70) = 'ntaneous Melting '
	1290	! #wH proc_1 = 'Qimlt g/kg'
	1291	! #wH procv1 = qiMELT
	1292	! #wH proc_2 = 'CiMelt /e15'
	1293	! #wH procv2 = CiMelt*1.e-18
	1294	! #wH proc_3 = ' '
	1295	! #wH procv3 = 0.
	1296	! #wH proc_4 = ' '
	1297	! #wH procv4 = 0.
	1298	! #wh include 'CMiPhy_Debug.h'
	1299	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1300	! #wH& debugV(k,04) = -qiMELT
	1301
	1302	END DO
	1303
	1304
	1305
	1306
	1307	!=============================================================================== CONDENSATION, Delobbe SCu
	1308	! +++++++++++++++++++++++++
	1309	! Water Vapor Condensation / Evaporation (Fractional Cloudiness)
	1310	! Reference: Laurent Delobbe Thesis (Ek &Mahrt 1991)
	1311	! --------------------------------------------------------------
	1312
	1313	DO k=mz1_CM,mzp ! Zeroing needed since
	1314	CFraCM(ikl,k) = 0.0 ! a maximization process
	1315	END DO
	1316
	1317	IF (Frac__Clouds.AND.fracSC) THEN
	1318
	1319	DO k=mz1_CM,mzp
	1320
	1321	! #wH dwMesh = 0.
	1322
	1323	! #hy IF (Ta_dgC(ikl,k).ge.TqwFrz) THEN
	1324
	1325	Flag_TqwFrz = max(zer0,sign(un_1,Ta_dgC(ikl,k) - TqwFrz))
	1326	! Flag_TqwFrz = 1.0 if Ta_dgC(ikl,k) > TqwFrz
	1327	! = 0.0 otherwise
	1328
	1329	t_qvqw = qv__DY(ikl,k) + qw__CM(ikl,k) ! Total Water Mixing Ratio
	1330	TLiqid = Ta__CM(ikl,k) - Lv_Cpd * qw__CM(ikl,k) ! Liquid Temperature
	1331
	1332	! Saturation specific humidity over water,
	1333	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ corresponding to liquid temperature
	1334	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1335	pa_hPa =(psa_DY(ikl) * sigma(k) + pt__DY) * 10.0d0 ! Dudhia (1989) JAS, (B1) and (B2) p.3103
	1336	es_hPa = 6.1078d0 * exp (ExpWat* log(WatIce /TLiqid)) &! (see also Pielke (1984), p.234, and
	1337	& * exp (ExpWa2*(un_1/WatIce-un_1/TLiqid)) ! Stull (1988), p.276
	1338
	1339	Qsat_L = .622d0es_hPa /(pa_hPa - .378d0es_hPa) ! Saturation Vapor Specific Concentration over Water
	1340	! (even for temperatures less than freezing point)
	1341
	1342	! Partial Condensation/Scheme
	1343	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1344	d_qvqw = qv__DY(ikl,MIN(k+1,mzp))-qv__DY(ikl,k) &
	1345	& + qwd_CM(ikl,MIN(k+1,mzp))-qw__CM(ikl,k)
	1346	Kdqvqw = Kzh_AT(ikl,k)*d_qvqw/(Z___DY(ikl,k+1)-Z___DY(ikl,k))
	1347
	1348	ww_TKE = 0.66d+0 * TKE_AT(ikl,k) ! Vertical Velocity Variance
	1349
	1350	coefC2 = Kdqvqw/(sqrt(ww_TKE)*Qsat_L)
	1351	RH_TKE = C1_EkM + C2_EkM * coefC2 * coefC2 ! Relative Humidity Variance
	1352	! (Ek and Mahrt, 1991, An. Geoph., 9, 716--724)
	1353
	1354	qt_TKE = sqrt(RH_TKE)*Qsat_L ! Total Water Variance
	1355
	1356	ARGerf = (t_qvqw-Qsat_L)/(1.414d+0*qt_TKE)
	1357	OUTerf = erf(ARGerf)
	1358
	1359	CFraCM(ikl,k) = 0.5d+0 * (1.d+0 + OUTerf) ! Cloud Fraction
	1360
	1361	CFrCoe = 1.d+0/(1.d+0+1.349d7Qsat_L/(TLiqidTLiqid)) !
	1362	CFr_t1 = qt_TKE/sqrt(piNmbr+piNmbr) &!
	1363	& * exp(-min(ARGerf*ARGerf,ea_MAX)) !
	1364	CFr_t2 = CFraCM(ikl,k) *(t_qvqw-Qsat_L) !
	1365
	1366	CFraOK = max(zer0,sign(un_1,CFraCM(ikl,k) - CFrMIN)) ! CFraOK = 1.0 if CFraCM(ikl,k) > CFrMIN
	1367	! = 0.0 otherwise
	1368
	1369	CFraCM(ikl,k) = CFraCM(ikl,k) * CFraOK * Flag_TqwFrz !
	1370	qwMesh = CFrCoe * (CFr_t1+CFr_t2) * CFraOK ! Mesh Averaged Liquid Water Mixing Ratio
	1371	dwMesh = qwMesh - qw__CM(ikl,k)
	1372
	1373	! Vectorisation of the Atmospheric Water Update
	1374	! ~~~~~~~~~~~~~+-------------------------------------------------+
	1375	! \| if (dwMesh.gt.0.d0) then \|
	1376	! \| dwMesh = min(qv__DY(ikl,k), dwMesh) \|
	1377	! \| else \|
	1378	! \| dwMesh =-min(qw__CM(ikl,k),-dwMesh) \|
	1379	! \| end if \|
	1380	! +-------------------------------------------------+
	1381
	1382	signdw = sign(un_1,dwMesh)
	1383	Flag_dqwPos = max(zer0,signdw)
	1384	updatw = Flag_dqwPos * qv__DY(ikl,k) &!
	1385	& + (1.d0 -Flag_dqwPos)* qw__CM(ikl,k) !
	1386	! #kk SCuLim = exp(min(0.,300.-Ta__CM(ikl,k))) ! SCu Lim.
	1387	dwMesh = signdw min(updatw, signdwdwMesh) &!
	1388	& *Flag_TqwFrz &!
	1389	! #kk& * SCuLim &! SCu
	1390	& + 0.0
	1391	! #kk CFraCM(ikl,k) = CFraCM(ikl,k) * SCuLim ! Limitor
	1392
	1393	! Update of qv__DY, qw__CM and Ta__CM
	1394	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1395	qw__CM(ikl,k) = qw__CM(ikl,k) + dwMesh
	1396	qwd_CM(ikl,k) = qwd_CM(ikl,k) + dwMesh
	1397	qv__DY(ikl,k) = qv__DY(ikl,k) - dwMesh
	1398	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Lv_Cpd * dwMesh
	1399
	1400	! Full Debug
	1401	! ~~~~~~~~~~
	1402	! #WQ write(6,*) 'QwEvp',dwMesh,it_EXP,ikl,k
	1403	! #WH if (ikl.eq.ikl0CM(1)) wwevp(k) = dwMesh
	1404
	1405	! #EW IF(Ta_dgC(ikl,k).ge.TqwFrz) THEN
	1406	! #EW mauxEW = mphyEW(ikl)
	1407	! #EW mauxEW(07:07) = 'W'
	1408	! #EW mphyEW(ikl) = mauxEW
	1409	! #EW END IF
	1410
	1411	! #hy END IF
	1412
	1413	! Debug
	1414	! ~~~~~
	1415	! #wH debugH( 1:35) = 'Delobbe: Condensation '
	1416	! #wH debugH(36:70) = ' '
	1417	! #wH proc_1 = 'dQw g/kg'
	1418	! #wH procv1 = dwMesh
	1419	! #wH proc_2 = ' '
	1420	! #wH procv2 = 0.
	1421	! #wH proc_3 = ' '
	1422	! #wH procv3 = 0.
	1423	! #wH proc_4 = ' '
	1424	! #wH procv4 = 0.
	1425	! #wh include 'CMiPhy_Debug.h'
	1426	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1427	! #wH& debugV(k,05) = dwMesh
	1428
	1429	END DO
	1430
	1431
	1432
	1433
	1434	!=============================================================================== CONDENSATION, NO SCu
	1435	! ++++++++++++++++++++
	1436	! Water Vapor Condensation / Evaporation
	1437	! Reference: Emde and Kahlig 1989, Ann.Geoph. 7, p.407 (7)
	1438	! --------------------------------------------------------
	1439
	1440	ELSE
	1441
	1442	DO k=mz1_CM,mzp
	1443
	1444	! #wH dwMesh = 0.
	1445
	1446	! #hy IF (Ta_dgC(ikl,k).ge.TqwFrz) THEN
	1447	Flag_TqwFrz = max(zer0, sign(un_1,Ta_dgC(ikl,k) - TqwFrz))
	1448	! Flag_TqwFrz = 1.0 if Ta_dgC(ikl,k) > TqwFrz
	1449	! = 0.0 otherwise
	1450
	1451	dwMesh = (qv__DY(ikl,k) -qvswCM(ikl,k)*RHcrit) &
	1452	& / (1.0d0+1.349d7*qvswCM(ikl,k) &
	1453	& /(Ta__CM(ikl,k)*Ta__CM(ikl,k)))
	1454	! 1.349e7=LvLv0.622/Cpa/Ra with Lv = 2500000 J/kg
	1455
	1456	! Vectorisation of the Atmospheric Water Update
	1457	! ~~~~~~~~~~~~~+-------------------------------------------------+
	1458	! \| if (dwMesh.gt.0.d0) then \|
	1459	! \| dwMesh = min(qv__DY(ikl,k), dwMesh) \|
	1460	! \| else \|
	1461	! \| dwMesh =-min(qw__CM(ikl,k),-dwMesh) \|
	1462	! \| end if \|
	1463	! +-------------------------------------------------+
	1464
	1465	signdw = sign(un_1,dwMesh)
	1466	Flag_dqwPos = max(zer0,signdw)
	1467	updatw = Flag_dqwPos * qv__DY(ikl,k) &
	1468	& + (1.d0 -Flag_dqwPos)* qw__CM(ikl,k)
	1469	dwMesh = signdw min(updatw,signdwdwMesh) &
	1470	& *Flag_TqwFrz
	1471
	1472	! Update of qv__DY, qw__CM and Ta__CM
	1473	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1474	qw__CM(ikl,k) = qw__CM(ikl,k) + dwMesh
	1475	qwd_CM(ikl,k) = qwd_CM(ikl,k) + dwMesh
	1476	qv__DY(ikl,k) = qv__DY(ikl,k) - dwMesh
	1477	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Lv_Cpd * dwMesh
	1478	! [Ls=2500000J/kg]/[Cp=1004J/kg/K]=2490.04
	1479
	1480	! #EW IF(Ta_dgC(ikl,k).ge.TqwFrz) THEN
	1481	! #EW mauxEW = mphyEW(ikl)
	1482	! #EW mauxEW(07:07) = 'W'
	1483	! #EW mphyEW(ikl) = mauxEW
	1484	! #EW END IF
	1485
	1486	! Full Debug
	1487	! ~~~~~~~~~~
	1488	! #WQ write(6,*) 'QwEvp',dwMesh,it_EXP,ikl,k
	1489	! #WH if (ikl.eq.ikl0CM(1)) wwevp(k) = dwMesh
	1490
	1491	! #hy END IF
	1492
	1493	! Debug
	1494	! ~~~~~
	1495	! #wH debugH( 1:35) = 'Emde and Kahlig: Water Vapor Conden'
	1496	! #wH debugH(36:70) = 'sation / Evaporation '
	1497	! #wH proc_1 = 'dQw g/kg'
	1498	! #wH procv1 = dwMesh
	1499	! #wH proc_2 = ' '
	1500	! #wH procv2 = 0.
	1501	! #wH proc_3 = ' '
	1502	! #wH procv3 = 0.
	1503	! #wH proc_4 = ' '
	1504	! #wH procv4 = 0.
	1505	! #wh include 'CMiPhy_Debug.h'
	1506	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1507	! #wH& debugV(k,05) = dwMesh
	1508
	1509	END DO
	1510
	1511	END IF
	1512
	1513
	1514
	1515
	1516	!=============================================================================== CONDENSATION, SCu added AFTER
	1517	! +++++++++++++++++++++++++++++
	1518	! Fractional Cloudiness ! Guess may be computed (Ek&Mahrt91 fracSC=.T.)
	1519	! ====================== ! Final value computed below
	1520
	1521	! #sc IF (Frac__Clouds.AND..NOT.fracSC) THEN
	1522	IF (Frac__Clouds) THEN
	1523	IF(fraCEP) THEN ! ECMWF Large Scale Cloudiness
	1524	! ----------------------------
	1525	DO k=mz1_CM,mzp
	1526	CFraCM(ikl,k) = (qi__CM(ikl,k) + qw__CM(ikl,k) &!
	1527	& +qs__CM(ikl,k) * 0.33 &!
	1528	& * (1.-min(un_1,exp((Ta__CM(ikl,k) -258.15)*0.1))))&!
	1529	& / (0.02 * qvswCM(ikl,k) ) !
	1530	CFraCM(ikl,k) =min(un_1 , CFraCM(ikl,k))
	1531	CFraCM(ikl,k) =max(0.001 , CFraCM(ikl,k)) &!
	1532	& * max(zer0,sign(un_1,qi__CM(ikl,k) + qw__CM(ikl,k) &!
	1533	& +qs__CM(ikl,k) -3.E-9 ))!
	1534	END DO
	1535	ELSE ! XU and Randall 1996, JAS 21, p.3099 (4)
	1536	! ----------------------------
	1537	DO k=mz1_CM,mzp
	1538	qvs_wi= qvswCM(ikl,k)
	1539	! #wi qvs_wi=max(epsn,((qi__CM(ikl,k)+qs__CM(ikl,k))*qvsiCM(ikl,k) &
	1540	! #wi& +qw__CM(ikl,k) *qvswCM(ikl,k)) &
	1541	! #wi& /max(epsn, qi__CM(ikl,k)+qs__CM(ikl,k) +qw__CM(ikl,k)))
	1542	RHumid= min( RH_MAX, max(qv__DY(ikl,k) ,qv_MIN) &
	1543	& / qvs_wi)
	1544	argEXP= ( (RH_MAX -RHumid) * qvs_wi) ** 0.49
	1545	argEXP= min(100.*(qi__CM(ikl,k)+qw__CM(ikl,k) &
	1546	& +qs__CM(ikl,k) * 0.33 &
	1547	& * (1.-min(1.,exp((Ta__CM(ikl,k) -258.15)*0.1))))&
	1548	& /max( epsn , argEXP ) &
	1549	& ,ea_MAX )
	1550
	1551	CFraCM(ikl,k) = ( RHumid ** 0.25 )&
	1552	& * (1. - exp(-argEXP) )
	1553	END DO
	1554	END IF
	1555
	1556	ELSE
	1557	! #sc ELSE IF ( .NOT.Frac__Clouds) THEN
	1558	! #sc IF (fracSC) stop 'fracSC set up when Frac__Clouds NOT'
	1559	DO k=mz1_CM,mzp
	1560	qCloud = qi__CM(ikl,k) + qw__CM(ikl,k)
	1561
	1562	! #hy IF (qCloud &gt.epsn) THEN
	1563	CFraCM(ikl,k) = max(zer0,sign(un_1, qCloud - epsn))
	1564	! CFraCM(ikl,k) = 1.0 if qCloud > epsn
	1565	! = 0.0 otherwise
	1566
	1567	! #hy END IF
	1568	END DO
	1569
	1570	END IF
	1571
	1572
	1573	! Debug
	1574	! ~~~~~
	1575	! #wH DO k=mz1_CM,mzp
	1576	! #wH debugH( 1:35) = 'Fractional Cloudiness (XU .OR. CEP)'
	1577	! #wH debugH(36:70) = ' '
	1578	! #wH proc_1 = ' '
	1579	! #wH procv1 = 0.
	1580	! #wH proc_2 = ' '
	1581	! #wH procv2 = 0.
	1582	! #wH proc_3 = ' '
	1583	! #wH procv3 = 0.
	1584	! #wH proc_4 = ' '
	1585	! #wH procv4 = 0.
	1586	! #wh include 'CMiPhy_Debug.h'
	1587	! #wH END DO
	1588
	1589
	1590
	1591
	1592	!=============================================================================== AUTO-CONVERSION, LIQUID
	1593	! +++++++++++++++++++++++
	1594	! Autoconversion (i.e., generation of precipitating particles), liquid water
	1595	! ==========================================================================
	1596
	1597	! Cloud Droplets Autoconversion
	1598	! Reference: Sundqvist 1988, Schlesinger, Reidel, p. 433)
	1599	! Reference: Lin et al. 1983, JCAM 22, p.1076 (50)
	1600	! ----------------------------------------------------------
	1601
	1602	DO k=mz1_CM,mzp
	1603
	1604	! #wH qr_AUT = 0.0
	1605
	1606	! #hy IF (qw__CM(ikl,k).gt.epsn) THEN
	1607	qw__OK = max(zer0,sign(un_1,qw__CM(ikl,k) - epsn))
	1608	! qw__OK = 1.0 if qw__CM(ikl,k) > epsn
	1609	! = 0.0 otherwise
	1610
	1611	! #hy IF (CFraCM(ikl,k).gt.CFrMIN) THEN
	1612	CFraOK = max(zer0,sign(un_1,CFraCM(ikl,k) - CFrMIN))
	1613	! CFraOK = 1.0 if CFraCM(ikl,k) > CFrMIN
	1614	! = 0.0 otherwise
	1615
	1616	qw__OK = qw__OK * CFraOK
	1617
	1618	! #EW IF(qw__OK.gt.eps6) THEN
	1619	! #EW mauxEW = mphyEW(ikl)
	1620	! #EW mauxEW(08:08) = 'r'
	1621	! #EW mphyEW(ikl) = mauxEW
	1622	! #EW END IF
	1623
	1624	! Sundqvist (1988, Schlesinger, Reidel, p. 433) Autoconversion Scheme
	1625	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1626	IF (AUTO_w_Sundqv) THEN
	1627	dwMesh = qw__OK *qw__CM(ikl,k)/qw_MAX &!
	1628	& /max(CFrMIN,CFraCM(ikl,k)) !
	1629	pr_AUT = qw__OK qw__CM(ikl,k)c_Sund &!
	1630	& (1.-exp(-min(dwMeshdwMesh,ea_MAX))) &!
	1631	& /max(CFrMIN,CFraCM(ikl,k)) !
	1632
	1633	! Liou and Ou (1989, JGR 94, p. 8599) Autoconversion Scheme
	1634	! Boucher et al. (1995, JGR 100, p.16395) ~~~~~~~~~~~~~~~~~~~~~
	1635	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1636	ELSE IF (AUTO_w_LiouOu) THEN
	1637	CCNwCM(ikl,k) = 1.2e+8 ! ASTEX (Duynkerke&al.1995, JAS 52, p.2763)
	1638	! #lo CCNwCM(ikl,k) = 1.e+11 ! (polluted air, Rogers&Yau 89, p.90)
	1639
	1640	qwCFra = qw__CM(ikl,k) / CFraCM(ikl,k)
	1641	dwTUR4 = 4.5 * qwTURB * qwTURB
	1642	dwTURi = qwCFra * roa_DY(ikl,k) &!
	1643	& * 6.0 /piNmbr / CCNwCM(ikl,k) /exp(dwTUR4)
	1644	dwTUR3 = exp(R_1by3*log(dwTURi))
	1645	dwTUR2 = dwTUR3 * dwTUR3
	1646	dwTUR8 = 8.0 * qwTURB * qwTURB
	1647	dwTURc = exp(dwTUR8) * dwTUR2 * dwTUR2
	1648	rwMEAN = 0.5 *sqrt(sqrt(dwTURc))
	1649
	1650	th_AUT = max(zer0,sign(un_1, rwMEAN -rwCrit)) ! Heaviside Function
	1651
	1652	pr_AUT = qw__OKCFraCM(ikl,k) th_AUT4.09d6piNmbr &!
	1653	& CCNwCM(ikl,k) dwTURc*qwCFra
	1654
	1655	! Lin et al.(1983) Autoconversion Scheme
	1656	! ~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~
	1657	ELSE IF (AUTO_w_LinAll) THEN
	1658	dwMesh = qw__OK * (qw__CM(ikl,k)-qw_MAX)
	1659	pr_AUT = dwMesh * dwMesh dwMesh/(cc1dwMesh+1000.d0*cc2/dd0)
	1660	ELSE
	1661	STOP 'AutoConversion of Cloud droplets is not defined'
	1662	END IF
	1663
	1664	qr_AUT = pr_AUT * dt__CM
	1665	qr_AUT = min(qr_AUT,qw__CM(ikl,k))
	1666	qw__CM(ikl,k) = qw__CM(ikl,k) - qr_AUT
	1667	qr__CM(ikl,k) = qr__CM(ikl,k) + qr_AUT
	1668
	1669	! #WQ write(6,*) 'QrAut',qr_AUT,it_EXP,ikl,k
	1670	! #WH if (ikl.eq.ikl0CM(1)) wraut(k) = qr_AUT
	1671
	1672	! #hy END IF
	1673	! #hy END IF
	1674
	1675	! Debug
	1676	! ~~~~~
	1677	! #wH debugH( 1:35) = 'Lin et al.(1983) Autoconversion Sch'
	1678	! #wH debugH(36:70) = 'eme '
	1679	! #wH proc_1 = 'Qraut g/kg'
	1680	! #wH procv1 = qr_AUT
	1681	! #wH proc_2 = ' '
	1682	! #wH procv2 = 0.
	1683	! #wH proc_3 = ' '
	1684	! #wH procv3 = 0.
	1685	! #wH proc_4 = ' '
	1686	! #wH procv4 = 0.
	1687	! #wh include 'CMiPhy_Debug.h'
	1688	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1689	! #wH& debugV(k,06) = qr_AUT
	1690
	1691	END DO
	1692
	1693
	1694
	1695
	1696	!=============================================================================== AUTO-CONVERSION, SOLID
	1697	! ++++++++++++++++++++++
	1698	! Autoconversion (i.e., generation of precipitating particles), Ice --> Snow
	1699	! ==========================================================================
	1700
	1701	! Conversion from Cloud Ice Crystals to Snow Flakes
	1702	! Reference: Levkov et al. 1992, Contr.Atm.Phys. 65, p.41
	1703	! ---------------------------------------------------------
	1704
	1705	IF (AUTO_i_Levkov) THEN
	1706
	1707
	1708	! Depositional Growth: Ice Crystals => Snow Flakes (BDEPIS)
	1709	! Reference: Levkov et al. 1992, Contr.Atm.Phys. 65, p.41 (28)
	1710	! --------------------------------------------------------------
	1711
	1712	IF (AUTO_i_LevkXX) THEN
	1713
	1714	DO k=mz1_CM,mzp
	1715
	1716	! #wH qs_AUT = 0.0
	1717
	1718	! #hy IF (qi__CM(ikl,k).gt.epsn) THEN
	1719	qi__OK = max(zer0, sign(un_1,qi__CM(ikl,k) - epsn))
	1720	! qi__OK = 1.0 if qi__CM(ikl,k) > epsn
	1721	! = 0.0 otherwise
	1722
	1723	! #hy IF (CCNiCM(ikl,k).gt.1.e0) THEN
	1724	CCNiOK = max(zer0, sign(un_1,CCNiCM(ikl,k) - 1.e0))
	1725	! CCNiOK = 1.0 if CCNiCM(ikl,k) > 1.e0
	1726	! = 0.0 otherwise
	1727
	1728	qi__OK = qi__OK * CCNiOK * qi__CM(ikl,k)
	1729
	1730	! Pristine Ice Crystals Diameter
	1731	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1732	Di_Pri = 0.156 exp(R_1by3log(R_1000*roa_DY(ikl,k) &! Pristine Ice Crystals Diameter
	1733	& *max(epsn ,qi__CM(ikl,k)) &! Levkov et al. 1992, Contr.Atm.Phys. 65, (5) p.37
	1734	& /max(un_1 ,CCNiCM(ikl,k)))) ! where 6/(piro_I)*1/3 ~ 0.156
	1735
	1736	! Deposition Time Scale
	1737	! ~~~~~~~~~~~~~~~~~~~~~
	1738	RH_Ice = max(epsq, qv__DY(ikl,k)) / qsiEFF(ikl,k)
	1739
	1740	dtsaut = 0.125 (qs__D0qs__D0-Di_Pri*Di_Pri) &!
	1741	& (0.702e12/(Ta__CM(ikl,k)Ta__CM(ikl,k)) &! 0.702e12 ~ 0.701987755e12 = (2.8345e+6)**2/0.0248/461.5
	1742	! Ls_H2O **2/Ka /Rw
	1743	& +1.0 /(2.36e-2 roa_DY(ikl,k) &! 2.36e-2 = 2.36e-5 1.e3
	1744	& max(epsq,qv__DY(ikl,k))RH_Ice)) ! Dv
	1745
	1746	! Deposition
	1747	! ~~~~~~~~~~
	1748	qs_AUT = dt__CM qi__OK(RH_Ice-1.)/dtsaut
	1749	qs_AUT = min( qi__CM(ikl,k) , qs_AUT)
	1750	qs_AUT = max(-qs__CM(ikl,k) , qs_AUT)
	1751	qi__CM(ikl,k) = qi__CM(ikl,k) - qs_AUT
	1752	qs__CM(ikl,k) = qs__CM(ikl,k) + qs_AUT
	1753
	1754	! #hy END IF
	1755	! #hy END IF
	1756
	1757	! Debug
	1758	! ~~~~~
	1759	! #wH debugH( 1:35) = 'Lin et al.(1983) Depositional Growt'
	1760	! #wH debugH(36:70) = 'h '
	1761	! #wH proc_1 = 'QsAUT g/kg'
	1762	! #wH procv1 = qs_AUT
	1763	! #wH proc_2 = ' '
	1764	! #wH procv2 = 0.
	1765	! #wH proc_3 = ' '
	1766	! #wH procv3 = 0.
	1767	! #wH proc_4 = ' '
	1768	! #wH procv4 = 0.
	1769	! #wh include 'CMiPhy_Debug.h'
	1770	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1771	! #wH& debugV(k,07) = qs_AUT
	1772
	1773	END DO
	1774
	1775	END IF
	1776
	1777
	1778	! Ice Crystals Aggregation => Snow Flakes (BAGRIS) AUTO-CONVERSION, SOLID
	1779	! Reference: Levkov et al. 1992, Contr.Atm.Phys. 65, p.41 (31) ++++++++++++++++++++++
	1780	! --------------------------------------------------------------
	1781
	1782	DO k=mz1_CM,mzp
	1783
	1784	! #wH qs_AUT = 0.0
	1785	! #wH dtsaut = 0.0
	1786
	1787	! #hy IF (qi__CM(ikl,k).gt.epsn) THEN
	1788	qi__OK = max(zer0,sign(un_1,qi__CM(ikl,k) - epsn))
	1789	! qi__OK = 1.0 if qi__CM(ikl,k) > epsn
	1790	! = 0.0 otherwise
	1791
	1792	! #hy IF (CCNiCM(ikl,k).gt.1.e0) THEN
	1793	CCNiOK = max(zer0,sign(un_1,CCNiCM(ikl,k) - 1.e0))
	1794	! CCNiOK = 1.0 if CCNiCM(ikl,k) > 1.e0
	1795	! = 0.0 otherwise
	1796
	1797	qi__OK = qi__OK * CCNiOK * qi__CM(ikl,k)
	1798
	1799	! #EW IF(qi__OK.gt.eps6) THEN
	1800	! #EW mauxEW = mphyEW(ikl)
	1801	! #EW mauxEW(09:09) = 's'
	1802	! #EW mphyEW(ikl) = mauxEW
	1803	! #EW END IF
	1804
	1805	! Pristine Ice Crystals Diameter
	1806	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1807	Di_Pri = 0.156 exp(R_1by3log(R_1000*roa_DY(ikl,k) &! Pristine Ice Crystals Diameter
	1808	& *max(epsn, qi__CM(ikl,k)) &! Levkov et al. 1992, Contr. Atm. Phys. 65, (5) p.37
	1809	& /max(un_1, CCNiCM(ikl,k)))) ! where [6/(piro_I)]*1/3 ~ 0.156
	1810
	1811	! Time needed for Ice Crystals Diameter to reach Snow Diameter Threshold
	1812	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1813	c1saut = max(epsn,qi__OK) roa_DY(ikl,k) 35.0 &!
	1814	& exp(R_1by3log(roa_DY(ikl,mzp) /roa_DY(ikl,k))) !
	1815
	1816	dtsaut =-6.d0*log(Di_Pri/qs__D0) /c1saut !
	1817	dtsaut = max(dt__CM, dtsaut) ! qi fully used if dtsaut<dt__CM
	1818
	1819	! Time needed for Ice Crystals Diameter to reach Snow Diameter Threshold
	1820	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(ALTERNATE PARAMETERIZATION)~
	1821	! #nt dtsaut =-2.0 (3.0log( Di_Pri /qs__D0) &!
	1822	! #nt& + log(max(qi__CM(ikl,k),epsn))) /c1saut !
	1823	! #nt dtsaut = max(epsn,dtsaut)
	1824
	1825	! Aggregation
	1826	! ~~~~~~~~~~~
	1827	qs_AUT = dt__CM*qi__OK / dtsaut
	1828	qs_AUT = min( qi__CM(ikl,k) , qs_AUT)
	1829	qs_AUT = max(-qs__CM(ikl,k) , qs_AUT)
	1830	qi__CM(ikl,k) = qi__CM(ikl,k) - qs_AUT
	1831	qs__CM(ikl,k) = qs__CM(ikl,k) + qs_AUT
	1832
	1833
	1834	! Decrease of Ice Crystals Number (BAGRII)
	1835	! Reference: Levkov et al. 1992, Contr.Atm.Phys. 65, p.41 (34)
	1836	! --------------------------------------------------------------
	1837
	1838	CCNiCM(ikl,k) = CCNiCM(ikl,k) * exp(-0.5c1sautdt__CM)
	1839
	1840	! #WQ write(6,*) 'QsAut', qs_AUT, &!
	1841	! #WQ& ' Qi' , qi__CM(ikl,k), &!
	1842	! #WQ& ' CcnI' ,CCNiCM(ikl,k),it_EXP,ikl,k !
	1843	! #WH if (ikl.eq.ikl0CM(1)) wsaut(k) = qs_AUT
	1844
	1845	! #hy END IF
	1846	! #hy END IF
	1847
	1848	! Debug
	1849	! ~~~~~
	1850	! #wH debugH( 1:35) = 'Lin et al.(1983) Ice Crystals Aggre'
	1851	! #wH debugH(36:70) = 'gation '
	1852	! #wH proc_1 = 'dtsaut sec'
	1853	! #wH procv1 = dtsaut
	1854	! #wH proc_2 = 'QsAUT g/kg'
	1855	! #wH procv2 = qs_AUT
	1856	! #wH proc_3 = ' '
	1857	! #wH procv3 = 0.
	1858	! #wH proc_4 = ' '
	1859	! #wH procv4 = 0.
	1860	! #wh include 'CMiPhy_Debug.h'
	1861	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1862	! #wH& debugV(k,07) = qs_AUT + debugV(k,07)
	1863
	1864	END DO
	1865
	1866
	1867	! Ice Crystals Autoconversion => Snow Flakes AUTO-CONVERSION, SOLID
	1868	! Reference: Lin et al. 1983, JCAM 22, p.1070 (21) ++++++++++++++++++++++
	1869	! Lin et al. 1983, JCAM 22, p.1074 (38)
	1870	! Emde and Kahlig 1989, Ann.Geoph. 7, p. 408 (18)
	1871	! ----------------------------------------------------------
	1872
	1873	ELSE IF (AUTO_i_EmdeKa) THEN
	1874
	1875	DO k=mz1_CM,mzp
	1876
	1877	! #wH qs_AUT = 0.0
	1878	! #wH cnsaut = 0.0
	1879
	1880	! #hy IF (qi__CM(ikl,k) .ge. qisMAX) THEN
	1881
	1882	ps_AUT = 0.001d0*(qi__CM(ikl,k)-qisMAX) &!
	1883	& exp(0.025d0 Ta_dgC(ikl,k))
	1884	qs_AUT = ps_AUT * dt__CM
	1885	qs_AUT = max(qs_AUT, zer0 )
	1886	qs_AUT = min(qs_AUT, qi__CM(ikl,k))
	1887	cnsaut = qs_AUT* CCNiCM(ikl,k) &!
	1888	& /max(qisMAX , qi__CM(ikl,k))
	1889	CCNiCM(ikl,k) = CCNiCM(ikl,k) - cnsaut
	1890	qi__CM(ikl,k) = qi__CM(ikl,k) - qs_AUT
	1891	qs__CM(ikl,k) = qs__CM(ikl,k) + qs_AUT
	1892	! #WQ write(6,*) 'QsAut',qs_AUT ,it_EXP,ikl,k
	1893	! #WH IF (ikl.eq.ikl0CM(1)) wsaut(k)= qs_AUT
	1894
	1895	! #EW IF (qi__CM(ikl,k) .ge. qisMAX) THEN
	1896	! #EW mauxEW = mphyEW(ikl)
	1897	! #EW mauxEW(09:09) = 's'
	1898	! #EW mphyEW(ikl) = mauxEW
	1899	! #EW END IF
	1900
	1901	! #hy END IF
	1902
	1903	! Debug
	1904	! ~~~~~
	1905	! #wH debugH( 1:35) = 'Emde and Kahlig Ice Crystals Autoc'
	1906	! #wH debugH(36:70) = 'onversion '
	1907	! #wH proc_1 = 'QsAUT g/kg'
	1908	! #wH procv1 = qs_AUT
	1909	! #wH proc_2 = 'cnsaut/e15'
	1910	! #wH procv2 = cnsaut*1.e-18
	1911	! #wH proc_3 = ' '
	1912	! #wH procv3 = 0.
	1913	! #wH proc_4 = ' '
	1914	! #wH procv4 = 0.
	1915	! #wh include 'CMiPhy_Debug.h'
	1916	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1917	! #wH& debugV(k,07) = qs_AUT
	1918
	1919	END DO
	1920
	1921
	1922	! Sundqvist (1988, Schlesinger, Reidel, p. 433) Autoconversion Scheme AUTO-CONVERSION, SOLID
	1923	! ------------------------------------------------------------------------- ++++++++++++++++++++++
	1924
	1925	ELSE IF (AUTO_i_Sundqv) THEN
	1926
	1927	DO k=mz1_CM,mzp
	1928
	1929	! #wH qs_AUT = 0.0
	1930	! #wH cnsaut = 0.0
	1931
	1932	! #hy IF (qi__CM(ikl,k).gt.epsn) THEN
	1933	qi__OK = max(zer0,sign(un_1,qi__CM(ikl,k) - epsn))
	1934	! qi__OK = 1.0 if qi__CM(ikl,k) > epsn
	1935	! = 0.0 otherwise
	1936
	1937	dqiDUM = qi__OK *qi__CM(ikl,k)/qi0_DC &!
	1938	! #mf& /max( CFrMIN ,CFraCM(ikl,k)) &!
	1939	& + 0. !
	1940	ps_AUT = qi__OK qi__CM(ikl,k)c_Sund &!
	1941	& (1.-exp(-dqiDUM dqiDUM)) &!
	1942	! #mf& *max( CFrMIN ,CFraCM(ikl,k)) &!
	1943	& + 0. !
	1944	qs_AUT = ps_AUT * dt__CM
	1945	qs_AUT = min(qi__CM(ikl,k) , qs_AUT)
	1946	qs_AUT = max(zer0 , qs_AUT)
	1947	cnsaut = CCNiCM(ikl,k) * qs_AUT &!
	1948	& /max(qi__CM(ikl,k) , epsn)
	1949	CCNiCM(ikl,k) = CCNiCM(ikl,k) - cnsaut
	1950	qi__CM(ikl,k) = qi__CM(ikl,k) - qs_AUT
	1951	qs__CM(ikl,k) = qs__CM(ikl,k) + qs_AUT
	1952
	1953	! #hy END IF
	1954
	1955	! Debug
	1956	! ~~~~~
	1957	! #wH debugH( 1:35) = 'Sundqvist (1988) Ice Crystals Autoc'
	1958	! #wH debugH(36:70) = 'onversion '
	1959	! #wH proc_1 = 'QsAUT g/kg'
	1960	! #wH procv1 = qs_AUT
	1961	! #wH proc_2 = 'cnsaut/e15'
	1962	! #wH procv2 = cnsaut*1.e-18
	1963	! #wH proc_3 = ' '
	1964	! #wH procv3 = 0.
	1965	! #wH proc_4 = ' '
	1966	! #wH procv4 = 0.
	1967	! #wh include 'CMiPhy_Debug.h'
	1968	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	1969	! #wH& debugV(k,07) = qs_AUT
	1970
	1971	END DO
	1972	ELSE
	1973	STOP 'AutoConversion of Cloud crystals is not defined'
	1974	END IF
	1975
	1976
	1977
	1978
	1979	!=============================================================================== AUTO-CONVERSION, SOLID
	1980	! ++++++++++++++++++++++
	1981	! Autoconversion (i.e., generation of precipitating particles), Ice --> Graupels
	1982	! ==============================================================================
	1983
	1984	! #qg DO k=mz1_CM,mzp
	1985
	1986	! #qg IF (qi__CM(ikl,k) .ge. qigMAX) THEN
	1987
	1988	! #qg pgaut = 0.001( qi__CM(ikl,k)-qigMAX)exp(0.090* Ta_dgC(ikl,k))
	1989	! #qg qgaut = pgaut * dt__CM
	1990	! #qg qgaut = max(qgaut,zer0 )
	1991	! #qg qgaut = min(qgaut,qi__CM(ikl,k))
	1992	! #qg qi__CM(ikl,k) = qi__CM(ikl,k) - qgaut
	1993	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qgaut
	1994
	1995	! #qg END IF
	1996
	1997	! #qg END DO
	1998
	1999
	2000
	2001
	2002	!=============================================================================== ACCRETION
	2003	! +++++++++
	2004	! Accretion Processes (i.e. increase in size of precipitating particles
	2005	! ==================== through a collision-coalescence process)===
	2006	! ==============================================
	2007
	2008	! Accretion of Cloud Droplets by Rain ACCRETION, o > .
	2009	! Reference: Lin et al. 1983, JCAM 22, p.1076 (51) +++++++++++++++++
	2010	! Emde and Kahlig 1989, Ann.Geoph. 7, p. 407 (10)
	2011	! ----------------------------------------------------------
	2012
	2013	DO k=mz1_CM,mzp
	2014
	2015	! #wH qr_ACW = 0.0
	2016
	2017	! #hy IF (qw__CM(ikl,k).gt.epsn) THEN
	2018	qw__OK = max(zer0,sign(un_1,qw__CM(ikl,k) - epsn))
	2019	! qw__OK = 1.0 if qw__CM(ikl,k) > epsn
	2020	! = 0.0 otherwise
	2021
	2022	! #hy IF (qr___0(ikl,k).gt.epsn) THEN
	2023	qr0_OK = max(zer0,sign(un_1,qr___0(ikl,k) - epsn))
	2024	! qr0_OK = 1.0 if qr___0(ikl,k) > epsn
	2025	! = 0.0 otherwise
	2026
	2027	Flag_qr_ACW = qw__OK * qr0_OK
	2028
	2029	! #EW IF(Flag_qr_ACW.gt.eps6) THEN
	2030	! #EW mauxEW = mphyEW(ikl)
	2031	! #EW mauxEW(10:10) = 'r'
	2032	! #EW mphyEW(ikl) = mauxEW
	2033	! #EW END IF
	2034
	2035	pr_ACW = 3104.28d0 * n0___r * sqrrro(ikl,k) &! 3104.28 = a pi Gamma[3+b] / 4
	2036	& qw__CM(ikl,k)/exp(3.8d0log(lamdaR(ikl,k))) ! where a = 842. and b = 0.8
	2037	qr_ACW = pr_ACWdt__CM Flag_qr_ACW
	2038	qr_ACW = min(qr_ACW,qw__CM(ikl,k))
	2039
	2040	qw__CM(ikl,k) = qw__CM(ikl,k) -qr_ACW
	2041	qr__CM(ikl,k) = qr__CM(ikl,k) +qr_ACW
	2042
	2043	! #WQ write(6,*) 'Qracw',qr_ACW,it_EXP,ikl,k
	2044	! #WH if (ikl.eq.ikl0CM(1)) wracw(k) = qr_ACW
	2045
	2046	! #hy END IF
	2047	! #hy END IF
	2048
	2049	! Debug
	2050	! ~~~~~
	2051	! #wH debugH( 1:35) = 'Lin et al.(1983): Accretion of Clou'
	2052	! #wH debugH(36:70) = 'd Droplets by Rain '
	2053	! #wH proc_1 = 'Qracw g/kg'
	2054	! #wH procv1 = qr_ACW
	2055	! #wH proc_2 = ' '
	2056	! #wH procv2 = 0.
	2057	! #wH proc_3 = ' '
	2058	! #wH procv3 = 0.
	2059	! #wH proc_4 = ' '
	2060	! #wH procv4 = 0.
	2061	! #wh include 'CMiPhy_Debug.h'
	2062	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2063	! #wH& debugV(k,08) = qr_ACW
	2064
	2065	END DO
	2066
	2067
	2068	! Accretion of Cloud Droplets by Snow Flakes ACCRETION, * > .
	2069	! Reference: Lin et al. 1983, JCAM 22, p.1070 (24) ++++++++++++++++
	2070	! ----------------------------------------------------------
	2071
	2072	DO k=mz1_CM,mzp
	2073
	2074	! #hy IF (qw__CM(ikl,k).gt.epsn) THEN
	2075	qw__OK = max(zer0,sign(un_1,qw__CM(ikl,k) - epsn))
	2076	! qw__OK = 1.0 if qw__CM(ikl,k) > epsn
	2077	! = 0.0 otherwise
	2078
	2079	! #hy IF (qs___0(ikl,k).gt.epsn) THEN
	2080	qs0_OK = max(zer0,sign(un_1,qs___0(ikl,k) - epsn))
	2081	! qs0_OK = 1.0 if qs___0(ikl,k) > epsn
	2082	! = 0.0 otherwise
	2083
	2084	Flag_qs_ACW = qw__OK * qs0_OK
	2085
	2086	! #EW IF(Flag_qs_ACW.gt.eps6) THEN
	2087	! #EW mauxEW = mphyEW(ikl)
	2088	! #EW mauxEW(11:11) = 's'
	2089	! #EW mphyEW(ikl) = mauxEW
	2090	! #EW END IF
	2091
	2092	! #cn n0___s = min(2.e8,2.e6exp(-.12min(0.,Ta_dgC(ikl,k))))
	2093
	2094	! ps_ACW is taken into account in the snow melting process (if positive temperatures)
	2095	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2096	IF (graupel_shape) THEN! Graupellike Snow Flakes of Hexagonal Type
	2097	ps_ACW(ikl,k)= 9.682d0 * n0___s * sqrrro(ikl,k) &! 9.682 = c pi Gamma[3+d] / 4
	2098	& qw__CM(ikl,k) /exp(3.25d0log(lamdaS(ikl,k))) ! where c = 4.836 and d = 0.25
	2099	! Ref.: Locatelli and Hobbs, 1974, JGR: table 1 p.2188
	2100
	2101	ELSE IF (planes__shape) THEN! Unrimed Side Plane
	2102	ps_ACW(ikl,k)= 3517. * n0___s * sqrrro(ikl,k) &! 3517. = c pi Gamma[3+d] / 4
	2103	& qw__CM(ikl,k) /exp(3.99d0log(lamdaS(ikl,k))) ! where c = 755.9 and d = 0.99
	2104
	2105	ELSE IF (aggrega_shape) THEN! Aggregates of unrimed radiating assemblages
	2106	ps_ACW(ikl,k)= 27.73 * n0___s * sqrrro(ikl,k) &! 27.73 = c pi Gamma[3+d] / 4
	2107	& qw__CM(ikl,k) /exp(3.41d0log(lamdaS(ikl,k))) ! where c = 11.718and d = 0.41
	2108
	2109	ELSE
	2110	STOP 'Snow Particles Shape is not defined'
	2111	END IF
	2112
	2113	qs_ACW = dt__CMps_ACW(ikl,k)Flag_qs_ACW
	2114	qs_ACW = min(qs_ACW,qw__CM(ikl,k))
	2115
	2116	Flag_Ta_Pos = max(zer0,sign(un_1,Ta__CM(ikl,k) - Tf_Sno))
	2117	! Flag_Ta_Pos = 1.0 if Ta__CM(ikl,k) > Tf_Sno
	2118	! = 0.0 otherwise
	2119
	2120	qw__CM(ikl,k) = qw__CM(ikl,k) - qs_ACW
	2121	qr__CM(ikl,k) = qr__CM(ikl,k) + Flag_Ta_Pos * qs_ACW
	2122	Flag_qs_ACW = (1.d0 - Flag_Ta_Pos) * qs_ACW
	2123	qs__CM(ikl,k) = qs__CM(ikl,k) + Flag_qs_ACW
	2124	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Lc_Cpd * Flag_qs_ACW
	2125	! Negative Temperatures => Latent Heat is released by Freezing
	2126
	2127	! Full Debug
	2128	! ~~~~~~~~~~
	2129	! #WQ write(6,*) 'Qsacw',qs_ACW,it_EXP,ikl,k
	2130	! #WH if (ikl.eq.ikl0CM(1)) wsacw(k) = qs_ACW
	2131
	2132	! #hy END IF
	2133	! #hy END IF
	2134
	2135	! Debug
	2136	! ~~~~~
	2137	! #wH debugH( 1:35) = 'Lin et al.(1983): Accretion of Clou'
	2138	! #wH debugH(36:70) = 'd Droplets by Snow Particles '
	2139	! #wH proc_1 = 'Qsacw g/kg'
	2140	! #wH procv1 = Flag_qs_ACW
	2141	! #wH proc_3 = ' '
	2142	! #wH procv2 = 0.
	2143	! #wH proc_2 = ' '
	2144	! #wH procv3 = 0.
	2145	! #wH proc_4 = ' '
	2146	! #wH procv4 = 0.
	2147	! #wh include 'CMiPhy_Debug.h'
	2148	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2149	! #wH& debugV(k,09) = Flag_qs_ACW
	2150
	2151	END DO
	2152
	2153
	2154	! Accretion of Cloud Droplets by Graupels (Dry Growth Mode) ACCRETION, # > . \| #
	2155	! Reference: Lin et al. 1983, JCAM 22, p.1075 (40) ++++++++++++++++++++
	2156	! Emde and Kahlig 1989, Ann.Geoph. 7, p. 407 (~20)
	2157	! -----------------------------------------------------------
	2158
	2159	! #qg DO k=mz1_CM,mzp
	2160
	2161	! #qg IF (qw__CM(ikl,k).gt.epsn) THEN
	2162	! #qg WbyG_w = max(zer0, sign(un_1,qw__CM(ikl,k) - epsn))
	2163	! WbyG_w = 1.0 if qw__CM(ikl,k) > epsn
	2164	! = 0.0 otherwise
	2165
	2166	! #qg IF (qg__CM(ikl,k).gt.epsn) THEN
	2167	! #qg WbyG_g = max(zer0, sign(un_1,qg__CM(ikl,k) - epsn))
	2168	! WbyG_g = 1.0 if qg__CM(ikl,k) > epsn
	2169	! = 0.0 otherwise
	2170
	2171	! #qg WbyGOK = WbyG_w * WbyG_g
	2172
	2173	! #qg IF (Ta__CM(ikl,k).lt.Tf_Sno) THEN
	2174	! #qg Fact_G = max(zer0,-sign(un_1,Ta__CM(ikl,k) - Tf_Sno))
	2175	! Fact_G = 1.0 if Ta__CM(ikl,k) > Tf_Sno
	2176	! = 0.0 otherwise
	2177
	2178	! #qg pgacw = PATATRAS
	2179	! #qg qgacw = pgacw * dt__CM * WbyGOK
	2180	! #qg qgacw = min(qgacw,qw__CM(ikl,k))
	2181
	2182	! #qg qw__CM(ikl,k) = qw__CM(ikl,k) - qgacw
	2183	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qgacw
	2184	! #qg Ta__CM(ikl,k) = Ta__CM(ikl,k) +Lc_Cpd gacw
	2185
	2186	! #qg END IF
	2187	! #qg END IF
	2188	! #qg END IF
	2189
	2190	! #qg END DO
	2191
	2192
	2193
	2194
	2195	! Accretion of Cloud Ice by Snow Particles ACCRETION, * > /
	2196	! Reference: Lin et al. 1983, JCAM 22, p.1070 (22) ++++++++++++++++
	2197	! ----------------------------------------------------------
	2198
	2199	DO k=mz1_CM,mzp
	2200
	2201	! #wH qs_ACI = 0.0
	2202	! #wH CNsACI = 0.0
	2203
	2204	! #hy IF (qi__CM(ikl,k).gt.epsn) THEN
	2205	qi__OK = max(zer0,sign(un_1,qi__CM(ikl,k) - epsn))
	2206	! qi__OK = 1.0 if qi__CM(ikl,k) > epsn
	2207	! = 0.0 otherwise
	2208
	2209	! #hy IF (qs___0(ikl,k).gt.epsn) THEN
	2210	qs0_OK = max(zer0,sign(un_1,qs___0(ikl,k) - epsn))
	2211	! qs0_OK = 1.0 if qs___0(ikl,k) > epsn
	2212	! = 0.0 otherwise
	2213
	2214	! #hy IF (Ta__CM(ikl,k).lt.Tf_Sno) THEN
	2215	Flag_Ta_Neg = max(zer0,-sign(un_1,Ta__CM(ikl,k) - Tf_Sno))
	2216	! Flag_Ta_Neg = 1.0 if Ta__CM(ikl,k) < Tf_Sno
	2217	! = 0.0 otherwise
	2218
	2219	Flag_qs_ACI = qi__OK * qs0_OK * Flag_Ta_Neg
	2220
	2221	! #EW IF(Flag_qs_ACI.gt.eps6) THEN
	2222	! #EW mauxEW = mphyEW(ikl)
	2223	! #EW mauxEW(12:12) = 's'
	2224	! #EW mphyEW(ikl) = mauxEW
	2225	! #EW END IF
	2226
	2227	effACI = exp(0.025d0*Ta_dgC(ikl,k)) ! Collection Efficiency
	2228	! Lin et al. 1983 JCAM 22 p.1070 (23)
	2229
	2230	! #cn n0___s = min(2.e8,2.e6exp(-.12min(0.,Ta_dgC(ikl,k))))
	2231
	2232	! ps_ACI
	2233	! ~~~~~~
	2234	IF (graupel_shape) THEN
	2235	ps_ACI = effACI * 9.682d0 * n0___s * sqrrro(ikl,k) &
	2236	& qi__CM(ikl,k) /exp(3.25d0log(lamdaS(ikl,k)))
	2237
	2238	ELSE IF (planes__shape) THEN
	2239	ps_ACI = effACI * 3517.d0 * n0___s * sqrrro(ikl,k) &
	2240	& qi__CM(ikl,k) /exp(3.99d0log(lamdaS(ikl,k)))
	2241
	2242	ELSE IF (aggrega_shape) THEN
	2243	ps_ACI = effACI * 27.73d0 * n0___s * sqrrro(ikl,k) &
	2244	& qi__CM(ikl,k) /exp(3.41d0log(lamdaS(ikl,k)))
	2245	ELSE
	2246	STOP 'Snow Particles Shape is not defined'
	2247	END IF
	2248
	2249	qs_ACI = ps_ACI * dt__CM * Flag_qs_ACI
	2250	qs_ACI = min(qs_ACI,qi__CM(ikl,k))
	2251
	2252	CNsACI = CCNiCM(ikl,k) * qs_ACI &!
	2253	& /max(qi__CM(ikl,k) , epsn)
	2254	CCNiCM(ikl,k) = CCNiCM(ikl,k) - CNsACI
	2255	qi__CM(ikl,k) = qi__CM(ikl,k) - qs_ACI
	2256	qs__CM(ikl,k) = qs__CM(ikl,k) + qs_ACI
	2257
	2258	! #WQ write(6,*) 'Qsaci',qs_ACI,it_EXP,ikl,k
	2259	! #WH if (ikl.eq.ikl0CM(1)) wsaci(k) = qs_ACI
	2260
	2261	! #hy END IF
	2262	! #hy END IF
	2263	! #hy END IF
	2264
	2265	! Debug
	2266	! ~~~~~
	2267	! #wH debugH( 1:35) = 'Lin et al.(1983): Accretion of Clou'
	2268	! #wH debugH(36:70) = 'd Ice by Snow Particles '
	2269	! #wH proc_1 = 'Qsaci g/kg'
	2270	! #wH procv1 = qs_ACI
	2271	! #wH proc_2 = 'CNsaci/e15'
	2272	! #wH procv2 = CNsACI*1.e-18
	2273	! #wH proc_3 = ' '
	2274	! #wH procv3 = 0.
	2275	! #wH proc_4 = ' '
	2276	! #wH procv4 = 0.
	2277	! #wh include 'CMiPhy_Debug.h'
	2278	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2279	! #wH& debugV(k,10) = qs_ACI
	2280
	2281	END DO
	2282
	2283
	2284
	2285
	2286	! Accretion of Cloud Ice by Graupel (Cloud Ice Sink) ACCRETION, # > /
	2287	! Reference: Lin et al. 1983, JCAM 22, p.1075 (41) ++++++++++++++++
	2288	! Emde and Kahlig 1989, Ann.Geoph. 7, p. 407 (~19)
	2289	! -----------------------------------------------------------
	2290
	2291	! #qg DO k=mz1_CM,mzp
	2292
	2293	! #qg IF (qi__CM(ikl,k).gt.epsn) THEN
	2294	! #qg CbyG_c = max(zer0, sign(un_1,qi__CM(ikl,k) - epsn))
	2295	! CbyG_c = 1.0 if qi__CM(ikl,k) > epsn
	2296	! = 0.0 otherwise
	2297
	2298	! #qg IF (qg__CM(ikl,k).gt.epsn) THEN
	2299	! #qg CbyG_g = max(zer0, sign(un_1,qg__CM(ikl,k) - epsn))
	2300	! CbyG_g = 1.0 if qg__CM(ikl,k) > epsn
	2301	! = 0.0 otherwise
	2302
	2303	! #qg IF (Ta__CM(ikl,k).lt.Tf_Sno) THEN
	2304	! #qg Fact_G = max(zer0,-sign(un_1,Ta__CM(ikl,k) - Tf_Sno))
	2305	! Fact_G = 1.0 if Ta__CM(ikl,k) < Tf_Sno
	2306	! = 0.0 otherwise
	2307
	2308	! #qg CbyGOK = CbyG_c * CbyG_g * Fact_G
	2309
	2310	! #qg pgaci = PATATRAS
	2311	! #qg qgaci = pgaci dt__CM CbyGOK
	2312	! #qg qgaci = min(qgaci,qi__CM(ikl,k))
	2313
	2314	! #qg qi__CM(ikl,k) = qi__CM(ikl,k) - qgaci
	2315	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qgaci
	2316
	2317	! #qg END IF
	2318	! #qg END IF
	2319	! #qg END IF
	2320
	2321	! #qg END DO
	2322
	2323
	2324	! Accretion of Cloud Ice by Rain (Cloud Ice Sink) ACCRETION, o > / \| o
	2325	! Reference: Lin et al. 1983, JCAM 22, p.1071 (25) ++++++++++++++++++++
	2326	! ----------------------------------------------------------
	2327
	2328	DO k=mz1_CM,mzp
	2329
	2330	! #wH qr_ACI = 0.0
	2331	! #wH qi_ACR = 0.0
	2332
	2333	! #hy IF (qi__CM(ikl,k).gt.epsn) THEN
	2334	qi__OK = max(zer0,sign(un_1,qi__CM(ikl,k) - epsn))
	2335	! qi__OK = 1.0 if qi__CM(ikl,k) > epsn
	2336	! = 0.0 otherwise
	2337
	2338	! #hy IF (qr___0(ikl,k).gt.epsn) THEN
	2339	qr0_OK = max(zer0,sign(un_1,qr___0(ikl,k) - epsn))
	2340	! qr0_OK = 1.0 if qr___0(ikl,k) > epsn
	2341	! = 0.0 otherwise
	2342
	2343	! #hy IF (Ta__CM(ikl,k).lt.Tf_Sno) THEN
	2344	Flag_Ta_Neg = max(zer0,-sign(un_1,Ta__CM(ikl,k) - Tf_Sno))
	2345	! Flag_Ta_Neg = 1.0 if Ta__CM(ikl,k) < Tf_Sno
	2346	! = 0.0 otherwise
	2347
	2348	Flag_qr_ACI = qi__OK * qr0_OK * Flag_Ta_Neg
	2349
	2350	! #EW IF(Flag_qr_ACI.gt.eps6) THEN
	2351	! #EW mauxEW = mphyEW(ikl)
	2352	! #EW IF (mauxEW(13:13).eq.'s'.or.mauxEW(13:13).eq.'A') THEN
	2353	! #EW mauxEW(13:13) = 'A'
	2354	! #EW ELSE
	2355	! #EW mauxEW(13:13) = 'r'
	2356	! #EW END IF
	2357	! #EW mphyEW(ikl) = mauxEW
	2358	! #EW END IF
	2359
	2360	pr_ACI = 3104.28d0 * n0___r * sqrrro(ikl,k) &!
	2361	& qi__CM(ikl,k) /exp(3.8d0log(lamdaR(ikl,k)))
	2362	qr_ACI = pr_ACIdt__CM Flag_qr_ACI
	2363	qr_ACI = min(qr_ACI,qi__CM(ikl,k))
	2364	CNrACI = CCNiCM(ikl,k)* qr_ACI/max(qi__CM(ikl,k),epsn)
	2365	CCNiCM(ikl,k) = CCNiCM(ikl,k)- CNrACI
	2366	qi__CM(ikl,k) = qi__CM(ikl,k)- qr_ACI
	2367
	2368	! #qg IF(qr__CM(ikl,k) .gt. 1.e-4 ) THEN
	2369	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qr_ACI
	2370	! CAUTION : Graupels Formation is not taken into account
	2371	! This could be a reasonable assumption for Antarctica
	2372
	2373	! #qg ELSE
	2374	qs__CM(ikl,k) = qs__CM(ikl,k) + qr_ACI
	2375	! #qg END IF
	2376
	2377	! #WQ write(6,*) 'Qraci',qr_ACI,it_EXP,ikl,k
	2378	! #WH if (ikl.eq.ikl0CM(1)) wraci(k) = qr_ACI
	2379
	2380
	2381	! Accretion of Rain by Cloud Ice (Rain Sink) ACCRETION, / > o \| *
	2382	! Reference: Lin et al. 1983, JCAM 22, p.1071 (26) ++++++++++++++++++++
	2383	! ----------------------------------------------------------
	2384
	2385	! #EW IF (Flag_qr_ACI.gt.eps6) THEN
	2386	! #EW mauxEW = mphyEW(ikl)
	2387	! #EW IF (mauxEW(13:13).eq.'r'.or.mauxEW(13:13).eq.'A') THEN
	2388	! #EW mauxEW(13:13) = 'A'
	2389	! #EW ELSE
	2390	! #EW mauxEW(13:13) = 's'
	2391	! #EW END IF
	2392	! #EW mphyEW(ikl) = mauxEW
	2393	! #EW END IF
	2394
	2395	pi_ACR = 4.1d20 * n0___r * sqrrro(ikl,k) &! 4.1e20 = a pi**2 rhow/mi Gamma[6+b] / 24
	2396	& qi__CM(ikl,k) /exp(6.8d0log(lamdaR(ikl,k))) ! where a=842., rhow=1000, mi=4.19e-13
	2397	! b = 0.8
	2398	! Lin et al, 1983, JAM,p1071: mi:Ice Crystal Mass
	2399	qi_ACR = pi_ACRdt__CM Flag_qr_ACI
	2400	qi_ACR = min(qi_ACR,qr__CM(ikl,k))
	2401	qr__CM(ikl,k) = qr__CM(ikl,k) - qi_ACR
	2402	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Lc_Cpd *qi_ACR
	2403
	2404	! #qg IF (qr__CM(ikl,k) .gt. 1.e-4 ) THEN
	2405	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qi_ACR
	2406	! CAUTION : Graupels Formation is not taken into account
	2407	! This could be a reasonable assumption for Antarctica
	2408
	2409	! #qg ELSE
	2410	qs__CM(ikl,k) = qs__CM(ikl,k) + qi_ACR
	2411	! #qg END IF
	2412
	2413	! Full Debug
	2414	! ~~~~~~~~~~
	2415	! #WQ write(6,*) 'Qiacr',qi_ACR,it_EXP,ikl,k
	2416	! #WH if (ikl.eq.ikl0CM(1)) wiacr(k) = qi_ACR
	2417
	2418	! #hy END IF
	2419	! #hy END IF
	2420	! #hy END IF
	2421
	2422	! Debug
	2423	! ~~~~~
	2424	! #wH debugH( 1:35) = 'Lin et al.(1983): Accretion of Clou'
	2425	! #wH debugH(36:70) = 'd Ice by Rain '
	2426	! #wH proc_1 = 'Qraci g/kg'
	2427	! #wH procv1 = qr_ACI
	2428	! #wH proc_2 = 'qi_ACR g/kg'
	2429	! #wH procv2 = qi_ACR
	2430	! #wH proc_3 = ' '
	2431	! #wH procv3 = 0.
	2432	! #wH proc_4 = ' '
	2433	! #wH procv4 = 0.
	2434	! #wh include 'CMiPhy_Debug.h'
	2435	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2436	! #wH& debugV(k,11) = qi_ACR
	2437
	2438	END DO
	2439
	2440
	2441
	2442
	2443	! Accretion of Rain by Snow Flakes ACCRETION o > , > o
	2444	! Accretion of Snow Flakes by Rain ++++++++++++++++++++++
	2445	! Reference: Lin et al. 1983, JCAM 22, p.1071 (27)
	2446	! Lin et al. 1983, JCAM 22, p.1071 (28)
	2447	! Emde and Kahlig 1989, Ann.Geoph. 7, p. 408 (~21)
	2448	! -----------------------------------------------------------
	2449
	2450	DO k=mz1_CM,mzp
	2451
	2452	ps_ACR(ikl,k) = 0.0
	2453	qs_ACR = 0.0
	2454	qs_ACR_S = 0.0
	2455	qs_ACR_R = 0.0
	2456
	2457	! #hy IF (qr___0(ikl,k).gt.epsn) THEN
	2458	qr0_OK = max(zer0,sign(un_1,qr___0(ikl,k) - epsn))
	2459	! qr0_OK = 1.0 if qr___0(ikl,k) > epsn
	2460	! = 0.0 otherwise
	2461
	2462	! #hy IF (qs___0(ikl,k).gt.epsn) THEN
	2463	qs0_OK = max(zer0,sign(un_1,qs___0(ikl,k) - epsn))
	2464	! qs0_OK = 1.0 if qs___0(ikl,k) > epsn
	2465	! = 0.0 otherwise
	2466
	2467	Flag_qr_ACS = qr0_OK * qs0_OK
	2468
	2469	! #EW IF(Flag_qr_ACI.gt.eps6) THEN
	2470	! #EW mauxEW = mphyEW(ikl)
	2471	! #EW mauxEW(14:14) = 'A'
	2472	! #EW mphyEW(ikl) = mauxEW
	2473	! #EW END IF
	2474
	2475	! Accretion of Rain by Snow --> Snow \| lamdaR : lambda_r
	2476	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \| lamdaS : lambda_s
	2477	coeACS=(5.0d0/(lamdaS(ikl,k)lamdaS(ikl,k)lamdaR(ikl,k)) &!
	2478	& +2.0d0/(lamdaS(ikl,k)lamdaR(ikl,k)lamdaR(ikl,k)) &!
	2479	& +0.5d0/(lamdaR(ikl,k)lamdaR(ikl,k)lamdaR(ikl,k))) &!
	2480	& /(lamdaS(ikl,k)lamdaS(ikl,k)lamdaS(ikl,k)*lamdaS(ikl,k)) !
	2481
	2482	! #cn n0___s = min(2.e8,2.e6exp(-.12 min(0.,Ta_dgC(ikl,k))))
	2483
	2484	pr_ACS = 986.96d-3(n0___rn0___s/roa_DY(ikl,k)) &! 986.96: pi*2 rhos
	2485	& * abs(FallVr(ikl,k)-FallVs(ikl,k))*coeACS ! (snow density assumed equal to 100 kg/m3)
	2486	qr_ACS = pr_ACS dt__CM Flag_qr_ACS
	2487	qr_ACS = min(qr_ACS, qr__CM(ikl,k))
	2488
	2489	! #WQ write(6,*) 'Qracs',qr_ACS,it_EXP,ikl,k
	2490	! #WH if (ikl.eq.ikl0CM(1)) wracs(k) = qr_ACS
	2491
	2492	! Accretion of Snow by Rain --> Rain
	2493	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2494	qr0_OK = max(zer0,sign(un_1,qr___0(ikl,k) - 1.e-4))
	2495	! qr0_OK = 1.0 if qr___0(ikl,k) > 1.e-4
	2496	! = 0.0 otherwise
	2497
	2498	qs0_OK = max(zer0,sign(un_1,qs___0(ikl,k) - 1.e-4))
	2499	! qs0_OK = 1.0 if qs___0(ikl,k) > 1.e-4
	2500	! = 0.0 otherwise
	2501
	2502	Flag_qs_ACR = max(qr0_OK,qs0_OK)
	2503
	2504	! #hy IF (Flag_qs_ACR.gt.eps6) THEN
	2505	coeACR=(5.0d0/(lamdaR(ikl,k)lamdaR(ikl,k)lamdaS(ikl,k)) &
	2506	& +2.0d0/(lamdaR(ikl,k)lamdaS(ikl,k)lamdaS(ikl,k)) &
	2507	& +0.5d0/(lamdaS(ikl,k)lamdaS(ikl,k)lamdaS(ikl,k))) &
	2508	& /(lamdaR(ikl,k) lamdaR(ikl,k)lamdaR(ikl,k)*lamdaR(ikl,k))
	2509
	2510	ps_ACR(ikl,k)=9869.6d-3(n0___rn0___s/roa_DY(ikl,k)) &! 9869.6: pi*2 rhow
	2511	& * abs(FallVr(ikl,k)-FallVs(ikl,k))*coeACR! (water density assumed equal to 1000 kg/m3)
	2512	qs_ACR = ps_ACR(ikl,k)dt__CM Flag_qr_ACS *Flag_qs_ACR
	2513	qs_ACR = min(qs_ACR,qs__CM(ikl,k))
	2514
	2515	! #WQ write(6,*) 'Qsacr',qs_ACR,it_EXP,ikl,k
	2516	! #WH if (ikl.eq.ikl0CM(1)) wsacr(k) = qs_ACR
	2517	! #hy ELSE
	2518	! #hy ps_ACR(ikl,k) = 0.d0
	2519	! #hy qs_ACR = 0.d0
	2520	! #hy END IF
	2521
	2522	Flag_Ta_Neg = max(zer0,-sign(un_1,Ta__CM(ikl,k) - Tf_Sno))
	2523	! Flag_Ta_Neg = 1.0 if Ta__CM(ikl,k) < Tf_Sno
	2524	! = 0.0 otherwise
	2525
	2526	qr_ACS_S = qr_ACS * Flag_Ta_Neg
	2527	qs_ACR_R = qs_ACR *(1.d0-Flag_Ta_Neg)
	2528	qr__CM(ikl,k) = qr__CM(ikl,k) - qr_ACS_S
	2529	! #qg IF (qr___0(ikl,k).lt.1.e-4 .and. qs___0(ikl,k).lt.1.e-4)THEN
	2530	! CAUTION : Graupel Formation is not taken into Account
	2531	qs__CM(ikl,k) = qs__CM(ikl,k) + qr_ACS_S
	2532	! #qg ELSE²
	2533	! #qg qs__CM(ikl,k) = qs__CM(ikl,k) - qr_ACS_S
	2534	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qs_ACR_S + qr_ACS_S
	2535	! #qg REND IF
	2536	Ta__CM(ikl,k) = Ta__CM(ikl,k) + qs_ACR_S * Lc_Cpd
	2537
	2538	qr__CM(ikl,k) = qr__CM(ikl,k) + qs_ACR_R
	2539	qs__CM(ikl,k) = qs__CM(ikl,k) - qs_ACR_R
	2540	Ta__CM(ikl,k) = Ta__CM(ikl,k) - qs_ACR_R * Lc_Cpd
	2541
	2542	! #hy END IF
	2543	! #hy END IF
	2544
	2545	! Debug
	2546	! ~~~~~
	2547	! #wH debugH( 1:35) = 'Lin et al.(1983): Accretion of Snow'
	2548	! #wH debugH(36:70) = '(Rain) by Rain(Snow) '
	2549	! #wH proc_1 = 'Qracs g/kg'
	2550	! #wH procv1 = qs_ACR_S
	2551	! #wH proc_2 = 'Qsacr g/kg'
	2552	! #wH procv2 = qs_ACR_R
	2553	! #wH proc_3 = ' '
	2554	! #wH procv3 = 0.
	2555	! #wH proc_4 = ' '
	2556	! #wH procv4 = 0.
	2557	! #wh include 'CMiPhy_Debug.h'
	2558	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2559	! #wH& debugV(k,12) = qs_ACR_S - qs_ACR_R
	2560
	2561	END DO
	2562
	2563
	2564
	2565
	2566	! Accretion of Snow by Graupels ACCRETION, # > *
	2567	! Reference: Lin et al. 1983, JCAM 22, p.1071 (29) ++++++++++++++++
	2568	! ----------------------------------------------------------
	2569
	2570	! #qg DO k=mz1_CM,mzp
	2571
	2572	! #qg IF (qg___0(ikl,k).gt.epsn) THEN
	2573	! #qg SbyG_g = max(zer0,sign(un_1,qg___0(ikl,k) - epsn))
	2574	! SbyG_g = 1.0 if qg___0(ikl,k) > epsn
	2575	! = 0.0 otherwise
	2576
	2577	! #qg IF (qs___0(ikl,k).gt.epsn) THEN
	2578	! #qg SbyG_s = max(zer0,sign(un_1,qs___0(ikl,k) - epsn))
	2579	! SbyG_s = 1.0 if qs___0(ikl,k) > epsn
	2580	! = 0.0 otherwise
	2581
	2582	! #qg SbyGOK = SbyG_g * SbyG_s
	2583	! #qg effACS = exp(0.090*Ta_dgC(ikl,k)) &! Collection Efficiency
	2584	! ! Lin et al. 1983 JCAM 22 p.1072 (30)
	2585
	2586	! #qg flg=exp(-6.0d0*log(lamdaS(ikl,k)) &!
	2587	! #qg& *(5.0/lamdaG(ikl,k) &!
	2588	! #qg& +2.0lamdaS(ikl,k)/(lamdaG(ikl,k)lamdaG(ikl,k)) &!
	2589	! #qg& +0.5lamdaS(ikl,k) lamdaS(ikl,k) &!
	2590	! #qg& /exp(3.0d0*log(lamdaG(ikl,k)))) !
	2591
	2592	! #cn n0___s = min(2.e8,2.e6exp(-.12min(0.,Ta_dgC(ikl,k))))
	2593
	2594	! #qg pgacs = 986.96d-3(n0___gn0___s/roa_DY(ikl,k)) &! 986.96: pi*2 rhog
	2595	! #qg& * abs(FallVg(ikl,k)-FallVs(ikl,k))flgeffACS !(graupel densitity assumed equal to snow density)
	2596	! #qg qgacs = pgacsdt__CM SbyGOK
	2597	! #qg qgacs = min(qgacs,qs__CM(ikl,k))
	2598	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qgacs
	2599	! #qg qs__CM(ikl,k) = qs__CM(ikl,k) - qgacs
	2600
	2601	! #qg END IF
	2602	! #qg END IF
	2603
	2604	! #qg END DO
	2605
	2606
	2607	! Accretion of Rain by Graupels (Dry Growth Mode) ACCRETION, # > o
	2608	! Reference: Lin et al. 1983, JCAM 22, p.1075 (42) ++++++++++++++++
	2609	! ----------------------------------------------------------
	2610
	2611	! #qg DO k=mz1_CM,mzp
	2612
	2613	! #qg IF (qg___0(ikl,k).gt.epsn) THEN
	2614	! #qg RbyG_g = max(zer0,sign(un_1,qg___0(ikl,k) - epsn))
	2615	! RbyG_g = 1.0 if qg___0(ikl,k) > epsn
	2616	! = 0.0 otherwise
	2617
	2618	! #qg IF (qr___0(ikl,k).gt.epsn) THEN
	2619	! #qg RbyG_r = max(zer0,sign(un_1,qr___0(ikl,k) - epsn))
	2620	! RbyG_r = 1.0 if qr___0(ikl,k) > epsn
	2621	! = 0.0 otherwise
	2622
	2623	! #qg IF (Ta__CM(ikl,k).lt.Tf_Sno) THEN
	2624	! #qg Fact_G = max(zer0,-sign(un_1,Ta__CM(ikl,k) - Tf_Sno))
	2625	! Fact_G = 1.0 if Ta__CM(ikl,k) < Tf_Sno
	2626	! = 0.0 otherwise
	2627
	2628	! #qg RbyGOK = RbyG_g * RbyG_s * Fact_G
	2629
	2630	! #qg flg=exp(-6.0d0*log(lamdaS(ikl,k)) &!
	2631	! #qg& *(5.0/lamdaG(ikl,k) &!
	2632	! #qg& +2.0lamdaS(ikl,k)/(lamdaG(ikl,k)lamdaG(ikl,k)) &!
	2633	! #qg& +0.5lamdaS(ikl,k) lamdaS(ikl,k) &!
	2634	! #qg& /exp(3.0d0*log(lamdaG(ikl,k)))) !
	2635
	2636	! #cn n0___s = min(2.e8,2.e6exp(-.12min(0.,Ta_dgC(ikl,k))))
	2637
	2638	! #qg pgacr = 986.96d-3(n0___gn0___s/roa_DY(ikl,k)) &!
	2639	! #qg& * abs(FallVg(ikl,k)-FallVr(ikl,k))*flg
	2640	! #qg qgacr = pgacr * dt__CM *RbyGOK
	2641	! #qg qgacr = min(qgacr,qr__CM(ikl,k))
	2642	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qgacr
	2643	! #qg qr__CM(ikl,k) = qr__CM(ikl,k) - qgacr
	2644	! #qg Ta__CM(ikl,k) = Ta__CM(ikl,k) + Lc_Cpd*qgacr
	2645
	2646	! #qg END IF
	2647	! #qg END IF
	2648	! #qg END IF
	2649
	2650	! #qg END DO
	2651
	2652
	2653	! Graupels Wet Growth Mode
	2654	! Reference: Lin et al. 1983, JCAM 22, p.1075 (43)
	2655	! ----------------------------------------------------------
	2656
	2657	! #qg ! TO BE ADDED !
	2658
	2659
	2660
	2661
	2662	! Microphysical Processes affecting Precipitating Cloud Particles
	2663	! ===================================================================
	2664
	2665
	2666	! Rain Drops Evaporation RAIN, EVAPORATION
	2667	! Reference: Lin et al. 1983, JCAM 22, p.1077 (52) +++++++++++++++++
	2668	! ----------------------------------------------------------
	2669
	2670	DO k=mz1_CM,mzp
	2671
	2672	! #wH qr_EVP = 0.0
	2673
	2674	! #hy IF (qr___0(ikl,k).gt.epsn) THEN
	2675	qr0_OK = max(zer0,sign(un_1,qr___0(ikl,k) - epsn))
	2676	! qr0_OK = 1.0 if qr___0(ikl,k) > epsn
	2677	! = 0.0 otherwise
	2678
	2679	! #EW IF(qr0_OK.gt.eps6) THEN
	2680	! #EW mauxEW = mphyEW(ikl)
	2681	! #EW mauxEW(15:15) = 'v'
	2682	! #EW mphyEW(ikl) = mauxEW
	2683	! #EW END IF
	2684
	2685	RH_Liq = qv__DY(ikl,k)/(RHcrit*qvswCM(ikl,k))
	2686	! RH_Liq : grid scale saturation humidity
	2687
	2688	! #hy IF (RH_Liq.lt.un_1) THEN
	2689	Flag_DryAir = max(zer0,-sign(un_1,RH_Liq - un_1))
	2690	! Flag_DryAir = 1.0 if RH_Liq < un_1
	2691	! = 0.0 otherwise
	2692
	2693	Flag_qr_EVP = qr0_OK * Flag_DryAir
	2694
	2695	lr_NUM = 0.78d0 /(lamdaR(ikl,k) *lamdaR(ikl,k)) &!
	2696	& + 3940.d0 * sqrt(sqrrro(ikl,k)) &! 3940.: 0.31 Sc*(1/3) (a/nu)*(1/2) Gamma[(b+5)/2]
	2697	& /exp(2.9d0 *log(lamdaR(ikl,k))) ! where Sc=0.8(Schm.) nu=1.5e-5 (Air Kinematic Viscosity)
	2698	lr_DEN = 5.423d11/(Ta__CM(ikl,k) Ta__CM(ikl,k)) &! 5.423e11 = [Lv=2500000J/kg] Lv / [kT=0.025W/m/K] / [Rv=461.J/kg/K]
	2699	& + 1.d0/(1.875d-2* roa_DY(ikl,k) *qvswCM(ikl,k)) ! kT: Air Thermal Conductivity
	2700
	2701	pr_EVP = 2. piNmbr(1.d0 -RH_Liq)n0___rlr_NUM/lr_DEN
	2702	qr_EVP = pr_EVP* dt__CM
	2703	qr_EVP = min(qr_EVP, qr__CM(ikl,k))
	2704
	2705	qr_EVP = min(qr_EVP,RHcrit*qvswCM(ikl,k) -qv__DY(ikl,k))
	2706	! supersaturation is not allowed to occur
	2707
	2708	qr_EVP = max(qr_EVP,zer0) * Flag_qr_EVP
	2709	! condensation is not allowed to occur
	2710
	2711	qr__CM(ikl,k) = qr__CM(ikl,k) - qr_EVP
	2712	qwd_CM(ikl,k) = qwd_CM(ikl,k) - qr_EVP
	2713	qv__DY(ikl,k) = qv__DY(ikl,k) + qr_EVP
	2714	Ta__CM(ikl,k) = Ta__CM(ikl,k) - Lv_Cpd *qr_EVP
	2715
	2716	! Full Debug
	2717	! ~~~~~~~~~~
	2718	! #WQ write(6,*) 'Qrevp',qr_EVP,it_EXP,ikl,k
	2719	! #WH if (ikl.eq.ikl0CM(1)) wrevp(k) = qr_EVP
	2720
	2721	! #hy END IF
	2722	! #hy END IF
	2723
	2724	! Debug
	2725	! ~~~~~
	2726	! #wH debugH( 1:35) = 'Lin et al.(1983): Rain Drops Evapor'
	2727	! #wH debugH(36:70) = 'ation '
	2728	! #wH proc_1 = 'Qrevp g/kg'
	2729	! #wH procv1 = qr_EVP
	2730	! #wH proc_2 = 'R.Hum [%]'
	2731	! #wH procv2 = RH_Liq*0.1
	2732	! #wH proc_3 = ' '
	2733	! #wH procv3 = 0.
	2734	! #wH proc_4 = ' '
	2735	! #wH procv4 = 0.
	2736	! #wh include 'CMiPhy_Debug.h'
	2737	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2738	! #wH& debugV(k,13) = -qr_EVP
	2739
	2740	END DO
	2741
	2742
	2743	! (Deposition on) Snow Flakes (Sublimation) SNOW, (SUBLI)/(DEPOSI)TION
	2744	! Reference: Lin et al. 1983, JCAM 22, p.1072 (31) ++++++++++++++++++++++++++
	2745	! ----------------------------------------------------------
	2746
	2747	DO k=mz1_CM,mzp
	2748
	2749	! #wH qs_SUB = 0.0
	2750
	2751	! #hy IF (qs___0(ikl,k).gt.epsn) THEN
	2752	qs0_OK = max(zer0,sign(un_1,qs___0(ikl,k) - epsn))
	2753	! qs0_OK = 1.0 if qs___0(ikl,k) > epsn
	2754	! = 0.0 otherwise
	2755
	2756	! #EW IF(qs0_OK.gt.eps6) THEN
	2757	! #EW mauxEW = mphyEW(ikl)
	2758	! #EW mauxEW(16:16) = 'V'
	2759	! #EW mphyEW(ikl) = mauxEW
	2760	! #EW END IF
	2761
	2762	RH_ICE = qv__DY(ikl,k)/qsiEFF(ikl,k)
	2763
	2764	ls_NUM = 0.78d0 /(lamdaS(ikl,k)*lamdaS(ikl,k)) &!
	2765	& + 238.d0 * sqrt(sqrrro(ikl,k)) &! 238.: 0.31 Sc*(1/3) (c/nu)*(1/2) Gamma[(d+5)/2]
	2766	& /exp(2.625d0*log(lamdaS(ikl,k))) ! where Sc=0.8(Schm.) nu=1.5e-5 (Air Kinematic Viscosity)
	2767	ls_DEN = 6.959d11/(Ta__CM(ikl,k)Ta__CM(ikl,k)) &! 6.959e11 = [Ls=2833600J/kg]Ls /[kT=0.025W/m/K] /[Rv=461.J/kg/K]
	2768	& + 1.d0/(1.875d-2roa_DY(ikl,k)qsiEFF(ikl,k)) ! kT: Air Thermal Conductivity
	2769
	2770	! #cn n0___s = min(2.e8,2.e6exp(-.12min(0.,Ta_dgC(ikl,k))))
	2771
	2772	ps_SUB = 2piNmbr(1.d0-RH_ICE)n0___sls_NUM &!
	2773	& /(1.d3roa_DY(ikl,k)ls_DEN)
	2774	qs_SUB = ps_SUB * dt__CM
	2775
	2776	dqsiqv = qsiEFF(ikl,k) -qv__DY(ikl,k)
	2777
	2778	Flag_SURSat = max(zer0,sign(un_1,RH_ICE - un_1))
	2779	! Flag_SURSat = 1.0 if RH_ICE > un_1
	2780	! = 0.0 otherwise
	2781
	2782	qs_SUB = max(qs_SUB ,dqsiqv)* Flag_SURSat &! qs_SUB < 0 ... Deposition
	2783	& + min(min(qs_SUB,qs__CM(ikl,k)),dqsiqv)*(1.-Flag_SURSat) ! > 0 ... Sublimation
	2784
	2785	qs_SUB = qs_SUB * qs0_OK
	2786
	2787	qs__CM(ikl,k) = qs__CM(ikl,k)- qs_SUB
	2788	qid_CM(ikl,k) = qid_CM(ikl,k)- qs_SUB
	2789	qv__DY(ikl,k) = qv__DY(ikl,k)+ qs_SUB
	2790	Ta__CM(ikl,k) = Ta__CM(ikl,k)-Ls_Cpd *qs_SUB
	2791
	2792	! Full Debug
	2793	! ~~~~~~~~~~
	2794	! #WQ write(6,*) 'Qssub',qs_SUB,it_EXP,ikl,k
	2795	! #WH if (ikl.eq.ikl0CM(1)) wssub(k) = -qs_SUB
	2796
	2797	! #hy END IF
	2798
	2799	! Debug
	2800	! ~~~~~
	2801	! #wH debugH( 1:35) = 'Lin et al.(1983): (Deposition on) S'
	2802	! #wH debugH(36:70) = 'now Particles (Sublimation) '
	2803	! #wH proc_1 = 'Qssub g/kg'
	2804	! #wH procv1 = qs_SUB
	2805	! #wH proc_2 = ' '
	2806	! #wH procv2 = 0.
	2807	! #wH proc_3 = ' '
	2808	! #wH procv3 = 0.
	2809	! #wH proc_4 = ' '
	2810	! #wH procv4 = 0.
	2811	! #wh include 'CMiPhy_Debug.h'
	2812	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2813	! #wH& debugV(k,14) = -qs_SUB
	2814
	2815	END DO
	2816
	2817
	2818	! Graupels Sublimation GRAUPEL, SUBLIMATION
	2819	! Reference: Lin et al. 1983, JCAM 22, p.1076 (46) ++++++++++++++++++++
	2820	! ----------------------------------------------------------
	2821
	2822	! #qg ! TO BE ADDED !
	2823
	2824
	2825	! Snow Flakes Melting PSMLT SNOW, MELT
	2826	! Reference: Lin et al. 1983, JCAM 22, p.1072 (32) ++++++++++
	2827	! ----------------------------------------------------------
	2828
	2829	DO k=mz1_CM,mzp
	2830
	2831	! #wH qsMELT = 0.0
	2832
	2833	! #hy IF (qs___0(ikl,k).gt.epsn) THEN
	2834	qs0_OK = max(zer0,sign(un_1,qs___0(ikl,k) - epsn))
	2835	! qs0_OK = 1.0 if qs___0(ikl,k) > epsn
	2836	! = 0.0 otherwise
	2837
	2838	! #hy IF (Ta_dgC(ikl,k).gt.0.) THEN! Ta_dgC : old Celsius Temperature
	2839	Flag_Ta_Pos = max(zer0,sign(un_1,Ta_dgC(ikl,k) - 0.))
	2840	! Flag_Ta_Pos = 1.0 if Ta_dgC(ikl,k) > 0.
	2841	! = 0.0 otherwise
	2842
	2843	Flag_qsMELT = qs0_OK * Flag_Ta_Pos
	2844
	2845	! #EW IF(Flag_qsMELT.gt.eps6) THEN
	2846	! #EW mauxEW = mphyEW(ikl)
	2847	! #EW mauxEW(17:17) = 'r'
	2848	! #EW mphyEW(ikl) = mauxEW
	2849	! #EW END IF
	2850
	2851	ls_NUM = 0.78 / (lamdaS(ikl,k) *lamdaS(ikl,k)) &
	2852	& + 238. * sqrt(sqrrro(ikl,k)) &
	2853	& / exp(2.625d0 *log(lamdaS(ikl,k)))
	2854
	2855	! #cn n0___s = min(2.e8,2.e6exp(-.12min(0.,Ta_dgC(ikl,k))))
	2856
	2857	xCoefM = 1.904d-8 n0___s ls_NUM *Lc_Cpd /roa_DY(ikl,k) ! 1.904e-8: 2 pi / Lc /[1.e3=rho Factor]
	2858
	2859	AcoefM = 0.025d00 *xCoefM &!
	2860	& +(ps_ACW(ikl,k) + ps_ACR(ikl,k)) *Lc_Cpd /78.8d0 ! 78.8 : Lc /[Cpw=4.187e3 J/kg/K]
	2861
	2862	BcoefM = 62.34d+3 roa_DY(ikl,k) &! 62.34 : Ls [psiv=2.200e-5 m2/s]
	2863	& (qv__DY(ikl,k)-qsiEFF(ikl,k)) &! 46.88 : Lv [psiv=1.875e-5 m2/s]
	2864	& *xCoefM
	2865	BcoefM = min(-epsn,BcoefM)
	2866
	2867	dTMELT = ( Ta__CM(ikl,k) -Tf_Sno -AcoefM/BcoefM) &!
	2868	& exp(-AcoefMdt__CM) !
	2869	qsMELT = ( Ta__CM(ikl,k) -Tf_Sno -dTMELT )/ Lc_Cpd !
	2870	qsMELT = max( qsMELT,0. ) *Flag_qsMELT !
	2871	qsMELT = min( qsMELT,qs__CM(ikl,k)) !
	2872
	2873	qs__CM(ikl,k) = qs__CM(ikl,k) - qsMELT
	2874	qr__CM(ikl,k) = qr__CM(ikl,k) + qsMELT
	2875	Ta__CM(ikl,k) = Ta__CM(ikl,k) - Lc_Cpd *qsMELT
	2876
	2877	! Full Debug
	2878	! ~~~~~~~~~~
	2879	! #WQ write(6,*) 'Qsmlt',qsMELT,it_EXP,ikl,k
	2880	! #WH if (ikl.eq.ikl0CM(1)) wsmlt(k) = qsMELT
	2881
	2882	! #hy END IF
	2883	! #hy END IF
	2884
	2885	! Debug
	2886	! ~~~~~
	2887	! #wH debugH( 1:35) = 'Lin et al.(1983): Snow Particles Me'
	2888	! #wH debugH(36:70) = 'lting '
	2889	! #wH proc_1 = 'Qsmlt g/kg'
	2890	! #wH procv1 = qsMELT
	2891	! #wH proc_2 = ' '
	2892	! #wH procv2 = 0.
	2893	! #wH proc_3 = ' '
	2894	! #wH procv3 = 0.
	2895	! #wH proc_4 = ' '
	2896	! #wH procv4 = 0.
	2897	! #wh include 'CMiPhy_Debug.h'
	2898	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2899	! #wH& debugV(k,15) = -qsMELT
	2900
	2901	END DO
	2902
	2903
	2904	! Graupels Melting GRAUPELS, MELT
	2905	! Reference: Lin et al. 1983, JCAM 22, p.1076 (47) ++++++++++++++
	2906	! ----------------------------------------------------------
	2907
	2908	! #qg ! TO BE ADDED !
	2909
	2910
	2911	! Rain Freezing RAIN, FREEZING
	2912	! Reference: Lin et al. 1983, JCAM 22, p.1075 (45) ++++++++++++++
	2913	! ----------------------------------------------------------
	2914
	2915	! CAUTION: Graupel Formation TO BE ADDED !
	2916
	2917	DO k=1,mzp
	2918
	2919	! #wH qs_FRZ = 0.0
	2920
	2921	! #hy IF (qr___0(ikl,k).gt.epsn) THEN
	2922	qr0_OK = max(zer0,sign(un_1,qr___0(ikl,k) - epsn))
	2923	! qr0_OK = 1.0 if qr___0(ikl,k) > epsn
	2924	! = 0.0 otherwise
	2925
	2926	! #hy IF (Ta_dgC(ikl,k).lt.0.e0)THEN! Ta_dgC : old Celsius Temperature
	2927	Flag_Ta_Neg = max(zer0,-sign(un_1,Ta_dgC(ikl,k) - 0.e0))
	2928	! Flag_Ta_Neg = 1.0 if Ta_dgC(ikl,k) < 0.e0
	2929	! = 0.0 otherwise
	2930
	2931	Flag_Freeze = qr0_OK * Flag_Ta_Neg
	2932
	2933	! #EW IF(Flag_Freeze.gt.eps6) THEN
	2934	! #EW mauxEW = mphyEW(ikl)
	2935	! #EW mauxEW(19:19) = 's'
	2936	! #EW mphyEW(ikl) = mauxEW
	2937	! #EW END IF
	2938
	2939	ps_FRZ = 1.974d4 *n0___r &
	2940	& /(roa_DY(ikl,k)exp(7.d0 log(lamdaR(ikl,k)))) &
	2941	& (exp(-0.66d0 Ta_dgC(ikl,k))-1.d0)
	2942	qs_FRZ = ps_FRZ * dt__CM *Flag_Freeze
	2943	qs_FRZ = min(qs_FRZ,qr__CM(ikl,k))
	2944
	2945	qr__CM(ikl,k) = qr__CM(ikl,k) - qs_FRZ
	2946	qs__CM(ikl,k) = qs__CM(ikl,k) + qs_FRZ
	2947	! CAUTION : graupel production is included into snow production
	2948	! proposed modification in line below.
	2949	! #qg qg__CM(ikl,k) = qg__CM(ikl,k) + qs_FRZ
	2950	Ta__CM(ikl,k) = Ta__CM(ikl,k) + Lc_Cpd *qs_FRZ
	2951
	2952	! Full Debug
	2953	! ~~~~~~~~~~
	2954	! #WQ write(6,*) 'Qsfre',qs_FRZ,it_EXP,ikl,k
	2955	! #WH if (ikl.eq.ikl0CM(1)) wsfre(kl) = qs_FRZ
	2956
	2957	! #hy END IF
	2958	! #hy END IF
	2959
	2960	! Debug
	2961	! ~~~~~
	2962	! #wH debugH( 1:35) = 'Lin et al.(1983): Rain Freezing '
	2963	! #wH debugH(36:70) = ' '
	2964	! #wH proc_1 = 'Qsfr g/kg'
	2965	! #wH procv1 = qs_FRZ
	2966	! #wH proc_2 = ' '
	2967	! #wH procv2 = 0.
	2968	! #wH proc_3 = ' '
	2969	! #wH procv3 = 0.
	2970	! #wH proc_4 = ' '
	2971	! #wH procv4 = 0.
	2972	! #wh include 'CMiPhy_Debug.h'
	2973	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1)) &
	2974	! #wH& debugV(k,16) = qs_FRZ
	2975
	2976	END DO
	2977
	2978
	2979
	2980
	2981	! Debug (Summary)
	2982	! ===============
	2983
	2984	! #wH DO k=mz1_CM,mzp
	2985
	2986	! #wH IF (ii__AP(ikl).EQ.i0__CM(1).AND.jj__AP(ikl).EQ.j0__CM(1))THEN
	2987	! #wH IF (k .EQ.mz1_CM) THEN
	2988	! #wH write(6,6022)
	2989	6022 format(/,'CMiPhy STATISTICS' &
	2990	& /,'=================')
	2991	! #wH write(6,6026)
	2992	6026 format( ' T_Air Qv Qw g/kg Qi g/kg CLOUDS % ' &
	2993	& , ' Qs g/kg Qr g/kg' &
	2994	& ,' Qi+ E.K.' &
	2995	& ,' Qi+ Mey.' &
	2996	& ,' Qi- Sub.' &
	2997	& ,' Qi- Mlt.' &
	2998	& ,' Qw+ Cds.' &
	2999	& ,' Qraut r+' &
	3000	& ,' QsAUT s+' &
	3001	& ,' Qracw r+')
	3002	! #wH END IF
	3003	! #wH write(6,6023) k &
	3004	! #wH& , Ta__CM(ikl,k)-Tf_Sno &
	3005	! #wH& ,1.e3* qv__DY(ikl,k) &
	3006	! #wH& ,1.e3* qw__CM(ikl,k) &
	3007	! #wH& ,1.e3* qi__CM(ikl,k) &
	3008	! #wH& ,1.e2* CFraCM(ikl,k) &
	3009	! #wH& ,1.e3* qs__CM(ikl,k) &
	3010	! #wH& ,1.e3* qr__CM(ikl,k) &
	3011	! #wH& ,(1.e3* debugV(k,kv),kv=1,08)
	3012	6023 format(i3,f6.1,f5.2,2f9.6,f9.1,2f9.3,8f9.6)
	3013	! #wH IF (k .EQ.mzp ) THEN
	3014	! #wH write(6,6026)
	3015	! #wH write(6,*) ' '
	3016	! #wH write(6,6024)
	3017	6024 format( 8x,'Z [km]' &
	3018	& ,' RH.w.[%]' &
	3019	& ,' RH.i.[%]' ,9x &
	3020	& ,' Vss cm/s' &
	3021	& ,' Vrr cm/s' &
	3022	& ,' Qsacw s+' &
	3023	& ,' Qsaci s+' &
	3024	& ,' Qiacr r+' &
	3025	& ,' Qracs ds' &
	3026	& ,' Qrevp w-' &
	3027	& ,' Qssub s-' &
	3028	& ,' Qsmlt s-' &
	3029	& ,' Qsfr s+')
	3030	! #wH DO nl=mz1_CM,mzp
	3031	! #wH write(6,6025) nl ,zsigma( nl)*1.e-3 &
	3032	! #wH& ,1.e2* qv__DY(ikl,nl)/qvswCM(ikl,nl) &
	3033	! #wH& ,1.e2* qv__DY(ikl,nl)/qvsiCM(ikl,nl) &
	3034	! #wH& ,1.e2* FallVs(ikl,nl) &
	3035	! #wH& ,1.e2* FallVr(ikl,nl) &
	3036	! #wH& ,(1.e3* debugV(nl,kv),kv=9,16)
	3037	6025 format(i3,f11.3, 2f9.1,9x, 2f9.1,8f9.6)
	3038	! #wH END DO
	3039	! #wH write(6,6024)
	3040	! #wH write(6,*) ' '
	3041	! #wH END IF
	3042
	3043	! #wH END IF
	3044
	3045	! #wH END DO
	3046
	3047
	3048
	3049
	3050	! Vertical Integrated Energy and Water Content
	3051	! ============================================
	3052
	3053	! #EW enr1EW(ikl) = 0.0d00
	3054	! #EW wat1EW(ikl) = 0.0d00
	3055
	3056	! #EW DO k=1,mzp
	3057	! #EW enr1EW(ikl) = enr1EW(ikl ) &
	3058	! #EW& +(Ta__CM(ikl,k) &
	3059	! #EW& -(qw__CM(ikl,k)+qr__CM(ikl,k)) *Lv_Cpd &
	3060	! #EW& -(qi__CM(ikl,k)+qs__CM(ikl,k)) Ls_Cpd)dsigmi(k)
	3061	! #EW wat1EW(ikl) = wat1EW(ikl ) &
	3062	! #EW& +(qv__DY(ikl,k) &
	3063	! #EW& + qw__CM(ikl,k)+qr__CM(ikl,k) &
	3064	! #EW& + qi__CM(ikl,k)+qs__CM(ikl,k) )*dsigmi(k)
	3065	! #EW END DO
	3066
	3067	! #ew enr1EW(ikl) = enr1EW(ikl ) * psa_DY(ikl) * Grav_I
	3068	! #EW wat1EW(ikl) = wat1EW(ikl ) * psa_DY(ikl) * Grav_I
	3069	! .. wat1EW [m] contains implicit factor 1.d3 [kPa-->Pa] /ro_Wat
	3070
	3071
	3072
	3073
	3074	! Precipitation
	3075	! =============
	3076
	3077	! Hydrometeors Fall Velocity
	3078	! --------------------------
	3079
	3080	! Pristine Ice Crystals Diameter and Fall Velocity
	3081	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	3082	DO k=mz1_CM,mzp
	3083
	3084	! #hy IF (qi__CM(ikl,k).gt.epsn) THEN
	3085	qi__OK = max(zer0,sign(un_1,qi__CM(ikl,k) - epsn))
	3086	! qi__OK = 1.0 if qi__CM(ikl,k) > epsn
	3087	! = 0.0 otherwise
	3088
	3089	! #hy IF (CCNiCM(ikl,k).gt.1.e0) THEN
	3090	CCNiOK = max(zer0,sign(un_1,CCNiCM(ikl,k) - 1.e0))
	3091	! CCNiOK = 1.0 if CCNiCM(ikl,k) > 1.e0
	3092	! = 0.0 otherwise
	3093
	3094	Flag_Fall_i = qi__OK * CCNiOK
	3095
	3096	Di_Pri = 0.16d0 exp(R_1by3log(R_1000roa_DY(ikl,k) &! Pristine Ice Crystals Diameter, where 6/(piro_I)**1/3 ~ 0.16
	3097	& *max(epsn,qi__CM(ikl,k))/max(un_1 ,CCNiCM(ikl,k)))) ! REF.: Levkov et al. 1992, Contr. Atm. Phys. 65, (5) p.37
	3098
	3099	FallVi(ikl,k) = Flag_Fall_i * 7.d2*Di_Pri &! Terminal Fall Velocity for Pristine Ice Crystals
	3100	& exp( 0.35d0 log(roa_DY(ikl,mzp) / roa_DY(ikl,k))) ! REF.: Levkov et al. 1992, Contr. Atm. Phys. 65, (4) p.37
	3101	! #hy ELSE
	3102	! #hy FallVi(ikl,k) = 0.0d00
	3103
	3104	! #hy END IF
	3105	! #hy END IF
	3106
	3107	END DO
	3108
	3109	! Set Up of the Numerical Scheme
	3110	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	3111	! #EW watfEW(ikl) = 0.d0 ! Water Flux (Atmosphere --> Surface)
	3112
	3113	! Snow and Rain Fall Velocity (Correction)
	3114	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	3115	DO k=mz1_CM,mzp
	3116	FallVi(ikl,k) = FallVi(ikl,k) *qi__CM(ikl,k)/max(qi__CM(ikl,k),epsn)
	3117	FallVs(ikl,k) = FallVs(ikl,k) *qs__CM(ikl,k)/max(qs__CM(ikl,k),epsn)
	3118	! #VW FallVw(ikl,k) = FallVw(ikl,k) *qw__CM(ikl,k)/max(qw__CM(ikl,k),epsn)
	3119	FallVr(ikl,k) = FallVr(ikl,k) *qr__CM(ikl,k)/max(qr__CM(ikl,k),epsn)
	3120	END DO
	3121
	3122
	3123	! -----------------------------------------------------
	3124	! Droplets Precipitation (Implicit Scheme)
	3125	! Pristine Ice Crystals Precipitation (Implicit Scheme)
	3126	! Snow Particles Precipitation (Implicit Scheme)
	3127	! Rain Drops Precipitation (Implicit Scheme)
	3128	! -----------------------------------------------------
	3129
	3130	qwLoss(mz1_CM-1) = 0.
	3131	qiLoss(mz1_CM-1) = 0.
	3132	qsLoss(mz1_CM-1) = 0.
	3133	qrLoss(mz1_CM-1) = 0.
	3134
	3135	! Precipitation Mass & Flux
	3136	! ~~~~~~~~~~~~~~~~~~~~~~~~~
	3137	DO k= mz1_CM,mzp
	3138	qwLoss(k) = 0.
	3139
	3140	a_rodz = Grav_I* psa_DY( ikl) *dsigmi(k) ! Air Mass
	3141	! #VW qwFlux = dt__CM* FallVw(ikl,k) *roa_DY(ikl,k) ! Flux Fact. (droplets)
	3142	qiFlux = dt__CM* FallVi(ikl,k) *roa_DY(ikl,k) ! Flux Fact. (crystals)
	3143	qsFlux = dt__CM* FallVs(ikl,k) *roa_DY(ikl,k) ! Flux Fact. (snow)
	3144	qrFlux = dt__CM* FallVr(ikl,k) *roa_DY(ikl,k) ! Flux Fact. (rain)
	3145
	3146	! #VW qwrodz = qw__CM(ikl,k) *a_rodz &! Droplets Mass
	3147	! #VW& + 0.5 *qwLoss(k-1) ! From abov.
	3148	qirodz = qi__CM(ikl,k) *a_rodz &! Crystals Mass
	3149	& + 0.5 *qiLoss(k-1) ! From abov.
	3150	qsrodz = qs__CM(ikl,k) *a_rodz &! Snow Mass
	3151	& + 0.5 *qsLoss(k-1) ! From abov.
	3152	qrrodz = qr__CM(ikl,k) *a_rodz &! Rain Mass
	3153	& + 0.5 *qrLoss(k-1) ! From abov.
	3154
	3155	! #VW wRatio = &! Var. Fact.
	3156	! #VW& min(2.,qwFlux /a_rodz ) ! Flux Limi.
	3157	iRatio = &! Var. Fact.
	3158	& min(2.,qiFlux /a_rodz ) ! Flux Limi.
	3159	sRatio = &! Var. Fact.
	3160	& min(2.,qsFlux /a_rodz ) ! Flux Limi.
	3161	rRatio = &! Var. Fact.
	3162	& min(2.,qrFlux /a_rodz ) ! Flux Limi.
	3163
	3164	! #VW qwLoss(k) = qwrodz *wRatio &! Mass Loss
	3165	! #VW& /(1.+wRatio *0.5) !
	3166	qiLoss(k) = qirodz *iRatio &! Mass Loss
	3167	& /(1.+iRatio *0.5) !
	3168	qsLoss(k) = qsrodz *sRatio &! Mass Loss
	3169	& /(1.+sRatio *0.5) !
	3170	qrLoss(k) = qrrodz *rRatio &! Mass Loss
	3171	& /(1.+rRatio *0.5) !
	3172
	3173	! #VW qwrodz = qwrodz -qwLoss(k) &!
	3174	! #VW& + 0.5 *qwLoss(k-1) ! From abov.
	3175	qirodz = qirodz -qiLoss(k) &!
	3176	& + 0.5 *qiLoss(k-1) ! From abov.
	3177	qsrodz = qsrodz -qsLoss(k) &!
	3178	& + 0.5 *qsLoss(k-1) ! From abov.
	3179	qrrodz = qrrodz -qrLoss(k) &!
	3180	& + 0.5 *qrLoss(k-1) ! From abov.
	3181
	3182	! Cooling from above precipitating flux
	3183	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	3184	Ta__CM(ikl,k) = &!
	3185	& (Ta__CM(ikl,k ) * a_rodz &!
	3186	& +Ta__CM(ikl,k-1) *(qwLoss(k-1)+qiLoss(k-1)+qsLoss(k-1)+qrLoss(k-1))) &!
	3187	& /(a_rodz +qwLoss(k-1)+qiLoss(k-1)+qsLoss(k-1)+qrLoss(k-1))
	3188
	3189	! #VW qw__CM(ikl,k) = qwrodz /a_rodz
	3190	qi__CM(ikl,k) = qirodz /a_rodz
	3191	qs__CM(ikl,k) = qsrodz /a_rodz
	3192	qr__CM(ikl,k) = qrrodz /a_rodz
	3193	ENDDO
	3194
	3195	! Precipitation reaching the Surface
	3196	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	3197	d_rain = qrLoss(mzp) + qwLoss(mzp) ! d_rain contains implicit factor 1.e3[kPa->Pa]/ro_Wat[kg/m2->m w.e.]
	3198	RainCM(ikl) = RainCM(ikl) + d_rain ! RainCM: rain precipitation height since start of run [m]
	3199	d_snow = qiLoss(mzp) ! d_snow contains implicit factor 1.e3[kPa->Pa]/ro_Wat[kg/m2->m w.e.]
	3200	Ice_CM(ikl) = Ice_CM(ikl) + d_snow ! Ice_CM: ice precipitation height since start of run [m]
	3201	d_snow = qsLoss(mzp) + qiLoss(mzp) ! d_snow contains implicit factor 1.e3[kPa->Pa]/ro_Wat[kg/m2->m w.e.]
	3202	SnowCM(ikl) = SnowCM(ikl) + d_snow ! SnowCM: ice + snow precipitation height since start of run [m]
	3203
	3204	! #EW watfEW(ikl ) = watfEW(ikl ) - d_rain - d_snow
	3205
	3206	d_rain = 0.0
	3207	d_snow = 0.0
	3208
	3209
	3210
	3211
	3212	! Fractional Cloudiness ! Guess may be computed (Ek&Mahrt91 fracSC=.T.)
	3213	! ====================== ! Final value computed below
	3214
	3215	! #sc IF (Frac__Clouds.AND..NOT.fracSC) THEN
	3216	IF (Frac__Clouds) THEN
	3217	IF(fraCEP) THEN ! ECMWF Large Scale Cloudiness
	3218	! ----------------------------
	3219	DO k=mz1_CM,mzp
	3220	CFraCM(ikl,k) = (qi__CM(ikl,k) + qw__CM(ikl,k)&
	3221	& +qs__CM(ikl,k) * 0.33 &
	3222	& * (1.-min(1.,exp((Ta__CM(ikl,k) -258.15)*0.1)))) &
	3223	& / (0.02 * qvswCM(ikl,k) )
	3224	CFraCM(ikl,k) =min(1.000 , CFraCM(ikl,k))
	3225	CFraCM(ikl,k) =max(0.001 , CFraCM(ikl,k)) &
	3226	& *max(zer0,sign(un_1,qi__CM(ikl,k) + qw__CM(ikl,k)&
	3227	& +qs__CM(ikl,k) -3.E-9 ))
	3228	END DO
	3229	ELSE ! XU and Randall 1996, JAS 21, p.3099 (4)
	3230	! ----------------------------
	3231	DO k=mz1_CM,mzp
	3232	qvs_wi= qvswCM(ikl,k)
	3233	! #wi qvs_wi=max(epsn,((qi__CM(ikl,k)+qs__CM(ikl,k))*qvsiCM(ikl,k) &
	3234	! #wi& +qw__CM(ikl,k) *qvswCM(ikl,k)) &
	3235	! #wi& /max(epsn,qi__CM(ikl,k)+qs__CM(ikl,k) +qw__CM(ikl,k)))
	3236	RHumid=min(RH_MAX,max(qv__DY(ikl,k),qv_MIN) / qvs_wi)
	3237	argEXP= ((RH_MAX -RHumid) * qvs_wi)**0.49
	3238	argEXP=min(100. *(qi__CM(ikl,k)+qw__CM(ikl,k) &
	3239	& +qs__CM(ikl,k) * 0.33 &
	3240	& * (1.-min(1.,exp((Ta__CM(ikl,k)-258.15)*0.1)))) &
	3241	& /max(epsn ,argEXP),ea_MAX)
	3242
	3243	CFraCM(ikl,k) = ( RHumid ** 0.25 ) * ( 1. - exp(-argEXP) )
	3244	END DO
	3245	END IF
	3246
	3247	ELSE
	3248	! #sc ELSE IF (.NOT.Frac__Clouds) THEN
	3249	DO k=mz1_CM,mzp
	3250	qCloud = qi__CM(ikl,k) + qw__CM(ikl,k)
	3251	! #hy IF (qCloud.gt.epsn) THEN
	3252	CFraCM(ikl,k) = max(zer0,sign(un_1, qCloud - epsn))
	3253	! CFraCM(ikl,k) = 1.0 if qCloud > epsn
	3254	! = 0.0 otherwise
	3255
	3256	! #hy END IF
	3257	END DO
	3258
	3259	END IF
	3260
	3261
	3262
	3263
	3264	! Vertically Integrated Energy and Water Content
	3265	! ==============================================
	3266
	3267	! #EW enr2EW(ikl) = 0.0d00
	3268	! #EW wat2EW(ikl) = 0.0d00
	3269	! .. Vertical Integrated Energy and Water Content
	3270
	3271	! #EW DO k=1,mzp
	3272	! #EW enr2EW(ikl) = enr2EW(ikl) &
	3273	! #EW& + (Ta__CM(ikl,k) &
	3274	! #EW& - (qw__CM(ikl,k)+qr__CM(ikl,k))*Lv_Cpd &
	3275	! #EW& - (qi__CM(ikl,k)+qs__CM(ikl,k))Ls_Cpd) dsigmi(k)
	3276	! #EW wat2EW(ikl) = wat2EW(ikl) &
	3277	! #EW& + (qv__DY(ikl,k) &
	3278	! #EW& + qw__CM(ikl,k)+qr__CM(ikl,k) &
	3279	! #EW& + qi__CM(ikl,k)+qs__CM(ikl,k) ) *dsigmi(k)
	3280	! #EW END DO
	3281
	3282	! #ew enr2EW(ikl) = enr2EW(ikl) * psa_DY(ikl) * Grav_I
	3283	! #EW wat2EW(ikl) = wat2EW(ikl) * psa_DY(ikl) * Grav_I
	3284	! .. wat2EW [m] contains implicit factor 1.d3 [kPa-->Pa] /ro_Wat
	3285
	3286
	3287
	3288
	3289	! Limits on Microphysical Variables
	3290	! =================================
	3291
	3292	DO k=1,mz1_CM
	3293	qv__DY(ikl,k)=max(qv__DY(ikl,k),qv_MIN)
	3294	qv__DY(ikl,k)=min(qv__DY(ikl,k),qvsiCM(ikl,k))
	3295	qw__CM(ikl,k)= zer0
	3296	qi__CM(ikl,k)= zer0
	3297	CCNiCM(ikl,k)= zer0
	3298	qr__CM(ikl,k)= zer0
	3299	qs__CM(ikl,k)= zer0
	3300	END DO
	3301
	3302	DO k=mz1_CM,mzp
	3303	qw__CM(ikl,k)=max(zer0, qw__CM(ikl,k))
	3304	qi__CM(ikl,k)=max(zer0, qi__CM(ikl,k))
	3305	CCNiCM(ikl,k)=max(zer0, CCNiCM(ikl,k))
	3306	qr__CM(ikl,k)=max(zer0, qr__CM(ikl,k))
	3307	qs__CM(ikl,k)=max(zer0, qs__CM(ikl,k))
	3308	END DO
	3309
	3310
	3311
	3312
	3313	! +++++++++++++++++++
	3314	ENDDO ! ikl=1,kcolp
	3315	! +++++++++++++++++++
	3316
	3317
	3318
	3319
	3320	! OUTPUT
	3321	! ======
	3322
	3323	! #WH IF (mod(MinuTU,6).eq.0.and.Sec_TU.eq.0.and.ikl0CM(1).gt.0) THEN
	3324	! #WH write(6,1030) HourTU,MinuTU,Sec_TU,it_EXP,i0__CM(1),j0__CM(1)
	3325	1030 format(//,i4,'UT',i2,'m',i2,'s (iter.',i6,') / Pt.(',2i4,')' &
	3326	& ,/,' ==========================================')
	3327	! #WH write(6,1031)(k, Z___DY(ikl0CM(1),k),qv__DY(ikl0CM(1),k), &
	3328	! #WH& 1.d3qi_io0(k), 1.d3qi__CM(ikl0CM(1),k), &
	3329	! #WH& 1.d3* wihm1(k),1.d3* wihm2(k), 1.d3* wicnd(k), &
	3330	! #WH& 1.d3* widep(k),1.d3* wisub(k), 1.d3* wimlt(k),k=mz1_CM,mzp)
	3331	1031 format(/, &
	3332	& ' \| Water Vapor \| Cloud Ice, Time n & n+1', &
	3333	& ' Cloud Ice Nucleation Processes \|', &
	3334	& ' Bergeron Sublimation Melting ', &
	3335	& /,' k z[m] \| qv [g/kg] \| qi_n [g/kg] qi_n+[g/kg]', &
	3336	& ' QiHm1[g/kg] QiHm2[g/kg] QiCnd[g/kg] \|', &
	3337	& ' QiDep[g/kg] QiSub[g/kg] QiMlt[q/kg]', &
	3338	& /,'------------+--------------+-------------------------', &
	3339	& '-------------------------------------+', &
	3340	& '-------------------------------------', &
	3341	& /,(i3,f8.1,' \| ',f12.6,' \| ',2f12.6,3d12.4,' \| ',3d12.4))
	3342
	3343	! #WH write(6,1032)(k,Z___DY(ikl0CM(1),k) &
	3344	! #WH& ,1.d3qs___0(ikl0CM(1),k),1.d3qs__CM(ikl0CM(1),k) &
	3345	! #WH& ,1.d3* wsaut(k), 1.d3* wsaci(k) &
	3346	! #WH& ,1.d3* wsacw(k), 1.d3* wiacr(k) &
	3347	! #WH& ,1.d3* wsacr(k), 1.d3* wssub(k) &
	3348	! #WH& , FallVs(k,ikl0CM(1)),k=mz1_CM,mzp)
	3349	1032 format(/, &
	3350	& ' \| Snow Flakes, Time n&n+1 Autoconver. \|', &
	3351	& ' Accretion Processes ===> Snow Flakes \|', &
	3352	& ' Sublimation \| Term.F.Vel', &
	3353	& /,' k z[m] \| qs_n [g/kg] qs_n+[g/kg] QsAUT[g/kg] \|', &
	3354	& ' Qsaci[g/kg] Qsacw[g/kg] Qiacr[g/kg] Qsacr[g/kg] \|', &
	3355	& ' QsSub[g/kg] \| vs [m/s]', &
	3356	& /,'------------+--------------------------------------+', &
	3357	& '--------------------------------------------------+', &
	3358	& '--------------+-----------', &
	3359	& /,(i3,f8.1,' \| ',2f12.6,e12.4,' \| ',4d12.4,' \| ',e12.4, &
	3360	& ' \| ',f10.6))
	3361
	3362	! #WH write(6,1033)(k,Z___DY(ikl0CM(1),k),Ta__CM(ikl0CM(1),k)
	3363	! #WH& ,1.d3qw_io0(k) ,1.d3qw__CM(ikl0CM(1),k)
	3364	! #WH& ,1.d3* wwevp(k) ,1.d2*CFraCM(ikl0CM(1),k),k=mz1_CM,mzp)
	3365	1033 format(/, &
	3366	& /,' \| Temperat.\| Cloud Water, Time n&n+1', &
	3367	& ' Condens/Evp \| Cloud ', &
	3368	& /,' k z[m] \| T [K] \| qw_n [g/kg] qw_n+[g/kg]', &
	3369	& ' QwEvp[g/kg] \| Fract.', &
	3370	& /,'------------+----------+-------------------------', &
	3371	& '-------------+-------', &
	3372	& /,(i3,f8.1,' \| ',f8.3,' \| ',2f12.6,e12.4,' \| ',f5.1))
	3373
	3374	! #WH write(6,1034)(k,Z___DY(ikl0CM(1),k), &
	3375	! #WH& ,1.d3qr___0(k,ikl0CM(1)),1.d3qr__CM(ikl0CM(1),k), &
	3376	! #WH& ,1.d3* wraut(k) ,1.d3* wracw(k) &
	3377	! #WH& ,1.d3* wraci(k) ,1.d3* wracs(k) &
	3378	! #WH& ,1.d3* wrevp(k) ,1.d3* wsfre(k) &
	3379	! #WH& , FallVr(k,ikl0CM(1)), k=mz1_CM,mzp)
	3380	1034 format(/, &
	3381	& /,' \| Rain Drops, Time n&n+1 Autoconver. \|', &
	3382	& ' Accretion Processes ===> Rain Drops \|', &
	3383	& ' Evaporation Freezing \| Term.F.Vel', &
	3384	& /,' k z[m] \| qr_n [g/kg] qr_n+[g/kg] Qraut[g/kg] \|', &
	3385	& ' Qracw[g/kg] Qraci[g/kg] Qracs[g/kg] \|', &
	3386	& ' QrEvp[g/kg] QsFre[g/kg] \| vr [m/s]', &
	3387	& /,'------------+--------------------------------------+', &
	3388	& '--------------------------------------+', &
	3389	& '--------------------------+-----------', &
	3390	& /,(i3,f8.1,' \| ',2f12.6,e12.4,' \| ',3d12.4,' \| ',2d12.4, &
	3391	& ' \| ',f10.6))
	3392
	3393	! #WH DO k=mz1_CM,mzp
	3394	! #WH wihm1(k) = 0.d0
	3395	! #WH wihm2(k) = 0.d0
	3396	! #WH wicnd(k) = 0.d0
	3397	! #WH widep(k) = 0.d0
	3398	! #WH wisub(k) = 0.d0
	3399	! #WH wimlt(k) = 0.d0
	3400	! #WH wwevp(k) = 0.d0
	3401	! #WH wraut(k) = 0.d0
	3402	! #WH wsaut(k) = 0.d0
	3403	! #WH wracw(k) = 0.d0
	3404	! #WH wsacw(k) = 0.d0
	3405	! #WH wsaci(k) = 0.d0
	3406	! #WH wraci(k) = 0.d0
	3407	! #WH wiacr(k) = 0.d0
	3408	! #WH wsacr(k) = 0.d0
	3409	! #WH wracs(k) = 0.d0
	3410	! #WH wrevp(k) = 0.d0
	3411	! #WH wssub(k) = 0.d0
	3412	! #WH wsmlt(k) = 0.d0
	3413	! #WH wsfre(k) = 0.d0
	3414	! #WH END DO
	3415	! #WH END IF
	3416
	3417
	3418	! Vertical Integrated Energy and Water Content: OUTPUT
	3419	! ====================================================
	3420
	3421	! #EW IF (ikl0CM(1).gt.0) THEN
	3422	! #EW WaterB = wat2EW(ikl0CM(1))-wat1EW(ikl0CM(1))-watfEW(ikl0CM(1))
	3423	! #EW write(6,606) it_EXP, &
	3424	! #EW& enr0EW(ikl0CM(1)),1.d3*wat0EW(ikl0CM(1)), &
	3425	! #EW& mphyEW(ikl0CM(1)), &
	3426	! #EW& enr1EW(ikl0CM(1)),1.d3*wat1EW(ikl0CM(1)), &
	3427	! #EW& enr2EW(ikl0CM(1)),1.d3*wat2EW(ikl0CM(1)), &
	3428	! #EW& 1.d3*watfEW(ikl0CM(1)), &
	3429	! #EW& 1.d3*WaterB
	3430	606 format(i9,' Before mPhy: E0 =',f12.6,' W0 = ',f9.6,3x,a20 &
	3431	& ,3x,/,9x,' Before Prec: E1 =',f12.6,' W1 = ',f9.6 &
	3432	& , /,9x,' After Prec: E2 =',f12.6,' W2 = ',f9.6 &
	3433	& , ' W Flux =',f9.6 &
	3434	& , ' Div(W) =',e9.3)
	3435	! #EW END IF
	3436
	3437	IF (MinuTU.eq.0.and.Sec_TU.eq.0.and.mod(HourTU,3).eq.0) THEN
	3438	DO ipt_CM=1,npt_CM
	3439	ikl=ikl0CM(ipt_CM)
	3440
	3441	write(4,1037) HourTU,MinuTU, &
	3442	& i0__CM(ipt_CM),j0__CM(ipt_CM)
	3443	1037 format(/,' Ice-Crystal mPhy ', &
	3444	& 2x,' ',2x,1x,i2,'h',i2,'UT', &
	3445	& ' -- Grid Point (',i5,',',i5,')', &
	3446	& /,' =========================================================='&
	3447	& , /,' \| z [m] \| T [K] \| qi[g/kg] \|' &
	3448	& , ' Ni [m-3] \| Ni0[m-3] \| vi [m/s] \| qs[g/kg] \|' &
	3449	& , /,'-----+---------+--------+----------+' &
	3450	& , '----------+----------+----------+----------+')
	3451	write(4,1038)(k,Z___DY(ikl,k),Ta__CM(ikl,k) &
	3452	& ,qi__CM(ikl,k)*1.d3 &
	3453	& ,CCNiCM(ikl,k),Fletch(ikl,k) &
	3454	& ,FallVi(ikl,k),qs__CM(ikl,k)*1.d3 &
	3455	& ,k=mz1_CM,mzp)
	3456	1038 format((i4,' \|' , f8.1,' \|',f7.2,' \|',f9.6,' \|', &
	3457	& 2(d9.3,' \|'),2(f9.6,' \|')))
	3458
	3459	END DO
	3460	END IF
	3461
	3462
	3463
	3464
	3465	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3466	! !
	3467	! DE-ALLOCATION !
	3468	! ============= !
	3469	! !
	3470	IF (FlagDALLOC) THEN !
	3471
	3472	deallocate ( qw_io0 ) ! Droplets Concentration entering CMiPhy [kg/kg]
	3473	deallocate ( qi_io0 ) ! Ice Part. Concentration entering CMiPhy [kg/kg]
	3474	deallocate ( qs___0 ) ! Snow Part. Concentration entering CMiPhy [kg/kg]
	3475	! #qg deallocate ( qg___0 ) ! Graupels Concentration entering CMiPhy [kg/kg]
	3476	deallocate ( qr___0 ) ! Rain Drops Concentration entering CMiPhy [kg/kg]
	3477	deallocate ( Ta_dgC ) ! Air Temperature [dgC]
	3478	deallocate ( sqrrro ) ! sqrt(roa(mzp)/roa(k)) [-]
	3479	deallocate ( qsiEFF ) ! EFFective Saturation Specific Humidity over Ice [kg/kg]
	3480	deallocate ( Fletch ) ! Monodisperse Nb of hexagonal Plates, Fletcher (1962) [-]
	3481	deallocate ( lamdaS ) ! Marshall-Palmer distribution parameter for Snow Particl.
	3482	! #qg deallocate ( lamdaG ) ! Marshall-Palmer distribution parameter for Graupels
	3483	deallocate ( lamdaR ) ! Marshall-Palmer distribution parameter for Rain Drops
	3484	deallocate ( ps_ACR ) ! Accretion of Snow by Rain Rate [kg/kg/s]
	3485	deallocate ( ps_ACW ) ! Accretion of Cloud Drop. by Snow Particl. Rate [kg/kg/s]
	3486	deallocate ( FallVw ) ! Sedimentation Velocity of Droplets
	3487	deallocate ( FallVi ) ! Sedimentation Velocity of Ice Particles
	3488	deallocate ( FallVs ) ! Sedimentation Velocity of Snow Particles
	3489	! #qg deallocate ( FallVg ) ! Sedimentation Velocity of Snow Particles
	3490	deallocate ( FallVr ) ! Sedimentation Velocity of Rain Drops
	3491	deallocate ( qwLoss ) ! Mass Loss related to Sedimentation of Rain Droplets
	3492	deallocate ( qiLoss ) ! Mass Loss related to Sedimentation of Ice Crystals
	3493	deallocate ( qsLoss ) ! Mass Loss related to Sedimentation of Snow Particles
	3494	deallocate ( qrLoss ) ! Mass Loss related to Sedimentation of Rain Drops
	3495	! #wH deallocate ( debugV ) ! Debug Variable (of 16 microphysical processes)
	3496
	3497	! #WH deallocate ( wihm1 ) ! Cloud Droplets Freezing
	3498	! #WH deallocate ( wihm2 ) ! Ice Crystals Homogeneous Sublimation
	3499	! #WH deallocate ( wicnd ) ! Ice Crystals Nucleation (Emde & Kahlig)
	3500	! #WH deallocate ( widep ) ! Ice Crystals Growth Bergeron Process (Emde & Kahlig)
	3501	! #WH deallocate ( wisub ) ! Ice Crystals Sublimation (Levkov)
	3502	! #WH deallocate ( wimlt ) ! Ice Crystals Melting
	3503	! #WH deallocate ( wwevp ) ! Water Vapor Condensation / Evaporation (Fractional Cloudiness)
	3504	! #WH deallocate ( wraut ) ! Cloud Droplets AUTO-Conversion
	3505	! #WH deallocate ( wsaut ) ! Ice Crystals AUTO-Conversion
	3506	! #WH deallocate ( wracw ) ! Accretion of Cloud Droplets by Rain, Ta > 0, --> Rain
	3507	! #WH deallocate ( wsacw ) ! Accretion of Cloud Droplets by Rain, Ta < 0, --> Snow
	3508	! #WH deallocate ( wsaci ) ! Accretion of Ice Crystals by Snow --> Snow
	3509	! #WH deallocate ( wraci ) ! Accretion of Ice Crystals by Rain --> Snow
	3510	! #WH deallocate ( wiacr ) ! Accretion of Rain by Ice Crystals --> Snow
	3511	! #WH deallocate ( wsacr ) ! Accretion of Rain by Snow --> Snow
	3512	! #WH deallocate ( wracs ) ! Accretion of Snow by Rain --> Snow, Rain
	3513	! #WH deallocate ( wrevp ) ! Rain Drops Evaporation
	3514	! #WH deallocate ( wssub ) ! Snow Particles Sublimation
	3515	! #WH deallocate ( wsmlt ) ! Snow Particles Melting
	3516	! #WH deallocate ( wsfre ) ! Rain Drops Freezing
	3517	! !
	3518	END IF !
	3519	! !
	3520	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3521
	3522
	3523
	3524
	3525	return
	3526	end subroutine CMiPhy

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format