LMDZ
CMiPhy.f90
Go to the documentation of this file.
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.1*qv__DY(ikl,k)/(RHcrit * qvswCM(ikl,k))
401 ! #wH proc_2 = 'R.Hum I[%]'
402 ! #wH procv2 = 0.1*qv__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-2*exp(-0.6*ta_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/(piNmbr*CCNwCM(ikl,k)*exp(dwTUR4))
579 ! #VW dwTUR5 = exp(R_5by3*log(dwTURi))
580 
581 ! #VW FallVw(ikl,k) = 1.19d8* piNmbr*CCNwCM(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.25*log((pinmbr*n0___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__ok*392. *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.e6*exp(-.12*min(0.,Ta_dgC(ikl,k))))
615 
616  lamdas(ikl,k) = exp(0.25*log((0.50*pinmbr*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__ok*2.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__ok*2976.*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__ok*20.06*sqrrro(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.250*log((piNmbr*n0___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.57*log(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.733e7*qsieff(ikl,k) &! 1.733e7=Ls*Ls*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.e6*exp(-.12*min(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.81*0.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+3*roa_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-2*roa_dy(ikl,k)*qsieff(ikl,k))
1187 ! 1.875d-5: Water Vapor Diffusivity in Air
1188  bdepvi = dt__cm *(1.-rh_ice)*4.0*di_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 = .622d0*es_hpa /(pa_hpa - .378d0*es_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.349d7*qsat_l/(tliqid*tliqid)) !
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, signdw*dwmesh) &!
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=Lv*Lv*0.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,signdw*dwmesh) &
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(dwmesh*dwmesh,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__ok*cfracm(ikl,k) *th_aut*4.09d6*pinmbr &!
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/(cc1*dwmesh+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_1by3*log(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/(pi*ro_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__d0*qs__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_1by3*log(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/(pi*ro_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_1by3*log(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.0*log( 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.5*c1saut*dt__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.8d0*log(lamdar(ikl,k))) ! where a = 842. and b = 0.8
2037  qr_acw = pr_acw*dt__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.e6*exp(-.12*min(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.25d0*log(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.99d0*log(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.41d0*log(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__cm*ps_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.e6*exp(-.12*min(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.25d0*log(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.99d0*log(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.41d0*log(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.8d0*log(lamdar(ikl,k)))
2362  qr_aci = pr_aci*dt__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.8d0*log(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_acr*dt__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.e6*exp(-.12 *min(0.,Ta_dgC(ikl,k))))
2483 
2484  pr_acs = 986.96d-3*(n0___r*n0___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___r*n0___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.0*lamdaS(ikl,k)/(lamdaG(ikl,k)*lamdaG(ikl,k)) &!
2589 ! #qg& +0.5*lamdaS(ikl,k)* lamdaS(ikl,k) &!
2590 ! #qg& /exp(3.0d0*log(lamdaG(ikl,k)))) !
2591 
2592 ! #cn n0___s = min(2.e8,2.e6*exp(-.12*min(0.,Ta_dgC(ikl,k))))
2593 
2594 ! #qg pgacs = 986.96d-3*(n0___g*n0___s/roa_DY(ikl,k)) &! 986.96: pi**2 * rhog
2595 ! #qg& * abs(FallVg(ikl,k)-FallVs(ikl,k))*flg*effACS !(graupel densitity assumed equal to snow density)
2596 ! #qg qgacs = pgacs*dt__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.0*lamdaS(ikl,k)/(lamdaG(ikl,k)*lamdaG(ikl,k)) &!
2633 ! #qg& +0.5*lamdaS(ikl,k)* lamdaS(ikl,k) &!
2634 ! #qg& /exp(3.0d0*log(lamdaG(ikl,k)))) !
2635 
2636 ! #cn n0___s = min(2.e8,2.e6*exp(-.12*min(0.,Ta_dgC(ikl,k))))
2637 
2638 ! #qg pgacr = 986.96d-3*(n0___g*n0___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___r*lr_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-2*roa_dy(ikl,k)*qsieff(ikl,k)) ! kT: Air Thermal Conductivity
2769 
2770 ! #cn n0___s = min(2.e8,2.e6*exp(-.12*min(0.,Ta_dgC(ikl,k))))
2771 
2772  ps_sub = 2*pinmbr*(1.d0-rh_ice)*n0___s*ls_num &!
2773  & /(1.d3*roa_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.e6*exp(-.12*min(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(-acoefm*dt__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_1by3*log(r_1000*roa_dy(ikl,k) &! Pristine Ice Crystals Diameter, where 6/(pi*ro_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.d3*qi_io0(k), 1.d3*qi__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.d3*qs___0(ikl0CM(1),k),1.d3*qs__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.d3*qw_io0(k) ,1.d3*qw__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.d3*qr___0(k,ikl0CM(1)),1.d3*qr__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, &
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
real(kind=real8), save r_1by3
real(kind=real8), save ea_min
real(kind=real8), save c_sund
real(kind=real8), dimension(31), save aa1
real(kind=real8), save rhcrit
real(kind=real8), dimension(:,:), allocatable, save qvswcm
real(kind=real8), save wa__hm
real(kind=real8), dimension(:,:), allocatable, save fallvr
real(kind=real8), save c2_ekm
real(kind=real8), dimension(:,:), allocatable, save fallvi
real(kind=real8), dimension(:,:), allocatable, save lamdas
real(kind=real8), save cc2
real(kind=real8), save un_1
real(kind=real8), dimension(:), allocatable, save qi_io0
real(kind=real8), dimension(:,:), allocatable, save qid_cm
real(kind=real8), dimension(:,:), allocatable, save kzh_at
real(kind=real8), dimension(:,:), allocatable, save qv__dy
real(kind=real8), dimension(:,:), allocatable, save qwd_cm
integer, dimension(npt_cm), save j0__cm
real(kind=real8), save qi0_dc
real(kind=real8), dimension(:,:), allocatable, save ta_dgc
real(kind=real8), dimension(:,:), allocatable, save qw__cm
real(kind=real8), dimension(:,:), allocatable, save ccnwcm
real(kind=real8), save tmaxhm
integer, parameter mz1_cm
subroutine cmiphy
Definition: CMiPhy.f90:2
real(kind=real8), save t_nuic
real(kind=real8), save c1_ekm
real(kind=real8), dimension(:,:), allocatable, save qr__cm
real(kind=real8), dimension(:), allocatable, save raincm
real(kind=real8), save qw_vol
real(kind=real8), save cfrmin
real(kind=real8), save pinmbr
real(kind=real8), dimension(:,:), allocatable, save ps_acr
real(kind=real8), dimension(:,:), allocatable, save cfracm
real(kind=real8), save b_nuic
real(kind=real8), save lc_cpd
real(kind=real8), save a_nuic
real(kind=real8), dimension(:), allocatable, save dsigmi
real(kind=real8), save a_nuid
real(kind=real8), save grav_i
real(kind=real8), dimension(:,:), allocatable, save ps_acw
real(kind=real8), save expwat
real(kind=real8), dimension(:,:), allocatable, save fletch
real(kind=real8), dimension(:), allocatable, save qrloss
real(kind=real8), dimension(:,:), allocatable, save qsieff
integer, save kcolp
real(kind=real8), save t_nuid
integer, dimension(npt_cm), save i0__cm
real(kind=real8), save qs__d0
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
real(kind=real8), save rh_max
real(kind=real8), dimension(:,:), allocatable, save wa__dy
real(kind=real8), save qwturb
real(kind=real8), dimension(:,:), allocatable, save qvsicm
real(kind=real8), save r_dair
integer, save ipt_cm
real(kind=real8), dimension(:,:), allocatable, save ccnicm
real(kind=real8), save n0___r
real(kind=real8), dimension(:), allocatable, save ice_cm
real(kind=real8), save tminhm
real(kind=real8), save rwcrit
real(kind=real8), dimension(:), allocatable, save qsloss
real(kind=real8), save watice
real(kind=real8), save di_hex
real(kind=real8), dimension(:,:), allocatable, save z___dy
real(kind=real8), save cc1
integer, parameter npt_cm
real(kind=real8), dimension(:,:), allocatable, save fallvs
!$Header!c c INCLUDE fxyprim h c c c Fonctions in line c c REAL fyprim REAL rj c c il faut la calculer avant d appeler ces fonctions c c c Fonctions a changer selon x(x) et y(y) choisis.c-----------------------------------------------------------------c c.....ici
real(kind=real8), save tqwfrz
real(kind=real8), dimension(:), allocatable, save qw_io0
real(kind=real8), save qv_min
real(kind=real8), dimension(:,:), allocatable, save roa_dy
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
real(kind=real8), dimension(:), allocatable, save snowcm
real(kind=real8), save r_1000
real(kind=real8), dimension(:), allocatable, save psa_dy
integer, save minutu
real(kind=real8), dimension(:,:), allocatable, save tke_at
real(kind=real8), save qismax
real(kind=real8), save epsq
integer, save mzp
real(kind=real8), save pt__dy
real(kind=real8), dimension(:), allocatable, save sigma
integer, save it_run
real(kind=real8), dimension(:,:), allocatable, save qi__cm
real(kind=real8), dimension(:), allocatable, save qiloss
integer, dimension(npt_cm), save ikl0cm
real(kind=real8), dimension(:,:), allocatable, save lamdar
real(kind=real8), dimension(:), allocatable, save qwloss
real(kind=real8), save tf_sno
real(kind=real8), save zer0
real(kind=real8), dimension(:,:), allocatable, save qs__cm
real(kind=real8), save b_nuid
real(kind=real8), dimension(:,:), allocatable, save ta__cm
real(kind=real8), save n0___s
real(kind=real8), save dd0
real(kind=real8), dimension(:,:), allocatable, save qr___0
real(kind=real8), dimension(:,:), allocatable, save qs___0
integer, save hourtu
real(kind=real8), dimension(31), save aa2
real(kind=real8), save qw_max
real(kind=real8), dimension(:,:), allocatable, save fallvw
integer, save sec_tu
real(kind=real8), save ssimax
real(kind=real8), save epsn
real(kind=real8), save ls_cpd
real(kind=real8), dimension(:,:), allocatable, save sqrrro
real(kind=real8), save dt__cm
real(kind=real8), save ea_max
real(kind=real8), save lv_cpd
real(kind=real8), save expwa2