Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

CMiPhy.f90 @ 5434

Last change on this file since 5434 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

Line
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