source: LMDZ5/branches/IPSLCM6.0.11.rc1/libf/phymar/CMiPhy.f90 @ 3152

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

Inclusion de la physique de MAR


Integration of MAR physics

File size: 175.1 KB
Line 
1      subroutine CMiPhy
2
3!------------------------------------------------------------------------------+
4!                                                         Mon 17-Jun-2013  MAR |
5!   MAR          CMiPhy                                                        |
6!     subroutine CMiPhy contains the MAR    Cloud Microphysical Scheme         |
7!                                                                              |
8!     version 3.p.4.1 created by H. Gallee,               Thu 21-Mar-2013      |
9!           Last Modification by H. Gallee,               Mon 17-Jun-2013      |
10!                                                                              |
11!------------------------------------------------------------------------------+
12!                                                                              |
13!   INPUT / OUTPUT: qv__DY(kcolp,mzp) : air   specific humidity        [kg/kg] |
14!   ^^^^^^^^^^^^^^  qw__CM(kcolp,mzp) : cloud drops                    [kg/kg] |
15!                   qi__CM(kcolp,mzp) : ice   crystals Concentration   [kg/kg] |
16!                   qs__CM(kcolp,mzp) : Snow  Particl. Concentration   [kg/kg] |
17!                   qr__CM(kcolp,mzp) : Rain  Drops    Concentration   [kg/kg] |
18!                                                                              |
19!   (to be added:   qg__CM(kcolp,mzp) : Graupels       Concentration   [kg/kg])|
20!                                                                              |
21!                   CCNwCM(kcolp,mzp) : cloud droplets number          [Nb/m3] |
22!                   CCNiCM(kcolp,mzp) : ice   crystals number          [Nb/m3] |
23!                                                                              |
24!                   CFraCM(kcolp,mzp) : cloud fraction                    [-]  |
25!                                                                              |
26!                   RainCM(kcolp    ) : rain  Precipitation           [m w.e.] |
27!                   SnowCM(kcolp    ) : snow  Precipitation           [m w.e.] |
28!                   Ice_CM(kcolp    ) : ice   Precipitation           [m w.e.] |
29!                                                                              |
30!                   qid_CM(kcolp,mzp) : Ice    Water Formation         [kg/kg] |
31!                   qwd_CM(kcolp,mzp) : Liquid Water Formation         [kg/kg] |
32!                                                                              |
33!   INPUT :         wa__DY(kcolp,mzp) : Vertical Wind Speed              [m/s] |
34!   ^^^^^           roa_DY(kcolp,mzp) : Air Density                     [T/m3] |
35!                   qvsiCM(kcolp,mzp) : Satur.specific humidity (ICE)  [kg/kg] |
36!                   qvswCM(kcolp,mzp) : Satur.specific humidity (LIQ.) [kg/kg] |
37!                                                                              |
38!   OUTPUT:                                                                    |
39!   ^^^^^^                                                                     |
40!                                                                              |
41!                                                                              |
42!   REFER. : 1) Ntezimana, unpubl.thes.LLN,          115 pp,     1993          |
43!   ^^^^^    2) Lin et al.       JCAM            22, 1065--1092, 1983          |
44!               (very similar, except that graupels are represented)           |
45!            3) Emde and Kahlig, Annal.Geophys.   7,  405-- 414, 1989          |
46!            4) Levkov et al.,   Contr.Atm.Phys. 65,   35--  57, 1992          |
47!            5) Meyers et al.,   JAM             31,  708-- 731, 1992          |
48!               (Primary Ice-Nucleation Parameterization)                      |
49!            6) Delobbe and Gallee, BLM          89,   75-- 107  1998          |
50!               (Partial Condensation Scheme)                                  |
51!                                                                              |
52! # OPTIONS: #hy  Additional IF/THEN/ELSE added precluding vectorization       |
53! # ^^^^^^^  #qg  Graupel Conservation Equation         (to  verify & include) |
54! #          #cn  Intercept Parameter / snow gamma distribution = fct(Ta)      |
55! #          #kk  Limitation of SCu fraction                                   |
56! #          #VW  Cloud Drop. Sediment.(Duynkerke .. 1995, JAS  52, p.2763     |
57! #          #pp  Emde & Kahlig Ice Crystal Deposition  (not included)         |
58! #          #wi  QSat modified by  qw/qi   Ratio       (not included)         |
59!                                                                              |
60! # DEBUG:   #WH  Additional Output (Each Process  is detailled)               |
61! # ^^^^^    #WQ  FULL       Output (Each Process  is detailled)               |
62! #          #wH  Additional Output                                            |
63! #          #wh  Additional Output (Include write CMiPhy_Debug.h)             |
64! #          #EW  Additional Output (Energy and Water Conservation)            |
65! #          #ew  Additional Output (Energy and Water Conservation)            |
66!                                                                              |
67!   REMARK : the sign '~' indicates that reference must be verified            |
68!   ^^^^^^^^                                                                   |
69!   CAUTION:     Partial Condensation Scheme NOT validated                     |
70!   ^^^^^^^      for SCu -- Cu Transition                                      |
71!                erf fonction is erroneous on HP                               |
72!                                                                              |
73!------------------------------------------------------------------------------+
74
75
76
77!  Global Variables
78!  ================
79
80      use Mod_Real
81      use Mod_PHY____dat
82      use Mod_PHY____grd
83      use Mod_PHY_CM_dat
84      use Mod_PHY_CM_grd
85      use Mod_PHY_CM_kkl
86      use Mod_PHY_DY_kkl
87      use Mod_PHY_AT_kkl
88
89
90
91!  Local  Variables
92!  ================
93
94      use Mod_CMiPhy_loc
95
96
97      IMPLICIT NONE
98
99
100      logical  ::  Heter_Freezng = .FALSE.  !  .TRUE. => Levkov et al. (1992)    Heterogeneous  Freezing of Cloud Droplets
101      logical  ::  Homog_Sublima = .FALSE.  !  .TRUE. => Levkov et al. (1992)    Homogeneous Sublimation of Cloud Ice Particles
102      logical  ::  Meyers        = .TRUE.   !  .TRUE. => Meyers et al. (1992)    Ice Nucleation I
103      logical  ::  AUTO_w_Sundqv = .TRUE.   !  .TRUE. => Sundqvist     (1988)    Autoconversion Cloud Droplets --> Rain
104      logical  ::  AUTO_w_LiouOu = .FALSE.  !  .TRUE. => Liou and Ou   (1989)    Autoconversion Cloud Droplets --> Rain (Tropical SCu ONLY)
105      logical  ::  AUTO_w_LinAll = .FALSE.  !  .TRUE. => Lin    et al. (1983)    Autoconversion Cloud Droplets --> Rain
106      logical  ::  AUTO_i_Levkov = .TRUE.   !  .TRUE. => Levkov et al. (1992)    Autoconversion Cloud Ice      --> Snow
107      logical  ::  AUTO_i_LevkXX = .TRUE.   !  .TRUE. => Levkov et al. (1992)    Autoconversion Cloud Ice      --> Snow (Deposition.Growth)
108      logical  ::  AUTO_i_EmdeKa = .FALSE.  !  .TRUE. => Emde & Kahlig (1989)    Autoconversion Cloud Ice      --> Snow
109      logical  ::  AUTO_i_Sundqv = .FALSE.  !  .TRUE. => Sundqvist               Autoconversion Cloud Ice      --> Snow
110                                            ! .FALSE. => Emde & Kahlig (1989)    Bergeron-Findeisen Process
111      logical  ::  fracSC        = .FALSE.  !  .TRUE. => Delobbe                 SCu Fractional Cloudiness may be set up if Frac__Clouds = .TRUE.
112      logical  ::  fraCEP        = .FALSE.  !  .TRUE. => ECMWF                   SCu Fractional Cloudiness
113      logical  ::  HalMos        = .TRUE.   !  .TRUE. => Levkov et al. (1992)    Ice Nucleation II (Hallet-Mossop Process)
114      logical  ::  graupel_shape = .TRUE.   !  .TRUE. => Snow Particles Shape:   Graupellike Snow Flakes of Hexagonal Type
115      logical  ::  planes__shape = .FALSE.  !  .TRUE. => Snow Particles Shape:   Unrimed Side Planes
116      logical  ::  aggrega_shape = .FALSE.  !  .TRUE. => Snow Particles Shape:   Aggregates of unrimed radiating assemblages
117
118      logical  ::  NO_Vec        = .TRUE.   !  .TRUE. => Preference of IF/THEN/ELSE to sign and max/min Functions
119
120      real(kind=real8)                             :: rad_ww            ! Cloud Droplets Radius                                     [...]
121      real(kind=real8)                             :: qvs_wi            ! Saturation Specific Humididy over Liquid Water          [kg/kg]
122                                                                        ! Ref.: Emde & Kahlig 1989, Ann.Geophys.    7, p.407 (5)
123      real(kind=real8)                             :: BNUCVI            ! Nucleation  I: Deposition & Condensation-Freez. Nucl.   [kg/kg]
124      real(kind=real8)                             :: SSat_I            ! Nucleation  I: Sursaturation % ICE                          [%]
125
126      real(kind=real8)                             :: Flag_T_NuId       ! Flag: T < T_NuId  => Nucleation I (Dep./Cond.)  may occur   [-]
127      real(kind=real8)                             :: CCNiId            !                                                         [Nb/m3]
128
129      real(kind=real8)                             :: Flag___NuIc       ! Flag: 1           => Nucleation I (Cont.Frz.)   may occur   [-]
130      real(kind=real8)                             :: Flag_T_NuIc       ! Flag: T < T_NuIc  => Nucleation I (Cont.Frz.)   may occur   [-]
131      real(kind=real8)                             :: qw__OK            ! Flag: 1           => Non-zero cloud droplets Concentration
132      real(kind=real8)                             :: qi__OK            ! Flag: 1           => Non-zero cloud ice      Concentration
133      real(kind=real8)                             :: CCNiOK            ! Flag: 1           => Non-zero cloud ice      Particles     
134      real(kind=real8)                             :: CCNiIc            !                                                         [Nb/m3]
135
136      real(kind=real8)                             :: Flag_TmaxHM       ! Flag: T < TmaxHM  => Nucleation II (Hall-Mossop) may occur  [-]
137      real(kind=real8)                             :: Flag_TminHM       ! Flag: T > TminHM  => Nucleation II (Hall-Mossop) may occur  [-]
138      real(kind=real8)                             :: Flag_wa__HM       ! Flag: w > wa__HM  => Nucleation II (Hall-Mossop) may occur  [-]
139      real(kind=real8)                             :: SplinJ            ! Hallet-Mossop Theory (Levkov et al., 1992,
140      real(kind=real8)                             :: SplinP            !                       Contr.Atm.Phy. 65, p.40)
141
142      real(kind=real8)                             :: Flag_Ta_Neg       ! Flag: T < 0.0dgC                                            [-]
143      real(kind=real8)                             :: Flag_TqwFrz       ! Flag: T < Temper. => Instantaneous Freezing may occur       [-]
144
145      real(kind=real8)                             :: RHumid            ! Relative Humidity                                           [-]
146
147      real(kind=real8)                             :: qi_Nu1,qi_Nu2     ! Nucleation ( 0 > Ta > -35 dgC ), Generations 1 and 2    [kg/kg]
148      real(kind=real8)                             :: qi_Nuc            ! Nucleation ( 0 > Ta > -35 dgC ), Generation (effective) [kg/kg]
149      real(kind=real8)                             :: NuIdOK            ! Nucleation  I:              Contact     -Freez. Nucl.       [-]
150      real(kind=real8)                             :: BSPRWI            ! Nucleation  I:              Contact     -Freez. Nucl.   [kg/kg]
151      real(kind=real8)                             :: BHAMWI            ! Nucleation II: Hallett-Mossop Ice-Multiplic. Process    [kg/kg]
152      real(kind=real8)                             :: BNUFWI            ! Heterogeneous Freezing of Cloud Droplets                [kg/kg]
153      real(kind=real8)                             :: BDEPVI            ! Ice   Crystals Sublim.                                  [kg/kg]
154      real(kind=real8)                             :: Flag_SURSat       ! Flag = 1/0 if (SUR/sub)Saturation                  FLAG     [-]
155      real(kind=real8)                             :: Flag_SUBSat       ! Flag = 1/0 if (SUB/sur)Saturation                  FLAG     [-]
156      real(kind=real8)                             :: Flag_Sublim       ! Flag = 1/0 for Sublimation Occurence or not                 [-]
157      real(kind=real8)                             :: RH_Ice            ! Relative Humidity   vs   Ice                                [-]
158      real(kind=real8)                             :: RH_Liq            ! Relative Humidity   vs   Water                              [-]
159      real(kind=real8)                             :: DenDi1,DenDi2     ! dqi/dt|sublimation:  terms of the denominator             [...]
160      real(kind=real8)                             :: dqsiqv            ! SUBsaturation  qsi - qv                                 [kg/kg]
161      real(kind=real8)                             :: dqvSUB            ! SURsaturation  qv  - qsiEFF                             [kg/kg]
162      real(kind=real8)                             :: dqvDUM            ! Water Vapor    Variation, dummy variable                [kg/kg]
163      real(kind=real8)                             :: dqiDUM            ! Ice   Crystals Variation, dummy variable                [kg/kg]
164
165      real(kind=real8)                             :: qCloud            ! qw + qi                                                 [kg/kg]
166      real(kind=real8)                             :: coefC2            ! Coefficient of Ek and Mahrt parameter C2_EkM            [../..]
167      real(kind=real8)                             :: pa_hPa,es_hPa     ! Pressure,     Pressure of Vapor at Saturat., over Water [hPa]
168      real(kind=real8)                             :: Qsat_L            ! Specific Concentration of Vapor at Saturat., over Water [kg/kg]
169                                                                        !       (even for temperatures smaller than freezing pt)
170      real(kind=real8)                             :: t_qvqw            ! qv+qw mixing ratio used in Partial Condensation Scheme  [kg/kg]
171      real(kind=real8)                             :: d_qvqw            ! qv+qw mixing ratio variation                            [kg/kg]
172      real(kind=real8)                             :: Kdqvqw            ! qv+qw vertical turbulent Flux                           [kg/kg m/s]
173      real(kind=real8)                             :: ww_TKE            !       vertical Velocity  Variance                       [m2/s2]
174      real(kind=real8)                             :: RH_TKE            !       Relative Humidity  Variance                       [../..]
175      real(kind=real8)                             :: qt_TKE            ! qv+qw                    Variance                       [../..]
176      real(kind=real8)                             :: TLiqid            ! Liquid Temperature                                          [K]
177      real(kind=real8)                             :: CFr_t1,CFr_t2     !                                                             [-]
178      real(kind=real8)                             :: CFrCoe            !                                                             [-]
179      real(kind=real8)                             :: CFraOK            !                                                             [-]
180      real(kind=real8)                             :: qwCFra            ! CloudAveraged Liquid Water Mixing Ratio                 [kg/kg]
181      real(kind=real8)                             :: qwMesh            ! Mesh Averaged Liquid Water Mixing Ratio                 [kg/kg]
182      real(kind=real8)                             :: dwMesh            ! Mesh Averaged Liquid Water Mixing Ratio Variation       [kg/kg]
183      real(kind=real8)                             :: signdw            ! Sign       of Liquid Water Mixing Ratio Variation  (-1/1)   [-]
184      real(kind=real8)                             :: Flag_dqwPos       ! Flag = 1/0 if Liquid Water Mixing Ratio Variation    >/< 0  [-]
185      real(kind=real8)                             :: updatw            !                                                             [-]
186      real(kind=real8)                             :: SCuLim            ! Fraction      Limit                                         [-]
187      real(kind=real8)                             :: ARGerf            ! Argument of erf function (used in partial condensation Scheme)
188      real(kind=real8)                             :: OUTerf            ! OUTPUT   of erf function
189      real(kind=real8)                             ::    erf            !             erf function
190      real(kind=real8)                             :: dwTUR4,dwTURi     !
191      real(kind=real8)                             :: dwTUR3,dwTUR2     !
192      real(kind=real8)                             :: dwTUR8,dwTURc     !
193      real(kind=real8)                             :: dwTUR5,dwTUR1     !
194
195      real(kind=real8)                             :: Di_Pri            ! Pristine Ice Diameter
196      real(kind=real8)                             :: c1saut,cnsaut     ! Ice   Crystals AUToconv.  Parameters
197      real(kind=real8)                             :: dtsaut            ! Ice   Crystals AUToconv.  Time Scale
198      real(kind=real8)                             :: ps_AUT            ! Ice   Crystals AUToconv. (BDEPIS, BAGRIS,...)      Rate [kg/kg/s]
199      real(kind=real8)                             :: qs_AUT            ! Ice   Crystals AUToconv. (BDEPIS, BAGRIS,...)           [kg/kg]
200
201      real(kind=real8)                             :: Flag_Ta_Pos       ! Flag = 1/0 if  Ta >/< 273.15 K                     FLAG     [-]
202      real(kind=real8)                             :: Flag_qiMELT       ! Flag = 1/0 for qi Melting or not                   FLAG     [-]
203      real(kind=real8)                             :: qxMelt            ! Potential qi to Melt                                    [kg/kg]
204      real(kind=real8)                             :: qiMELT            ! Effective qi to Melt                                    [kg/kg]
205      real(kind=real8)                             :: CiMelt            ! Effective CCNi removed   by qi   Melting                [nb/m3]
206
207      real(kind=real8)                             :: Flag_qsMELT       ! Flag = 1/0 for qs Melting or not                   FLAG     [-]
208      real(kind=real8)                             :: xCoefM            !
209      real(kind=real8)                             :: AcoefM,BcoefM     !
210      real(kind=real8)                             :: dTMELT            ! Temperature Offset before   Snow Particles Melting          [K]
211      real(kind=real8)                             :: qsMELT            ! Melt                     of Snow Particles              [kg/kg]
212
213      real(kind=real8)                             :: qs_ACW            ! Accretion of Cloud Drop. by Snow Particl.               [kg/kg]
214      real(kind=real8)                             :: Flag_qs_ACW       ! Accretion of Cloud Drop. by Snow Particl.          FLAG     [-]
215      real(kind=real8)                             :: Flag_qs_ACI       ! Accretion of Cloud Ice   by Snow Particl.          FLAG     [-]
216      real(kind=real8)                             :: effACI            ! Accretion of Cloud Ice   by Snow Particl.    Efficiency     [-]
217      real(kind=real8)                             :: ps_ACI            ! Accretion of Cloud Ice   by Snow Particl.          Rate [kg/kg/s]
218      real(kind=real8)                             :: qs_ACI            ! Accretion of Cloud Ice   by Snow Particl.               [kg/kg]
219      real(kind=real8)                             :: CNsACI            ! Accretion of Cloud Ice   by Snow Particl. CCNi decrease [nb/m3]
220      real(kind=real8)                             :: Flag_qs_ACR       ! Accretion of Snow        by Rain                   FLAG     [-]
221      real(kind=real8)                             :: coeACR            ! Accretion of Snow        by Rain Sedimentation Coeffic.   [...]
222      real(kind=real8)                             :: qs_ACR            ! Accretion of Snow        by Rain                        [kg/kg]
223      real(kind=real8)                             :: qs_ACR_R          ! Accretion of Snow        by Rain          => Rain       [kg/kg]
224      real(kind=real8)                             :: qs_ACR_S          !                                           => Graupels   [kg/kg]
225      real(kind=real8)                             :: Flag_qr_ACS       ! Accretion of Rain        by Snow Particl. => Snow  FLAG     [-]
226      real(kind=real8)                             :: coeACS            ! Accretion of Rain        by Snow Sedimentation Coeffic.   [...]
227      real(kind=real8)                             :: pr_ACS            ! Accretion of Rain        by Snow Particl. => Snow  Rate [kg/kg/s]
228      real(kind=real8)                             :: qr_ACS            ! Accretion of Rain        by Snow Particl. => Snow       [kg/kg]
229      real(kind=real8)                             :: qr_ACS_S          ! Accretion of Rain        by Snow Particl. => Snow       [kg/kg]
230      real(kind=real8)                             :: ps_SUB            ! Deposition/Sublim.    on/of Snow Particl.          Rate [kg/kg/s]
231      real(kind=real8)                             :: qs_SUB            ! Deposition/Sublim.    on/of Snow Particl.               [kg/kg]
232      real(kind=real8)                             :: ls_NUM            ! Sublimation of Snow: NUMerator   of Snow Distribut.Coef.  [...]
233      real(kind=real8)                             :: ls_DEN            ! Sublimation of Snow: DENominator of Snow Distribut.Coef.  [...]
234
235      real(kind=real8)                             :: Flag_Freeze       ! Freezing                 of Rain                   FLAG     [-]   
236      real(kind=real8)                             :: ps_FRZ            ! Freezing                 of Rain                   Rate [kg/kg/s]
237      real(kind=real8)                             :: qs_FRZ            ! Freezing                 of Rain                        [kg/kg]
238
239      real(kind=real8)                             :: rwMEAN            ! Droplets Autoconversion:    Mean     Radius
240
241      real(kind=real8)                             :: th_AUT            ! Cloud Droplets AUToconv. Threshold                          [-]
242      real(kind=real8)                             :: pr_AUT            ! Cloud Droplets AUToconv.                           Rate [kg/kg/s]
243      real(kind=real8)                             :: qr_AUT            ! Cloud Droplets AUToconv.                                [kg/kg]
244      real(kind=real8)                             :: Flag_qr_ACW       ! Accretion of Cloud Drop. by Rain                   FLAG     [-]
245      real(kind=real8)                             :: pr_ACW            ! Accretion of Cloud Drop. by Rain                   Rate [kg/kg/s]
246      real(kind=real8)                             :: qr_ACW            ! Accretion of Cloud Drop. by Rain                        [kg/kg]
247      real(kind=real8)                             :: Flag_qr_ACI       ! Accretion of Cloud Drop. by Rain          => Rain  FLAG     [-]
248      real(kind=real8)                             :: pr_ACI            ! Accretion of Cloud Ice   by Rain          => Rain  Rate [kg/kg/s]
249      real(kind=real8)                             :: qr_ACI            ! Accretion of Cloud Ice   by Rain          => Rain       [kg/kg]
250      real(kind=real8)                             :: CNrACI            ! Accretion of Cloud Ice   by Rain          CCNi decrease [nb/m3]
251      real(kind=real8)                             :: pi_ACR            ! Accretion of Cloud Ice   by Rain          => Snow  Rate [kg/kg/s]
252      real(kind=real8)                             :: qi_ACR            ! Accretion of Cloud Ice   by Rain          => Snow       [kg/kg]
253      real(kind=real8)                             :: Flag_DryAir       ! 1 => RH_Liq < 1                                    FLAG     [-]
254      real(kind=real8)                             :: Flag_qr_EVP       ! Evaporation of Rain                                FLAG     [-]
255      real(kind=real8)                             :: lr_NUM            ! Evaporation of Rain: NUMerator   of rain Distribut.Coef.  [...]
256      real(kind=real8)                             :: lr_DEN            ! Evaporation of Rain: DENominator of rain Distribut.Coef.  [...]
257      real(kind=real8)                             :: pr_EVP            ! Evaporation of Rain                                Rate [kg/kg/s]
258      real(kind=real8)                             :: qr_EVP            ! Evaporation of Rain                                     [kg/kg]
259
260      real(kind=real8)                             :: effACS            ! Accretion of Snow        by Graupels         Efficiency     [-]
261
262      real(kind=real8)                             :: a_rodz            ! Air               Mass
263      real(kind=real8)                             :: qwFlux            ! qw  Sedimentation Flux          Coefficient
264      real(kind=real8)                             :: qwrodz            ! Droplets          Mass
265      real(kind=real8)                             :: wRatio            ! Droplets          Mass          Ratio                       [-]
266      real(kind=real8)                             :: qiFlux            ! qi  Sedimentation Flux          Coefficient
267      real(kind=real8)                             :: qirodz            ! Crystals          Mass
268      real(kind=real8)                             :: iRatio            ! Crystals          Mass          Ratio                       [-]
269      real(kind=real8)                             :: qsFlux            ! qs  Sedimentation Flux          Coefficient
270      real(kind=real8)                             :: qsrodz            ! Snow              Mass
271      real(kind=real8)                             :: sRatio            ! Snow              Mass          Ratio                       [-]
272      real(kind=real8)                             :: qrFlux            ! qr  Sedimentation Flux          Coefficient
273      real(kind=real8)                             :: qrrodz            ! Rain              Mass
274      real(kind=real8)                             :: rRatio            ! Rain              Mass          Ratio                       [-]
275
276      real(kind=real8)                             :: Vw_MAX            ! MAX Sedimentation Velocity   of Droplets
277      real(kind=real8)                             :: Flag_Fall_i       !     Sedimentation            of Ice                FLAG     [-]
278      real(kind=real8)                             :: Vi_MAX            ! MAX Sedimentation Velocity   of Ice  Particles
279      real(kind=real8)                             :: Vs_MAX,VsMMAX     ! MAX Sedimentation Velocity   of Snow Particles
280      real(kind=real8)                             :: Vs__OK            !     Sedimentation Velocity   of Snow Particles
281      real(kind=real8)                             :: Vr_MAX,VrMMAX     ! MAX Sedimentation Velocity   of Rain Drops
282      real(kind=real8)                             :: Vr__OK            !     Sedimentation Velocity   of Rain Drops
283      real(kind=real8)                             :: dtPrec            !     Sedimentation Time Step
284      real(kind=real8)                             :: d_Snow            !     Sedimented                  Snow Part. over 1 time step [m]
285      real(kind=real8)                             :: d_Rain            !     Sedimented                  Rain Drops over 1 time step [m]
286
287      real(kind=real8)                             :: SatiOK            !
288      real(kind=real8)                             :: FlagNu            ! 1 =>Nucleation for -35 dgC < Ta < 0 dgc
289
290      real(kind=real8)                             :: Qw0_OK            ! HYDROMETEORS INPUT STATUS
291      real(kind=real8)                             :: Qi0_OK            !
292      real(kind=real8)                             :: Qi0qOK            !
293      real(kind=real8)                             :: Ci0_OK            !
294      real(kind=real8)                             :: Ci0cOK            !
295      real(kind=real8)                             :: Qs0_OK            !
296      real(kind=real8)                             :: Qr0_OK            !
297
298      real(kind=real8)                             :: qiBerg            ! Bergeron-Findeisen Process: avaiklable qi
299      real(kind=real8)                             :: qwBerg            ! Bergeron-Findeisen Process: avaiklable qw
300      real(kind=real8)                             :: qxBerg            ! Bergeron-Findeisen Process: potential qi Generation
301      real(kind=real8)                             :: a1Berg,a2Berg     ! Bergeron-Findeisen Process: parameters
302      real(kind=real8)                             :: a0Berg,afBerg     ! Bergeron-Findeisen Process: integration constants
303
304      real(kind=real8)                             :: WaterB            ! Vertically Integrated Water Budget
305
306      real(kind=real8)                             :: argEXP            !
307
308      integer                                      :: k     ,ikl        !
309      integer                                      :: i_Berg            !
310      integer                                      :: nItMAX,itFall     !
311      integer                                      :: ikl_io,io__Pt     !
312
313
314!  Debug Variables
315!  ---------------
316
317! #wH character(len=70)                            :: debugH            !
318! #wH character(len=10)                            :: proc_1,proc_2     !
319! #wH character(len=10)                            :: proc_3,proc_4     !
320! #wH real(kind=real8)                             :: procv1,procv2     !
321! #wH real(kind=real8)                             :: procv3,procv4     !
322! #wH integer                                      :: kv    ,nl         !
323
324
325
326
327!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
328!                                                                       !
329! ALLOCATION                                                            !
330! ==========                                                            !
331
332      IF (it_RUN.EQ.1 .OR. FlagDALLOC)                             THEN !
333
334          allocate ( qw_io0(      mzp) )                                ! Droplets   Concentration entering CMiPhy                  [kg/kg]
335          allocate ( qi_io0(      mzp) )                                ! Ice  Part. Concentration entering CMiPhy                  [kg/kg]
336          allocate ( qs___0(kcolp,mzp) )                                ! Snow Part. Concentration entering CMiPhy                  [kg/kg]
337! #qg     allocate ( qg___0(kcolp,mzp) )                                ! Graupels   Concentration entering CMiPhy                  [kg/kg]
338          allocate ( qr___0(kcolp,mzp) )                                ! Rain Drops Concentration entering CMiPhy                  [kg/kg]
339          allocate ( Ta_dgC(kcolp,mzp) )                                ! Air   Temperature                                           [dgC]
340          allocate ( sqrrro(kcolp,mzp) )                                ! sqrt(roa(mzp)/roa(k))                                         [-]
341          allocate ( qsiEFF(kcolp,mzp) )                                ! EFFective Saturation Specific Humidity over Ice           [kg/kg]
342          allocate ( Fletch(kcolp,mzp) )                                ! Monodisperse Nb of hexagonal Plates, Fletcher (1962)          [-]
343          allocate ( lamdaS(kcolp,mzp) )                                ! Marshall-Palmer distribution parameter for Snow Particl.
344! #qg     allocate ( lamdaG(kcolp,mzp) )                                ! Marshall-Palmer distribution parameter for Graupels
345          allocate ( lamdaR(kcolp,mzp) )                                ! Marshall-Palmer distribution parameter for Rain Drops
346          allocate ( ps_ACR(kcolp,mzp) )                                ! Accretion of Snow        by Rain                   Rate [kg/kg/s]
347          allocate ( ps_ACW(kcolp,mzp) )                                ! Accretion of Cloud Drop. by Snow Particl.          Rate [kg/kg/s]
348          allocate ( FallVw(kcolp,mzp) )                                ! Sedimentation Velocity   of Droplets
349          allocate ( FallVi(kcolp,mzp) )                                ! Sedimentation Velocity   of Ice  Particles
350          allocate ( FallVs(kcolp,mzp) )                                ! Sedimentation Velocity   of Snow Particles
351! #qg     allocate ( FallVg(kcolp,mzp) )                                ! Sedimentation Velocity   of Snow Particles
352          allocate ( FallVr(kcolp,mzp) )                                ! Sedimentation Velocity   of Rain Drops
353          allocate ( qwLoss(      mzp) )                                ! Mass Loss related to Sedimentation of Rain Droplets
354          allocate ( qiLoss(      mzp) )                                ! Mass Loss related to Sedimentation of Ice  Crystals
355          allocate ( qsLoss(      mzp) )                                ! Mass Loss related to Sedimentation of Snow Particles
356          allocate ( qrLoss(      mzp) )                                ! Mass Loss related to Sedimentation of Rain Drops
357! #wH     allocate ( debugV(      mzp,16) )                             ! Debug Variable (of 16 microphysical processes)
358
359! #WH     allocate ( wihm1(mzp) )                                       ! Cloud Droplets Freezing
360! #WH     allocate ( wihm2(mzp) )                                       ! Ice   Crystals Homogeneous Sublimation
361! #WH     allocate ( wicnd(mzp) )                                       ! Ice   Crystals Nucleation              (Emde & Kahlig)
362! #WH     allocate ( widep(mzp) )                                       ! Ice   Crystals Growth Bergeron Process (Emde & Kahlig)
363! #WH     allocate ( wisub(mzp) )                                       ! Ice   Crystals             Sublimation (Levkov)
364! #WH     allocate ( wimlt(mzp) )                                       ! Ice   Crystals Melting 
365! #WH     allocate ( wwevp(mzp) )                                       ! Water Vapor Condensation / Evaporation (Fractional Cloudiness)
366! #WH     allocate ( wraut(mzp) )                                       ! Cloud Droplets AUTO-Conversion
367! #WH     allocate ( wsaut(mzp) )                                       ! Ice   Crystals AUTO-Conversion
368! #WH     allocate ( wracw(mzp) )                                       ! Accretion of Cloud Droplets by Rain, Ta > 0, --> Rain
369! #WH     allocate ( wsacw(mzp) )                                       ! Accretion of Cloud Droplets by Rain, Ta < 0, --> Snow
370! #WH     allocate ( wsaci(mzp) )                                       ! Accretion of Ice   Crystals by Snow          --> Snow
371! #WH     allocate ( wraci(mzp) )                                       ! Accretion of Ice   Crystals by Rain          --> Snow
372! #WH     allocate ( wiacr(mzp) )                                       ! Accretion of Rain by Ice   Crystals          --> Snow
373! #WH     allocate ( wsacr(mzp) )                                       ! Accretion of Rain by Snow                    --> Snow
374! #WH     allocate ( wracs(mzp) )                                       ! Accretion of Snow by Rain                    --> Snow, Rain
375! #WH     allocate ( wrevp(mzp) )                                       ! Rain  Drops     Evaporation 
376! #WH     allocate ( wssub(mzp) )                                       ! Snow  Particles Sublimation
377! #WH     allocate ( wsmlt(mzp) )                                       ! Snow  Particles Melting
378! #WH     allocate ( wsfre(mzp) )                                       ! Rain  Drops     Freezing
379
380      END IF                                                            !
381!                                                                       !
382!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
383
384
385
386
387!     +++++++++++++++++++
388      DO      ikl=1,kcolp
389!     +++++++++++++++++++
390
391
392
393
394!  Debug
395!  ~~~~~
396! #wH   DO k=mz1_CM,mzp
397! #wH     debugH( 1:35) = 'CMiPhy: Debugged Variables: Initial'
398! #wH     debugH(36:70) = '                                   '
399! #wH     proc_1        = 'R.Hum W[%]'
400! #wH     procv1        =  0.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,               &
3442     &                      i0__CM(ipt_CM),j0__CM(ipt_CM)
3443 1037         format(/,' Ice-Crystal mPhy ',                           &
3444     &                   2x,' ',2x,1x,i2,'h',i2,'UT',                  &
3445     &                 ' -- Grid Point (',i5,',',i5,')',               &
3446     &  /,' =========================================================='&
3447     &   ,     /,'     |  z  [m] | T  [K] | qi[g/kg] |'                &
3448     &   ,            ' Ni [m-3] | Ni0[m-3] | vi [m/s] | qs[g/kg] |'   &
3449     &   ,     /,'-----+---------+--------+----------+'                &
3450     &   ,            '----------+----------+----------+----------+')
3451              write(4,1038)(k,Z___DY(ikl,k),Ta__CM(ikl,k)              &
3452     &                       ,qi__CM(ikl,k)*1.d3                       &
3453     &                       ,CCNiCM(ikl,k),Fletch(ikl,k)              &
3454     &                       ,FallVi(ikl,k),qs__CM(ikl,k)*1.d3         &
3455     &                     ,k=mz1_CM,mzp)
3456 1038         format((i4,' |' ,  f8.1,' |',f7.2,' |',f9.6,' |',        &
3457     &            2(d9.3,' |'),2(f9.6,' |')))
3458
3459        END DO
3460      END IF
3461
3462
3463
3464
3465!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3466!                                                                       !
3467! DE-ALLOCATION                                                         !
3468! =============                                                         !
3469!                                                                       !
3470      IF (FlagDALLOC)                                              THEN !
3471
3472          deallocate ( qw_io0 )                                         ! Droplets   Concentration entering CMiPhy                  [kg/kg]
3473          deallocate ( qi_io0 )                                         ! Ice  Part. Concentration entering CMiPhy                  [kg/kg]
3474          deallocate ( qs___0 )                                         ! Snow Part. Concentration entering CMiPhy                  [kg/kg]
3475! #qg     deallocate ( qg___0 )                                         ! Graupels   Concentration entering CMiPhy                  [kg/kg]
3476          deallocate ( qr___0 )                                         ! Rain Drops Concentration entering CMiPhy                  [kg/kg]
3477          deallocate ( Ta_dgC )                                         ! Air   Temperature                                           [dgC]
3478          deallocate ( sqrrro )                                         ! sqrt(roa(mzp)/roa(k))                                         [-]
3479          deallocate ( qsiEFF )                                         ! EFFective Saturation Specific Humidity over Ice           [kg/kg]
3480          deallocate ( Fletch )                                         ! Monodisperse Nb of hexagonal Plates, Fletcher (1962)          [-]
3481          deallocate ( lamdaS )                                         ! Marshall-Palmer distribution parameter for Snow Particl.
3482! #qg     deallocate ( lamdaG )                                         ! Marshall-Palmer distribution parameter for Graupels
3483          deallocate ( lamdaR )                                         ! Marshall-Palmer distribution parameter for Rain Drops
3484          deallocate ( ps_ACR )                                         ! Accretion of Snow        by Rain                   Rate [kg/kg/s]
3485          deallocate ( ps_ACW )                                         ! Accretion of Cloud Drop. by Snow Particl.          Rate [kg/kg/s]
3486          deallocate ( FallVw )                                         ! Sedimentation Velocity   of Droplets
3487          deallocate ( FallVi )                                         ! Sedimentation Velocity   of Ice  Particles
3488          deallocate ( FallVs )                                         ! Sedimentation Velocity   of Snow Particles
3489! #qg     deallocate ( FallVg )                                         ! Sedimentation Velocity   of Snow Particles
3490          deallocate ( FallVr )                                         ! Sedimentation Velocity   of Rain Drops
3491          deallocate ( qwLoss )                                         ! Mass Loss related to Sedimentation of Rain Droplets
3492          deallocate ( qiLoss )                                         ! Mass Loss related to Sedimentation of Ice  Crystals
3493          deallocate ( qsLoss )                                         ! Mass Loss related to Sedimentation of Snow Particles
3494          deallocate ( qrLoss )                                         ! Mass Loss related to Sedimentation of Rain Drops
3495! #wH     deallocate ( debugV )                                         ! Debug Variable (of 16 microphysical processes)
3496
3497! #WH     deallocate ( wihm1 )                                          ! Cloud Droplets Freezing
3498! #WH     deallocate ( wihm2 )                                          ! Ice   Crystals Homogeneous Sublimation
3499! #WH     deallocate ( wicnd )                                          ! Ice   Crystals Nucleation              (Emde & Kahlig)
3500! #WH     deallocate ( widep )                                          ! Ice   Crystals Growth Bergeron Process (Emde & Kahlig)
3501! #WH     deallocate ( wisub )                                          ! Ice   Crystals             Sublimation (Levkov)
3502! #WH     deallocate ( wimlt )                                          ! Ice   Crystals Melting 
3503! #WH     deallocate ( wwevp )                                          ! Water Vapor Condensation / Evaporation (Fractional Cloudiness)
3504! #WH     deallocate ( wraut )                                          ! Cloud Droplets AUTO-Conversion
3505! #WH     deallocate ( wsaut )                                          ! Ice   Crystals AUTO-Conversion
3506! #WH     deallocate ( wracw )                                          ! Accretion of Cloud Droplets by Rain, Ta > 0, --> Rain
3507! #WH     deallocate ( wsacw )                                          ! Accretion of Cloud Droplets by Rain, Ta < 0, --> Snow
3508! #WH     deallocate ( wsaci )                                          ! Accretion of Ice   Crystals by Snow          --> Snow
3509! #WH     deallocate ( wraci )                                          ! Accretion of Ice   Crystals by Rain          --> Snow
3510! #WH     deallocate ( wiacr )                                          ! Accretion of Rain by Ice   Crystals          --> Snow
3511! #WH     deallocate ( wsacr )                                          ! Accretion of Rain by Snow                    --> Snow
3512! #WH     deallocate ( wracs )                                          ! Accretion of Snow by Rain                    --> Snow, Rain
3513! #WH     deallocate ( wrevp )                                          ! Rain  Drops     Evaporation 
3514! #WH     deallocate ( wssub )                                          ! Snow  Particles Sublimation
3515! #WH     deallocate ( wsmlt )                                          ! Snow  Particles Melting
3516! #WH     deallocate ( wsfre )                                          ! Rain  Drops     Freezing
3517!                                                                       !
3518      END IF                                                            !
3519!                                                                       !
3520!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3521
3522
3523
3524
3525      return
3526      end subroutine CMiPhy
Note: See TracBrowser for help on using the repository browser.