Philippe 01/12/2015: added missing & at the beginning of output message lines
[MNH-git_open_source-lfs.git] / src / MNH / ini_rain_ice.f90
1 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
4 !MNH_LIC for details. version 1.
5 !-----------------------------------------------------------------
6 !--------------- special set of characters for RCS information
7 !-----------------------------------------------------------------
8 ! $Source$ $Revision$
9 ! masdev4_7 BUG1 2007/06/15 17:47:18
10 !-----------------------------------------------------------------
11 !      ########################
12        MODULE MODI_INI_RAIN_ICE
13 !      ########################
14 !
15 INTERFACE
16       SUBROUTINE INI_RAIN_ICE ( KLUOUT, PTSTEP, PDZMIN, KSPLITR, HCLOUD )
17 !
18 INTEGER,                 INTENT(IN) :: KLUOUT    ! Logical unit number for prints
19 INTEGER,                 INTENT(OUT):: KSPLITR   ! Number of small time step
20                                                  ! integration for  rain
21                                                  ! sedimendation
22 !
23 REAL,                    INTENT(IN) :: PTSTEP    ! Effective Time step
24 !
25 REAL,                    INTENT(IN) :: PDZMIN    ! minimun vertical mesh size
26 !
27 CHARACTER (LEN=4), INTENT(IN)       :: HCLOUD    ! Indicator of the cloud scheme
28 !
29 END SUBROUTINE INI_RAIN_ICE
30 !
31 END INTERFACE
32 !
33 END MODULE MODI_INI_RAIN_ICE
34 !     ###########################################################
35       SUBROUTINE INI_RAIN_ICE ( KLUOUT, PTSTEP, PDZMIN, KSPLITR, HCLOUD )
36 !     ###########################################################
37 !
38 !!****  *INI_RAIN_ICE * - initialize the constants necessary for the warm and
39 !!                        cold microphysical schemes.
40 !!
41 !!    PURPOSE
42 !!    -------
43 !!      The purpose of this routine is to initialize the constants used to 
44 !!    resolve the mixed phase microphysical scheme. The collection kernels of 
45 !!    the precipitating particles are recomputed if necessary if some parameters
46 !!    defining the ice categories have been modified. The number of small
47 !!    time steps leading to stable scheme for the rain, ice, snow and ggraupeln
48 !!    sedimentation is also computed (time-splitting technique).     
49 !!
50 !!**  METHOD
51 !!    ------
52 !!      The constants are initialized to their numerical values and the number
53 !!    of small time step is computed by dividing the 2* Deltat time interval of
54 !!    the Leap-frog scheme so that the stability criterion for the rain
55 !!    sedimentation is fulfilled for a Raindrop maximal fall velocity equal
56 !!    VTRMAX. The parameters defining the collection kernels are read and are
57 !!    checked against the new ones. If any change occurs, these kernels are
58 !!    recomputed and their numerical values are written in the output listiing. 
59 !!
60 !!    EXTERNAL
61 !!    --------
62 !!      GAMMA    :  gamma function
63 !!     
64 !!
65 !!    IMPLICIT ARGUMENTS
66 !!    ------------------
67 !!      Module MODD_CST
68 !!        XPI                  !
69 !!        XP00                 ! Reference pressure
70 !!        XRD                  ! Gaz constant for dry air
71 !!        XRHOLW               ! Liquid water density
72 !!      Module MODD_REF
73 !!        XTHVREFZ             ! Reference virtual pot.temp. without orography
74 !!      Module MODD_PARAMETERS
75 !!        JPVEXT               !
76 !!      Module MODD_RAIN_ICE_DESCR
77 !!      Module MODD_RAIN_ICE_PARAM
78 !!
79 !!    REFERENCE
80 !!    ---------
81 !!      Book2 of documentation ( routine INI_RAIN_ICE )
82 !!
83 !!    AUTHOR
84 !!    ------
85 !!      J.-P. Pinty      * Laboratoire d'Aerologie*
86 !!
87 !!    MODIFICATIONS
88 !!    -------------
89 !!      Original    04/12/95
90 !!      J.-P. Pinty 05/04/96 Add automatic control and regeneration of the
91 !!                           collection kernels
92 !!      J.-P. Pinty 10/05/96 Correction of ZRATE and computations of RIM
93 !!      J.-P. Pinty 24/11/97 Sedimentation of ice made for Columns and bug for XAG
94 !!      J.-P. Lafore 23/11/98 Back to Lin et al. 83 formulation for RIAUTS
95 !!                            with a Critical ice content set to .5 g/Kg
96 !!      N. Asencio  13/08/98 parallel code: PDZMIN is computed outside in ini_modeln
97 !!      J.-P. Lafore 12/8/98 In case of nesting microphysics constants of
98 !!                           MODD_RAIN_ICE_PARAM are computed only once.
99 !!                           Only KSPLTR is computed for each model.
100 !!      J. Stein    20/04/99 remove 2 unused local variables
101 !!      G Molinie   21/05/99 Bug in XEXRCFRI and XRCFRI
102 !!      J.-P. Pinty 24/06/00 Bug in RCRIMS
103 !!      J.-P. Pinty 24/12/00 Update hail case
104 !!      J.-P. Chaboureau & J.-P. Pinty
105 !!                  24/03/01 Update XCRIAUTI for cirrus cases
106 !!      J.-P. Pinty 24/11/01 Update ICE3/ICE4 options
107 !!
108 !-------------------------------------------------------------------------------
109 !
110 !*       0.    DECLARATIONS
111 !              ------------
112 !
113 USE MODE_FM
114 USE MODD_CST
115 USE MODD_LUNIT
116 USE MODD_PARAMETERS
117 USE MODD_PARAM_ICE
118 USE MODD_RAIN_ICE_DESCR
119 USE MODD_RAIN_ICE_PARAM
120 USE MODD_REF
121 !
122 USE MODI_GAMMA
123 USE MODI_GAMMA_INC
124 USE MODI_RRCOLSS
125 USE MODI_RZCOLX
126 USE MODI_RSCOLRG
127 USE MODI_READ_XKER_RACCS
128 USE MODI_READ_XKER_SDRYG
129 USE MODI_READ_XKER_RDRYG
130 USE MODI_READ_XKER_SWETH
131 USE MODI_READ_XKER_GWETH
132 !
133 IMPLICIT NONE
134 !
135 !*       0.1   Declarations of dummy arguments :
136 !
137 !
138 INTEGER,                 INTENT(IN) :: KLUOUT   ! Logical unit number for prints
139 INTEGER,                 INTENT(OUT):: KSPLITR   ! Number of small time step
140                                                  ! integration for  rain
141                                                  ! sedimendation
142 !
143 REAL,                    INTENT(IN) :: PTSTEP    ! Effective Time step 
144 !
145 REAL,                    INTENT(IN) :: PDZMIN    ! minimun vertical mesh size
146 !
147 CHARACTER (LEN=4), INTENT(IN)       :: HCLOUD    ! Indicator of the cloud scheme
148 !
149 !
150 !
151 !*       0.2   Declarations of local variables :
152 !
153 INTEGER :: IKB                ! Coordinates of the first physical 
154                               ! points along z
155 INTEGER :: J1,J2              ! Internal loop indexes
156 REAL :: ZT                    ! Work variable
157 REAL :: ZVTRMAX               ! Raindrop maximal fall velocity
158 REAL :: ZRHO00                ! Surface reference air density
159 REAL :: ZRATE                 ! Geometrical growth of Lbda in the tabulated
160                               ! functions and kernels
161 REAL :: ZBOUND                ! XDCSLIM*Lbda_s: upper bound for the partial
162                               ! integration of the riming rate of the aggregates
163 REAL :: ZEGS, ZEGR, ZEHS, ZEHG! Bulk collection efficiencies
164 !
165 INTEGER :: IND                ! Number of interval to integrate the kernels
166 REAL :: ZESR                  ! Mean efficiency of rain-aggregate collection
167 REAL :: ZFDINFTY              ! Factor used to define the "infinite" diameter
168 !
169 !
170 !
171 LOGICAL  :: GFLAG   ! Logical flag for printing the constatnts on the output
172                     ! listing
173 REAL     :: ZCONC_MAX ! Maximal concentration for snow
174 REAL     :: ZGAMC,ZGAMC2 ! parameters
175                     ! involving various moments of the generalized gamma law
176 REAL     :: ZFACT_NUCL! Amplification factor for the minimal ice concentration
177 REAL     :: ZXR     ! Value of x_r in N_r = C_r lambda_r ** x_r
178 !  
179 INTEGER  :: KND
180 INTEGER  :: KACCLBDAS,KACCLBDAR,KDRYLBDAG,KDRYLBDAS,KDRYLBDAR
181 INTEGER  :: KWETLBDAS,KWETLBDAG,KWETLBDAH
182 REAL     :: PALPHAR,PALPHAS,PALPHAG,PALPHAH
183 REAL     :: PNUR,PNUS,PNUG,PNUH
184 REAL     :: PBR,PBS,PBG,PBH
185 REAL     :: PCR,PCS,PCG,PCH
186 REAL     :: PDR,PDS,PDG,PDH
187 REAL     :: PESR,PEGS,PEGR,PEHS,PEHG
188 REAL     :: PFDINFTY
189 REAL     :: PACCLBDAS_MAX,PACCLBDAR_MAX,PACCLBDAS_MIN,PACCLBDAR_MIN
190 REAL     :: PDRYLBDAG_MAX,PDRYLBDAS_MAX,PDRYLBDAG_MIN,PDRYLBDAS_MIN
191 REAL     :: PDRYLBDAR_MAX,PDRYLBDAR_MIN
192 REAL     :: PWETLBDAS_MAX,PWETLBDAG_MAX,PWETLBDAS_MIN,PWETLBDAG_MIN
193 REAL     :: PWETLBDAH_MAX,PWETLBDAH_MIN
194 !-------------------------------------------------------------------------------
195 !
196 !
197 !*       0.     FUNCTION STATEMENTS
198 !               -------------------
199 !
200 !
201 !*       0.1    p_moment of the Generalized GAMMA function
202 !
203 !
204 !
205 !        1.     COMPUTE KSPLTR FOR EACH MODEL
206 !               ---------------------------------------------------------
207 !
208 !*       1.1    Set the hailstones maximum fall velocity
209 !
210 IF (CSEDIM == 'SPLI') THEN
211  IF (HCLOUD == 'ICE4') THEN
212   ZVTRMAX = 40.
213  ELSE IF (HCLOUD == 'ICE3') THEN
214   ZVTRMAX = 10.
215  END IF 
216 END IF
217 !
218 !*       1.2    Compute the number of small time step integration
219 !
220 KSPLITR = 1
221 IF (CSEDIM == 'SPLI') THEN
222  SPLIT : DO
223   ZT = PTSTEP / FLOAT(KSPLITR)
224   IF ( ZT * ZVTRMAX / PDZMIN .LT. 1.) EXIT SPLIT
225   KSPLITR = KSPLITR + 1
226  END DO SPLIT
227 END IF
228 !
229 IF (ALLOCATED(XRTMIN)) THEN       ! In case of nesting microphysics constants of
230                                   ! MODD_RAIN_ICE_PARAM are computed only once,
231                                   ! but if INI_RAIN_ICE has been called already
232                                   ! one must change the XRTMIN size.
233   DEALLOCATE(XRTMIN)
234 END IF
235 !
236 !-------------------------------------------------------------------------------
237 !
238 !*       2.     CHARACTERISTICS OF THE SPECIES
239 !               ------------------------------
240 !
241 !
242 !*       2.1    Cloud droplet and Raindrop characteristics
243 !
244 XAC = (XPI/6.0)*XRHOLW
245 XBC = 3.0
246 XCC = XRHOLW*XG/(18.0*1.7E-5) ! Stokes flow (Pruppacher p 322 for T=273K)
247 XDC = 2.0
248 !
249 !
250 XAR = (XPI/6.0)*XRHOLW
251 XBR = 3.0
252 XCR = 842.
253 XDR = 0.8
254 !
255 !XCCR = 1.E7    ! N0_r =  XCXR * lambda_r ** ZXR
256 XCCR = 8.E6    ! N0_r =  XCXR * lambda_r ** ZXR
257 ZXR  = -1.     !
258 !
259 XF0R = 1.00
260 XF1R = 0.26
261 !
262 XC1R = 1./2.
263 !
264 !
265 !*       2.2    Ice crystal characteristics
266 !
267 !
268 SELECT CASE (CPRISTINE_ICE)
269   CASE('PLAT')
270     XAI = 0.82      ! Plates
271     XBI = 2.5       ! Plates 
272     XC_I = 800.     ! Plates
273     XDI = 1.0       ! Plates
274     XC1I = 1./XPI   ! Plates
275   CASE('COLU')
276     XAI = 2.14E-3   ! Columns
277     XBI = 1.7       ! Columns
278     XC_I = 2.1E5    ! Columns
279     XDI = 1.585     ! Columns
280     XC1I = 0.8      ! Columns 
281   CASE('BURO')
282     XAI = 44.0      ! Bullet rosettes
283     XBI = 3.0       ! Bullet rosettes
284     XC_I = 4.3E5    ! Bullet rosettes
285     XDI = 1.663     ! Bullet rosettes
286     XC1I = 0.5      ! Bullet rosettes
287 END SELECT
288 !
289 !  Note that XCCI=N_i (a locally predicted value) and XCXI=0.0, implicitly
290 !
291 XF0I = 1.00
292 XF2I = 0.14
293 !
294 !
295 !*       2.3    Snowflakes/aggregates characteristics
296 !
297 !
298 XAS = 0.02
299 XBS = 1.9
300 XCS = 5.1
301 XDS = 0.27
302 !
303 XCCS = 5.0
304 XCXS = 1.0
305 !
306 XF0S = 0.86
307 XF1S = 0.28
308 !
309 XC1S = 1./XPI
310 !
311 !
312 !*       2.4    Graupel/Frozen drop characteristics
313 !
314 !
315 XAG = 19.6  ! Lump graupel case
316 XBG = 2.8   ! Lump graupel case
317 XCG = 124.  ! Lump graupel case
318 XDG = 0.66  ! Lump graupel case
319 !
320 XCCG = 5.E5
321 XCXG = -0.5
322 ! XCCG = 4.E4 ! Test of Ziegler (1988)
323 ! XCXG = -1.0 ! Test of Ziegler (1988)
324 !
325 XF0G = 0.86
326 XF1G = 0.28
327 !
328 XC1G = 1./2.
329 !
330 !
331 !*       2.5    Hailstone characteristics
332 !
333 !
334 XAH = 470.
335 XBH = 3.0
336 XCH = 207.
337 XDH = 0.64
338 !
339 !XCCH = 5.E-4
340 !XCXH = 2.0
341 !!!!!!!!!!!!
342    XCCH = 4.E4 ! Test of Ziegler (1988)
343    XCXH = -1.0 ! Test of Ziegler (1988)
344 !!!    XCCH = 5.E5 ! Graupel_like
345 !!!    XCXH = -0.5 ! Graupel_like
346 !!!!!!!!!!!!
347 !
348 XF0H = 0.86
349 XF1H = 0.28
350 !
351 XC1H = 1./2.
352 !
353 !-------------------------------------------------------------------------------
354 !
355 !*       3.     DIMENSIONAL DISTRIBUTIONS OF THE SPECIES
356 !                   ----------------------------------------
357 !
358 !
359 !        3.1    Cloud droplet distribution
360 !
361 ! Over land
362 XALPHAC = 1.0  ! Gamma law of the Cloud droplet (here volume-like distribution)
363 XNUC    = 3.0  ! Gamma law with little dispersion
364 !
365 !
366 ! Over sea
367 XALPHAC2 = 3.0  ! Gamma law of the Cloud droplet (here volume-like distribution)
368 XNUC2    = 1.0  ! Gamma law with little dispersion
369 !
370 !*       3.2    Raindrops distribution
371 !
372 XALPHAR = 1.0  ! Exponential law
373 XNUR    = 1.0  ! Exponential law
374 !
375 !*       3.3    Ice crystal distribution
376 !
377 XALPHAI = 3.0  ! Gamma law for the ice crystal volume
378 XNUI    = 3.0  ! Gamma law with little dispersion
379 !
380 XALPHAS = 1.0  ! Exponential law
381 XNUS    = 1.0  ! Exponential law
382 !
383 XALPHAG = 1.0  ! Exponential law
384 XNUG    = 1.0  ! Exponential law
385 !
386 XALPHAH = 1.0  ! Gamma law
387 XNUH    = 8.0  ! Gamma law with little dispersion
388 !
389 !*       3.4    Constants for shape parameter
390 !
391 ZGAMC = MOMG(XALPHAC,XNUC,3.)
392 ZGAMC2 = MOMG(XALPHAC2,XNUC2,3.)
393 XLBC(1)   = XAR*ZGAMC
394 XLBC(2)   = XAR*ZGAMC2
395 XLBEXC = 1.0/XBC
396 !
397 XLBEXR = 1.0/(-1.0-XBR)
398 XLBR   = ( XAR*XCCR*MOMG(XALPHAR,XNUR,XBR) )**(-XLBEXR)
399 !
400 XLBEXI = 1.0/(-XBI)
401 XLBI   = ( XAI*MOMG(XALPHAI,XNUI,XBI) )**(-XLBEXI)
402 !
403 XLBEXS = 1.0/(XCXS-XBS)
404 XLBS   = ( XAS*XCCS*MOMG(XALPHAS,XNUS,XBS) )**(-XLBEXS)
405 !
406 XLBEXG = 1.0/(XCXG-XBG)
407 XLBG   = ( XAG*XCCG*MOMG(XALPHAG,XNUG,XBG) )**(-XLBEXG)
408 !
409 XLBEXH = 1.0/(XCXH-XBH)
410 XLBH   = ( XAH*XCCH*MOMG(XALPHAH,XNUH,XBH) )**(-XLBEXH)
411 !
412 !*       3.5    Minimal values allowed for the mixing ratios
413 !
414 XLBDAR_MAX = 100000.0
415 XLBDAS_MAX = 100000.0
416 XLBDAG_MAX = 100000.0
417 !
418 ZCONC_MAX  = 1.E6 ! Maximal concentration for falling particules set to 1 per cc
419 XLBDAS_MAX = ( ZCONC_MAX/XCCS )**(1./XCXS)
420 !
421 IF (HCLOUD == 'ICE4') THEN
422   ALLOCATE( XRTMIN(7) )
423 ELSE IF (HCLOUD == 'ICE3') THEN
424   ALLOCATE( XRTMIN(6) )
425 END IF
426 !
427 XRTMIN(1) = 1.0E-20
428 XRTMIN(2) = 1.0E-20
429 XRTMIN(3) = 1.0E-20
430 XRTMIN(4) = 1.0E-20
431 XRTMIN(5) = 1.0E-15
432 XRTMIN(6) = 1.0E-15
433 IF (HCLOUD == 'ICE4') XRTMIN(7) = 1.0E-15
434 !
435 XCONC_SEA=1E8 ! 100/cm3
436 XCONC_LAND=3E8 ! 300/cm3
437 XCONC_URBAN=5E8 ! 500/cm3
438 !
439 !-------------------------------------------------------------------------------
440 !
441 !*       4.     CONSTANTS FOR THE SEDIMENTATION
442 !               -------------------------------
443 !
444 !
445 !*       4.1    Exponent of the fall-speed air density correction
446 !
447 XCEXVT = 0.4
448 !
449 IKB = 1 + JPVEXT
450 ZRHO00 = XP00/(XRD*XTHVREFZ(IKB))
451 !
452 !*       4.2    Constants for sedimentation
453 !
454 XFSEDC(1)  = GAMMA(XNUC+(XDC+3.)/XALPHAC)/GAMMA(XNUC+3./XALPHAC)*     &
455             (ZRHO00)**XCEXVT
456 XFSEDC(2)  = GAMMA(XNUC2+(XDC+3.)/XALPHAC2)/GAMMA(XNUC2+3./XALPHAC2)*     &
457             (ZRHO00)**XCEXVT
458 !
459 XEXSEDR = (XBR+XDR+1.0)/(XBR+1.0)
460 XFSEDR  = XCR*XAR*XCCR*MOMG(XALPHAR,XNUR,XBR+XDR)*                         &
461             (XAR*XCCR*MOMG(XALPHAR,XNUR,XBR))**(-XEXSEDR)*(ZRHO00)**XCEXVT
462 !
463 XEXRSEDI = (XBI+XDI)/XBI
464 XEXCSEDI = 1.0-XEXRSEDI
465 XFSEDI   = (4.*XPI*900.)**(-XEXCSEDI) *                         &
466            XC_I*XAI*MOMG(XALPHAI,XNUI,XBI+XDI) *                &
467            ((XAI*MOMG(XALPHAI,XNUI,XBI)))**(-XEXRSEDI) *        &
468            (ZRHO00)**XCEXVT
469 !
470 !  Computations made for Columns
471 !
472 XEXRSEDI = 1.9324
473 XEXCSEDI =-0.9324
474 XFSEDI   = 3.89745E11*MOMG(XALPHAI,XNUI,3.285)*                          &
475                       MOMG(XALPHAI,XNUI,1.7)**(-XEXRSEDI)*(ZRHO00)**XCEXVT
476 XEXCSEDI =-0.9324*3.0
477 WRITE (KLUOUT,FMT=*)' PRISTINE ICE SEDIMENTATION for columns XFSEDI =',XFSEDI
478 !
479 !
480 XEXSEDS = (XBS+XDS-XCXS)/(XBS-XCXS)
481 XFSEDS  = XCS*XAS*XCCS*MOMG(XALPHAS,XNUS,XBS+XDS)*                         &
482             (XAS*XCCS*MOMG(XALPHAS,XNUS,XBS))**(-XEXSEDS)*(ZRHO00)**XCEXVT
483 !
484 XEXSEDG = (XBG+XDG-XCXG)/(XBG-XCXG)
485 XFSEDG  = XCG*XAG*XCCG*MOMG(XALPHAG,XNUG,XBG+XDG)*                         &
486             (XAG*XCCG*MOMG(XALPHAG,XNUG,XBG))**(-XEXSEDG)*(ZRHO00)**XCEXVT
487 !
488 XEXSEDH = (XBH+XDH-XCXH)/(XBH-XCXH)
489 XFSEDH  = XCH*XAH*XCCH*MOMG(XALPHAH,XNUH,XBH+XDH)*                         &
490             (XAH*XCCH*MOMG(XALPHAH,XNUH,XBH))**(-XEXSEDH)*(ZRHO00)**XCEXVT
491 !
492 !
493 !-------------------------------------------------------------------------------
494 !
495 !*       5.     CONSTANTS FOR THE SLOW COLD PROCESSES
496 !               -------------------------------------
497 !
498 !
499 !*       5.1    Constants for ice nucleation
500 !
501 SELECT CASE (CPRISTINE_ICE)
502   CASE('PLAT')
503     ZFACT_NUCL =  1.0    ! Plates
504   CASE('COLU')
505     ZFACT_NUCL = 25.0    ! Columns
506   CASE('BURO')
507     ZFACT_NUCL = 17.0    ! Bullet rosettes
508 END SELECT
509 !
510 XNU10 = 50.*ZFACT_NUCL
511 XALPHA1 = 4.5
512 XBETA1 = 0.6
513 !
514 XNU20 = 1000.*ZFACT_NUCL
515 XALPHA2 = 12.96
516 XBETA2 = 0.639
517 !
518 XMNU0 = 6.88E-13
519 !
520 GFLAG = .TRUE.
521 IF (GFLAG) THEN
522   WRITE(UNIT=KLUOUT,FMT='("      Heterogeneous nucleation")')
523   WRITE(UNIT=KLUOUT,FMT='(" NU10=",E13.6," ALPHA1=",E13.6," BETA1=",E13.6)') &
524                                                       XNU10,XALPHA1,XBETA1
525   WRITE(UNIT=KLUOUT,FMT='(" NU20=",E13.6," ALPHA2=",E13.6," BETA2=",E13.6)') &
526                                                       XNU20,XALPHA2,XBETA2
527   WRITE(UNIT=KLUOUT,FMT='(" mass of embryo XMNU0=",E13.6)') XMNU0
528 END IF
529 !
530 XALPHA3 = -3.075
531 XBETA3 = 81.00356
532 XHON = (XPI/6.)*((2.0*3.0*4.0*5.0*6.0)/(2.0*3.0))*(1.1E5)**(-3.0) !
533                                        ! Pi/6 * (G_c(6)/G_c(3)) * (1/Lbda_c**3)
534                                        ! avec Lbda_c=1.1E5 m^-1
535                                        !     the formula is equivalent to
536                                        !        rho_dref * r_c     G(6)
537                                        ! Pi/6 * -------------- * ---------
538                                        !         rho_lw * N_c    G(3)*G(3)
539 !
540 GFLAG = .TRUE.
541 IF (GFLAG) THEN
542   WRITE(UNIT=KLUOUT,FMT='("      Homogeneous nucleation")')
543   WRITE(UNIT=KLUOUT,FMT='(" ALPHA3=",E13.6," BETA3=",E13.6)') XALPHA3,XBETA3
544   WRITE(UNIT=KLUOUT,FMT='(" constant XHON=",E13.6)') XHON
545 END IF
546 !
547 !
548 !*       5.2    Constants for vapor deposition on ice
549 !
550 XSCFAC = (0.63**(1./3.))*SQRT((ZRHO00)**XCEXVT) ! One assumes Sc=0.63
551 !
552 X0DEPI = (4.0*XPI)*XC1I*XF0I*MOMG(XALPHAI,XNUI,1.)
553 X2DEPI = (4.0*XPI)*XC1I*XF2I*XC_I*MOMG(XALPHAI,XNUI,XDI+2.0)
554 !
555 X0DEPS = (4.0*XPI)*XCCS*XC1S*XF0S*MOMG(XALPHAS,XNUS,1.)
556 X1DEPS = (4.0*XPI)*XCCS*XC1S*XF1S*SQRT(XCS)*MOMG(XALPHAS,XNUS,0.5*XDS+1.5)
557 XEX0DEPS = XCXS-1.0
558 XEX1DEPS = XCXS-0.5*(XDS+3.0)
559 !
560 X0DEPG = (4.0*XPI)*XCCG*XC1G*XF0G*MOMG(XALPHAG,XNUG,1.)
561 X1DEPG = (4.0*XPI)*XCCG*XC1G*XF1G*SQRT(XCG)*MOMG(XALPHAG,XNUG,0.5*XDG+1.5)
562 XEX0DEPG = XCXG-1.0
563 XEX1DEPG = XCXG-0.5*(XDG+3.0)
564 !
565 X0DEPH = (4.0*XPI)*XCCH*XC1H*XF0H*MOMG(XALPHAH,XNUH,1.)
566 X1DEPH = (4.0*XPI)*XCCH*XC1H*XF1H*SQRT(XCH)*MOMG(XALPHAH,XNUH,0.5*XDH+1.5)
567 XEX0DEPH = XCXH-1.0
568 XEX1DEPH = XCXH-0.5*(XDH+3.0)
569 !
570 !*       5.3    Constants for pristine ice autoconversion
571 !
572 XTIMAUTI = 1.E-3  !  Time constant at T=T_t
573 XTEXAUTI = 0.015  !  Temperature factor of the I+I collection efficiency
574 XCRIAUTI = 0.2E-4 !  Critical ice content for the autoconversion to occur
575                   !  Revised value by Chaboureau and Pinty (2006)
576 XACRIAUTI = 0.06
577 XBCRIAUTI = -3.5
578 XT0CRIAUTI= (LOG10(XCRIAUTI)-XBCRIAUTI)/0.06
579 !
580 GFLAG = .TRUE.
581 IF (GFLAG) THEN
582   WRITE(UNIT=KLUOUT,FMT='("      pristine ice autoconversion")')
583   WRITE(UNIT=KLUOUT,FMT='(" Time constant   XTIMAUTI=",E13.6)') XTIMAUTI
584   WRITE(UNIT=KLUOUT,FMT='(" Temp. factor    XTEXAUTI=",E13.6)') XTEXAUTI
585   WRITE(UNIT=KLUOUT,FMT='(" Crit. ice cont. XCRIAUTI=",E13.6)') XCRIAUTI
586   WRITE(UNIT=KLUOUT,FMT='(" A Coef. for cirrus law XACRIAUTI=",E13.6)') XACRIAUTI
587   WRITE(UNIT=KLUOUT,FMT='(" B Coef. for cirrus law XBCRIAUTI=",E13.6)') XBCRIAUTI
588   WRITE(UNIT=KLUOUT,FMT='(" Temp degC at which cirrus law starts to be &
589                             & used=",E13.6)') XT0CRIAUTI
590 END IF
591 !
592 !
593 !*       5.4    Constants for snow aggregation
594 !
595 XCOLIS   = 0.25 ! Collection efficiency of I+S
596 XCOLEXIS = 0.05 ! Temperature factor of the I+S collection efficiency
597 XFIAGGS  = (XPI/4.0)*XCOLIS*XCCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
598 XEXIAGGS = XCXS-XDS-2.0
599 !
600 GFLAG = .TRUE.
601 IF (GFLAG) THEN
602   WRITE(UNIT=KLUOUT,FMT='("      snow aggregation")')
603   WRITE(UNIT=KLUOUT,FMT='(" Coll. efficiency XCOLIS=",E13.6)') XCOLIS
604   WRITE(UNIT=KLUOUT,FMT='(" Temp. factor     XCOLEXIS=",E13.6)') XCOLEXIS
605 END IF
606 !
607 !
608 !-------------------------------------------------------------------------------
609 !
610 !*       6.     CONSTANTS FOR THE SLOW WARM PROCESSES
611 !               -------------------------------------
612 !
613 !
614 !*       6.1    Constants for the cloud droplets autoconversion
615 !
616 XTIMAUTC = 1.E-3
617 XCRIAUTC = 0.5E-3 
618 !
619 GFLAG = .TRUE.
620 IF (GFLAG) THEN
621   WRITE(UNIT=KLUOUT,FMT='("      cloud droplets autoconversion")')
622   WRITE(UNIT=KLUOUT,FMT='(" Time constant   XTIMAUTC=",E13.6)') XTIMAUTC
623   WRITE(UNIT=KLUOUT,FMT='(" Crit. ice cont. XCRIAUTC=",E13.6)') XCRIAUTC
624 END IF
625 !
626 !*       6.2    Constants for the accretion of cloud droplets by raindrops
627 !
628 XFCACCR  = (XPI/4.0)*XCCR*XCR*(ZRHO00**XCEXVT)*MOMG(XALPHAR,XNUR,XDR+2.0)
629 XEXCACCR = -XDR-3.0
630 !
631 !*       6.3    Constants for the evaporation of the raindrops
632 !
633 X0EVAR = (4.0*XPI)*XCCR*XC1R*XF0R*MOMG(XALPHAR,XNUR,1.)
634 X1EVAR = (4.0*XPI)*XCCR*XC1R*XF1R*SQRT(XCR)*MOMG(XALPHAR,XNUR,0.5*XDR+1.5)
635 XEX0EVAR = -2.0
636 XEX1EVAR = -1.0-0.5*(XDR+3.0)
637 !
638 !
639 !-------------------------------------------------------------------------------
640 !
641 !*       7.     CONSTANTS FOR THE FAST COLD PROCESSES FOR THE AGGREGATES
642 !               --------------------------------------------------------
643 !
644 !
645 !*       7.1    Constants for the riming of the aggregates
646 !
647 XDCSLIM  = 0.007 ! D_cs^lim = 7 mm as suggested by Farley et al. (1989)
648 XCOLCS   = 1.0
649 XEXCRIMSS= XCXS-XDS-2.0
650 XCRIMSS  = (XPI/4.0)*XCOLCS*XCCS*XCS*(ZRHO00**XCEXVT)*MOMG(XALPHAS,XNUS,XDS+2.0)
651 XEXCRIMSG= XEXCRIMSS
652 XCRIMSG  = XCRIMSS
653 XSRIMCG  = XCCS*XAS*MOMG(XALPHAS,XNUS,XBS)
654 XEXSRIMCG= XCXS-XBS
655 !
656 GFLAG = .TRUE.
657 IF (GFLAG) THEN
658   WRITE(UNIT=KLUOUT,FMT='("      riming of the aggregates")')
659   WRITE(UNIT=KLUOUT,FMT='(" D_cs^lim (Farley et al.) XDCSLIM=",E13.6)') XDCSLIM
660   WRITE(UNIT=KLUOUT,FMT='(" Coll. efficiency          XCOLCS=",E13.6)') XCOLCS
661 END IF
662 !
663 NGAMINC = 80
664 XGAMINC_BOUND_MIN = 1.0E-1 ! Minimal value of (Lbda * D_cs^lim)**alpha
665 XGAMINC_BOUND_MAX = 1.0E7  ! Maximal value of (Lbda * D_cs^lim)**alpha
666 ZRATE = EXP(LOG(XGAMINC_BOUND_MAX/XGAMINC_BOUND_MIN)/FLOAT(NGAMINC-1))
667 !
668 IF( .NOT.ALLOCATED(XGAMINC_RIM1) ) ALLOCATE( XGAMINC_RIM1(NGAMINC) )
669 IF( .NOT.ALLOCATED(XGAMINC_RIM2) ) ALLOCATE( XGAMINC_RIM2(NGAMINC) )
670 !
671 DO J1=1,NGAMINC
672   ZBOUND = XGAMINC_BOUND_MIN*ZRATE**(J1-1)
673   XGAMINC_RIM1(J1) = GAMMA_INC(XNUS+(2.0+XDS)/XALPHAS,ZBOUND)
674   XGAMINC_RIM2(J1) = GAMMA_INC(XNUS+XBS/XALPHAS      ,ZBOUND)
675 END DO
676 !
677 XRIMINTP1 = XALPHAS / LOG(ZRATE)
678 XRIMINTP2 = 1.0 + XRIMINTP1*LOG( XDCSLIM/(XGAMINC_BOUND_MIN)**(1.0/XALPHAS) )
679 !
680 !*       7.2    Constants for the accretion of raindrops onto aggregates
681 !
682 XFRACCSS = ((XPI**2)/24.0)*XCCS*XCCR*XRHOLW*(ZRHO00**XCEXVT)
683 !
684 XLBRACCS1   =    MOMG(XALPHAS,XNUS,2.)*MOMG(XALPHAR,XNUR,3.)
685 XLBRACCS2   = 2.*MOMG(XALPHAS,XNUS,1.)*MOMG(XALPHAR,XNUR,4.)
686 XLBRACCS3   =                          MOMG(XALPHAR,XNUR,5.)
687 !
688 XFSACCRG = (XPI/4.0)*XAS*XCCS*XCCR*(ZRHO00**XCEXVT)
689 !
690 XLBSACCR1   =    MOMG(XALPHAR,XNUR,2.)*MOMG(XALPHAS,XNUS,XBS)
691 XLBSACCR2   = 2.*MOMG(XALPHAR,XNUR,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
692 XLBSACCR3   =                          MOMG(XALPHAS,XNUS,XBS+2.)
693 !
694 !*       7.2.1  Defining the ranges for the computation of the kernels
695 !
696 ! Notice: One magnitude of lambda discretized over 10 points for rain
697 ! Notice: One magnitude of lambda discretized over 10 points for snow
698 !
699 NACCLBDAS = 40
700 XACCLBDAS_MIN = 5.0E1 ! Minimal value of Lbda_s to tabulate XKER_RACCS
701 XACCLBDAS_MAX = 5.0E5 ! Maximal value of Lbda_s to tabulate XKER_RACCS
702 ZRATE = LOG(XACCLBDAS_MAX/XACCLBDAS_MIN)/FLOAT(NACCLBDAS-1)
703 XACCINTP1S = 1.0 / ZRATE
704 XACCINTP2S = 1.0 - LOG( XACCLBDAS_MIN ) / ZRATE
705 NACCLBDAR = 40
706 XACCLBDAR_MIN = 1.0E3 ! Minimal value of Lbda_r to tabulate XKER_RACCS
707 XACCLBDAR_MAX = 1.0E7 ! Maximal value of Lbda_r to tabulate XKER_RACCS
708 ZRATE = LOG(XACCLBDAR_MAX/XACCLBDAR_MIN)/FLOAT(NACCLBDAR-1)
709 XACCINTP1R = 1.0 / ZRATE
710 XACCINTP2R = 1.0 - LOG( XACCLBDAR_MIN ) / ZRATE
711 !
712 !*       7.2.2  Computations of the tabulated normalized kernels
713 !
714 IND      = 50    ! Interval number, collection efficiency and infinite diameter
715 ZESR     = 1.0   ! factor used to integrate the dimensional distributions when
716 ZFDINFTY = 20.0  ! computing the kernels XKER_RACCSS, XKER_RACCS and XKER_SACCRG
717 !
718 IF( .NOT.ALLOCATED(XKER_RACCSS) ) ALLOCATE( XKER_RACCSS(NACCLBDAS,NACCLBDAR) )
719 IF( .NOT.ALLOCATED(XKER_RACCS ) ) ALLOCATE( XKER_RACCS (NACCLBDAS,NACCLBDAR) )
720 IF( .NOT.ALLOCATED(XKER_SACCRG) ) ALLOCATE( XKER_SACCRG(NACCLBDAR,NACCLBDAS) )
721 !
722 CALL READ_XKER_RACCS (KACCLBDAS,KACCLBDAR,KND,                                &
723                       PALPHAS,PNUS,PALPHAR,PNUR,PESR,PBS,PBR,PCS,PDS,PCR,PDR, &
724                       PACCLBDAS_MAX,PACCLBDAR_MAX,PACCLBDAS_MIN,PACCLBDAR_MIN,&
725                       PFDINFTY                                                )
726 IF( (KACCLBDAS/=NACCLBDAS) .OR. (KACCLBDAR/=NACCLBDAR) .OR. (KND/=IND) .OR. &
727     (PALPHAS/=XALPHAS) .OR. (PNUS/=XNUS)                               .OR. &
728     (PALPHAR/=XALPHAR) .OR. (PNUR/=XNUR)                               .OR. &
729     (PESR/=ZESR) .OR. (PBS/=XBS) .OR. (PBR/=XBR)                       .OR. &
730     (PCS/=XCS) .OR. (PDS/=XDS) .OR. (PCR/=XCR) .OR. (PDR/=XDR)         .OR. &
731     (PACCLBDAS_MAX/=XACCLBDAS_MAX) .OR. (PACCLBDAR_MAX/=XACCLBDAR_MAX) .OR. &
732     (PACCLBDAS_MIN/=XACCLBDAS_MIN) .OR. (PACCLBDAR_MIN/=XACCLBDAR_MIN) .OR. &
733     (PFDINFTY/=ZFDINFTY)                                               ) THEN
734   CALL RRCOLSS ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
735                  ZESR, XBR, XCS, XDS, XCR, XDR,                              &
736                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
737                  ZFDINFTY, XKER_RACCSS, XAG, XBS, XAS                        )
738   CALL RZCOLX  ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
739                  ZESR, XBR, XCS, XDS, XCR, XDR,                              &
740                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
741                  ZFDINFTY, XKER_RACCS                                        )
742   CALL RSCOLRG ( IND, XALPHAS, XNUS, XALPHAR, XNUR,                          &
743                  ZESR, XBS, XCS, XDS, XCR, XDR,                              &
744                  XACCLBDAS_MAX, XACCLBDAR_MAX, XACCLBDAS_MIN, XACCLBDAR_MIN, &
745                  ZFDINFTY, XKER_SACCRG,  XAG, XBS, XAS                       )
746   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
747   WRITE(UNIT=KLUOUT,FMT='("**** UPDATE NEW SET OF RACSS KERNELS ****")')
748   WRITE(UNIT=KLUOUT,FMT='("**** UPDATE NEW SET OF RACS  KERNELS ****")')
749   WRITE(UNIT=KLUOUT,FMT='("**** UPDATE NEW SET OF SACRG KERNELS ****")')
750   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
751   WRITE(UNIT=KLUOUT,FMT='("!")')
752   WRITE(UNIT=KLUOUT,FMT='("KND=",I3)') IND
753   WRITE(UNIT=KLUOUT,FMT='("KACCLBDAS=",I3)') NACCLBDAS
754   WRITE(UNIT=KLUOUT,FMT='("KACCLBDAR=",I3)') NACCLBDAR
755   WRITE(UNIT=KLUOUT,FMT='("PALPHAS=",E13.6)') XALPHAS
756   WRITE(UNIT=KLUOUT,FMT='("PNUS=",E13.6)') XNUS
757   WRITE(UNIT=KLUOUT,FMT='("PALPHAR=",E13.6)') XALPHAR
758   WRITE(UNIT=KLUOUT,FMT='("PNUR=",E13.6)') XNUR
759   WRITE(UNIT=KLUOUT,FMT='("PESR=",E13.6)') ZESR
760   WRITE(UNIT=KLUOUT,FMT='("PBS=",E13.6)') XBS
761   WRITE(UNIT=KLUOUT,FMT='("PBR=",E13.6)') XBR
762   WRITE(UNIT=KLUOUT,FMT='("PCS=",E13.6)') XCS
763   WRITE(UNIT=KLUOUT,FMT='("PDS=",E13.6)') XDS
764   WRITE(UNIT=KLUOUT,FMT='("PCR=",E13.6)') XCR
765   WRITE(UNIT=KLUOUT,FMT='("PDR=",E13.6)') XDR
766   WRITE(UNIT=KLUOUT,FMT='("PACCLBDAS_MAX=",E13.6)') &
767                                                     XACCLBDAS_MAX
768   WRITE(UNIT=KLUOUT,FMT='("PACCLBDAR_MAX=",E13.6)') &
769                                                     XACCLBDAR_MAX
770   WRITE(UNIT=KLUOUT,FMT='("PACCLBDAS_MIN=",E13.6)') &
771                                                     XACCLBDAS_MIN
772   WRITE(UNIT=KLUOUT,FMT='("PACCLBDAR_MIN=",E13.6)') &
773                                                     XACCLBDAR_MIN
774   WRITE(UNIT=KLUOUT,FMT='("PFDINFTY=",E13.6)') ZFDINFTY
775   WRITE(UNIT=KLUOUT,FMT='("!")')
776   WRITE(UNIT=KLUOUT,FMT='("IF( PRESENT(PKER_RACCSS) ) THEN")')
777   DO J1 = 1 , NACCLBDAS
778     DO J2 = 1 , NACCLBDAR
779     WRITE(UNIT=KLUOUT,FMT='("  PKER_RACCSS(",I3,",",I3,") = ",E13.6)') &
780                         J1,J2,XKER_RACCSS(J1,J2)
781     END DO
782   END DO
783   WRITE(UNIT=KLUOUT,FMT='("END IF")')
784   WRITE(UNIT=KLUOUT,FMT='("!")')
785   WRITE(UNIT=KLUOUT,FMT='("IF( PRESENT(PKER_RACCS ) ) THEN")')
786   DO J1 = 1 , NACCLBDAS
787     DO J2 = 1 , NACCLBDAR
788     WRITE(UNIT=KLUOUT,FMT='("  PKER_RACCS (",I3,",",I3,") = ",E13.6)') &
789                         J1,J2,XKER_RACCS (J1,J2)
790     END DO
791   END DO
792   WRITE(UNIT=KLUOUT,FMT='("END IF")')
793   WRITE(UNIT=KLUOUT,FMT='("!")')
794   WRITE(UNIT=KLUOUT,FMT='("IF( PRESENT(PKER_SACCRG) ) THEN")')
795   DO J1 = 1 , NACCLBDAR
796     DO J2 = 1 , NACCLBDAS
797     WRITE(UNIT=KLUOUT,FMT='("  PKER_SACCRG(",I3,",",I3,") = ",E13.6)') &
798                         J1,J2,XKER_SACCRG(J1,J2)
799     END DO
800   END DO
801   WRITE(UNIT=KLUOUT,FMT='("END IF")')
802   ELSE
803   CALL READ_XKER_RACCS (KACCLBDAS,KACCLBDAR,KND,                               &
804                        PALPHAS,PNUS,PALPHAR,PNUR,PESR,PBS,PBR,PCS,PDS,PCR,PDR, &
805                        PACCLBDAS_MAX,PACCLBDAR_MAX,PACCLBDAS_MIN,PACCLBDAR_MIN,&
806                        PFDINFTY,XKER_RACCSS,XKER_RACCS,XKER_SACCRG             )
807   WRITE(UNIT=KLUOUT,FMT='(" Read XKER_RACCSS")')
808   WRITE(UNIT=KLUOUT,FMT='(" Read XKER_RACCS ")')
809   WRITE(UNIT=KLUOUT,FMT='(" Read XKER_SACCRG")')
810 END IF
811 !
812 !*       7.3    Constant for the conversion-melting rate
813 !
814 XFSCVMG = 2.0
815 !
816 GFLAG = .TRUE.
817 IF (GFLAG) THEN
818   WRITE(UNIT=KLUOUT,FMT='("      conversion-melting of the aggregates")')
819   WRITE(UNIT=KLUOUT,FMT='(" Conv. factor XFSCVMG=",E13.6)') XFSCVMG
820 END IF
821 !
822 !
823 !-------------------------------------------------------------------------------
824 !
825 !*       8.     CONSTANTS FOR THE FAST COLD PROCESSES FOR THE GRAUPELN
826 !               ------------------------------------------------------
827 !
828 !
829 !*       8.1    Constants for the rain contact freezing
830 !
831 XCOLIR    = 1.0
832 !
833 XEXRCFRI  = -XDR-5.0+ZXR
834 XRCFRI    = ((XPI**2)/24.0)*XCCR*XRHOLW*XCOLIR*XCR*(ZRHO00**XCEXVT)     &
835                                                      *MOMG(XALPHAR,XNUR,XDR+5.0)
836 XEXICFRR  = -XDR-2.0+ZXR
837 XICFRR    = (XPI/4.0)*XCOLIR*XCR*(ZRHO00**XCEXVT)          &
838                                    *XCCR*MOMG(XALPHAR,XNUR,XDR+2.0)
839 !
840 GFLAG = .TRUE.
841 IF (GFLAG) THEN
842   WRITE(UNIT=KLUOUT,FMT='("      rain contact freezing")')
843   WRITE(UNIT=KLUOUT,FMT='(" Coll. efficiency          XCOLIR=",E13.6)') XCOLIR
844 END IF
845 !
846 !
847 !*       8.2    Constants for the dry growth of the graupeln
848 !
849 !*       8.2.1  Constants for the cloud droplet collection by the graupeln
850 !
851 XFCDRYG = (XPI/4.0)*XCCG*XCG*(ZRHO00**XCEXVT)*MOMG(XALPHAG,XNUG,XDG+2.0)
852 !
853 !*       8.2.2  Constants for the cloud ice collection by the graupeln
854 !
855 XCOLIG    = 0.25 ! Collection efficiency of I+G
856 XCOLEXIG  = 0.05 ! Temperature factor of the I+G collection efficiency
857 XCOLIG   = 0.01 ! Collection efficiency of I+G
858 XCOLEXIG = 0.1  ! Temperature factor of the I+G collection efficiency
859 WRITE (KLUOUT, FMT=*) ' NEW Constants for the cloud ice collection by the graupeln'
860 WRITE (KLUOUT, FMT=*) ' XCOLIG, XCOLEXIG  = ',XCOLIG,XCOLEXIG
861 XFIDRYG = (XPI/4.0)*XCOLIG*XCCG*XCG*(ZRHO00**XCEXVT)*MOMG(XALPHAG,XNUG,XDG+2.0)
862 !
863 GFLAG = .TRUE.
864 IF (GFLAG) THEN
865   WRITE(UNIT=KLUOUT,FMT='("      cloud ice collection by the graupeln")')
866   WRITE(UNIT=KLUOUT,FMT='(" Coll. efficiency XCOLIG=",E13.6)') XCOLIG
867   WRITE(UNIT=KLUOUT,FMT='(" Temp. factor     XCOLEXIG=",E13.6)') XCOLEXIG
868 END IF
869 !
870 !*       8.2.3  Constants for the aggregate collection by the graupeln
871 !
872 XCOLSG    = 0.25 ! Collection efficiency of S+G
873 XCOLEXSG  = 0.05 ! Temperature factor of the S+G collection efficiency
874 XCOLSG   = 0.01 ! Collection efficiency of S+G
875 XCOLEXSG = 0.1  ! Temperature factor of the S+G collection efficiency
876 WRITE (KLUOUT, FMT=*) ' NEW Constants for the aggregate collection by the graupeln'
877 WRITE (KLUOUT, FMT=*) ' XCOLSG, XCOLEXSG  = ',XCOLSG,XCOLEXSG
878 XFSDRYG = (XPI/4.0)*XCOLSG*XCCG*XCCS*XAS*(ZRHO00**XCEXVT)
879 !
880 XLBSDRYG1   =    MOMG(XALPHAG,XNUG,2.)*MOMG(XALPHAS,XNUS,XBS)
881 XLBSDRYG2   = 2.*MOMG(XALPHAG,XNUG,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
882 XLBSDRYG3   =                          MOMG(XALPHAS,XNUS,XBS+2.)
883 !
884 GFLAG = .TRUE.
885 IF (GFLAG) THEN
886   WRITE(UNIT=KLUOUT,FMT='("      aggregate collection by the graupeln")')
887   WRITE(UNIT=KLUOUT,FMT='(" Coll. efficiency XCOLSG=",E13.6)') XCOLSG
888   WRITE(UNIT=KLUOUT,FMT='(" Temp. factor     XCOLEXSG=",E13.6)') XCOLEXSG
889 END IF
890 !
891 !*       8.2.4  Constants for the raindrop collection by the graupeln
892 !
893 XFRDRYG = ((XPI**2)/24.0)*XCCG*XCCR*XRHOLW*(ZRHO00**XCEXVT)
894 !
895 XLBRDRYG1   =    MOMG(XALPHAG,XNUG,2.)*MOMG(XALPHAR,XNUR,3.)
896 XLBRDRYG2   = 2.*MOMG(XALPHAG,XNUG,1.)*MOMG(XALPHAR,XNUR,4.)
897 XLBRDRYG3   =                          MOMG(XALPHAR,XNUR,5.)
898 !
899 ! Notice: One magnitude of lambda discretized over 10 points
900 !
901 NDRYLBDAR = 40
902 XDRYLBDAR_MIN = 1.0E3 ! Minimal value of Lbda_r to tabulate XKER_RDRYG
903 XDRYLBDAR_MAX = 1.0E7 ! Maximal value of Lbda_r to tabulate XKER_RDRYG
904 ZRATE = LOG(XDRYLBDAR_MAX/XDRYLBDAR_MIN)/FLOAT(NDRYLBDAR-1)
905 XDRYINTP1R = 1.0 / ZRATE
906 XDRYINTP2R = 1.0 - LOG( XDRYLBDAR_MIN ) / ZRATE
907 NDRYLBDAS = 80
908 XDRYLBDAS_MIN = 2.5E1 ! Minimal value of Lbda_s to tabulate XKER_SDRYG
909 XDRYLBDAS_MAX = 2.5E9 ! Maximal value of Lbda_s to tabulate XKER_SDRYG
910 ZRATE = LOG(XDRYLBDAS_MAX/XDRYLBDAS_MIN)/FLOAT(NDRYLBDAS-1)
911 XDRYINTP1S = 1.0 / ZRATE
912 XDRYINTP2S = 1.0 - LOG( XDRYLBDAS_MIN ) / ZRATE
913 NDRYLBDAG = 40
914 XDRYLBDAG_MIN = 1.0E3 ! Min value of Lbda_g to tabulate XKER_SDRYG,XKER_RDRYG
915 XDRYLBDAG_MAX = 1.0E7 ! Max value of Lbda_g to tabulate XKER_SDRYG,XKER_RDRYG
916 ZRATE = LOG(XDRYLBDAG_MAX/XDRYLBDAG_MIN)/FLOAT(NDRYLBDAG-1)
917 XDRYINTP1G = 1.0 / ZRATE
918 XDRYINTP2G = 1.0 - LOG( XDRYLBDAG_MIN ) / ZRATE
919 !
920 !*       8.2.5  Computations of the tabulated normalized kernels
921 !
922 IND      = 50    ! Interval number, collection efficiency and infinite diameter
923 ZEGS     = 1.0   ! factor used to integrate the dimensional distributions when
924 ZFDINFTY = 20.0  ! computing the kernels XKER_SDRYG
925 !
926 IF( .NOT.ALLOCATED(XKER_SDRYG) ) ALLOCATE( XKER_SDRYG(NDRYLBDAG,NDRYLBDAS) )
927 !
928 CALL READ_XKER_SDRYG (KDRYLBDAG,KDRYLBDAS,KND,                              &
929                    PALPHAG,PNUG,PALPHAS,PNUS,PEGS,PBS,PCG,PDG,PCS,PDS,      &
930                    PDRYLBDAG_MAX,PDRYLBDAS_MAX,PDRYLBDAG_MIN,PDRYLBDAS_MIN, &
931                    PFDINFTY                                                 )
932 IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAS/=NDRYLBDAS) .OR. (KND/=IND) .OR. &
933     (PALPHAG/=XALPHAG) .OR. (PNUG/=XNUG)                               .OR. &
934     (PALPHAS/=XALPHAS) .OR. (PNUS/=XNUS)                               .OR. &
935     (PEGS/=ZEGS) .OR. (PBS/=XBS)                                       .OR. &
936     (PCG/=XCG) .OR. (PDG/=XDG) .OR. (PCS/=XCS) .OR. (PDS/=XDS)         .OR. &
937     (PDRYLBDAG_MAX/=XDRYLBDAG_MAX) .OR. (PDRYLBDAS_MAX/=XDRYLBDAS_MAX) .OR. &
938     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAS_MIN/=XDRYLBDAS_MIN) .OR. &
939     (PFDINFTY/=ZFDINFTY)                                               ) THEN
940   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAS, XNUS,                          &
941                 ZEGS, XBS, XCG, XDG, XCS, XDS,                              &
942                 XDRYLBDAG_MAX, XDRYLBDAS_MAX, XDRYLBDAG_MIN, XDRYLBDAS_MIN, &
943                 ZFDINFTY, XKER_SDRYG                                        )
944   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
945   WRITE(UNIT=KLUOUT,FMT='("**** UPDATE NEW SET OF SDRYG KERNELS ****")')
946   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
947   WRITE(UNIT=KLUOUT,FMT='("!")')
948   WRITE(UNIT=KLUOUT,FMT='("KND=",I3)') IND
949   WRITE(UNIT=KLUOUT,FMT='("KDRYLBDAG=",I3)') NDRYLBDAG
950   WRITE(UNIT=KLUOUT,FMT='("KDRYLBDAS=",I3)') NDRYLBDAS
951   WRITE(UNIT=KLUOUT,FMT='("PALPHAG=",E13.6)') XALPHAG
952   WRITE(UNIT=KLUOUT,FMT='("PNUG=",E13.6)') XNUG
953   WRITE(UNIT=KLUOUT,FMT='("PALPHAS=",E13.6)') XALPHAS
954   WRITE(UNIT=KLUOUT,FMT='("PNUS=",E13.6)') XNUS
955   WRITE(UNIT=KLUOUT,FMT='("PEGS=",E13.6)') ZEGS
956   WRITE(UNIT=KLUOUT,FMT='("PBS=",E13.6)') XBS
957   WRITE(UNIT=KLUOUT,FMT='("PCG=",E13.6)') XCG
958   WRITE(UNIT=KLUOUT,FMT='("PDG=",E13.6)') XDG
959   WRITE(UNIT=KLUOUT,FMT='("PCS=",E13.6)') XCS
960   WRITE(UNIT=KLUOUT,FMT='("PDS=",E13.6)') XDS
961   WRITE(UNIT=KLUOUT,FMT='("PDRYLBDAG_MAX=",E13.6)') &
962                                                     XDRYLBDAG_MAX
963   WRITE(UNIT=KLUOUT,FMT='("PDRYLBDAS_MAX=",E13.6)') &
964                                                     XDRYLBDAS_MAX
965   WRITE(UNIT=KLUOUT,FMT='("PDRYLBDAG_MIN=",E13.6)') &
966                                                     XDRYLBDAG_MIN
967   WRITE(UNIT=KLUOUT,FMT='("PDRYLBDAS_MIN=",E13.6)') &
968                                                     XDRYLBDAS_MIN
969   WRITE(UNIT=KLUOUT,FMT='("PFDINFTY=",E13.6)') ZFDINFTY
970   WRITE(UNIT=KLUOUT,FMT='("!")')
971   WRITE(UNIT=KLUOUT,FMT='("IF( PRESENT(PKER_SDRYG) ) THEN")')
972   DO J1 = 1 , NDRYLBDAG
973     DO J2 = 1 , NDRYLBDAS
974     WRITE(UNIT=KLUOUT,FMT='("PKER_SDRYG(",I3,",",I3,") = ",E13.6)') &
975                         J1,J2,XKER_SDRYG(J1,J2)
976     END DO
977   END DO
978   WRITE(UNIT=KLUOUT,FMT='("END IF")')
979   ELSE
980   CALL READ_XKER_SDRYG (KDRYLBDAG,KDRYLBDAS,KND,                              &
981                      PALPHAG,PNUG,PALPHAS,PNUS,PEGS,PBS,PCG,PDG,PCS,PDS,      &
982                      PDRYLBDAG_MAX,PDRYLBDAS_MAX,PDRYLBDAG_MIN,PDRYLBDAS_MIN, &
983                      PFDINFTY,XKER_SDRYG                                      )
984   WRITE(UNIT=KLUOUT,FMT='(" Read XKER_SDRYG")')
985 END IF
986 !
987 !
988 IND      = 50    ! Number of interval used to integrate the dimensional
989 ZEGR     = 1.0   ! distributions when computing the kernel XKER_RDRYG
990 ZFDINFTY = 20.0
991 !
992 IF( .NOT.ALLOCATED(XKER_RDRYG) ) ALLOCATE( XKER_RDRYG(NDRYLBDAG,NDRYLBDAR) )
993 !
994 CALL READ_XKER_RDRYG (KDRYLBDAG,KDRYLBDAR,KND,                              &
995                    PALPHAG,PNUG,PALPHAR,PNUR,PEGR,PBR,PCG,PDG,PCR,PDR,      &
996                    PDRYLBDAG_MAX,PDRYLBDAR_MAX,PDRYLBDAG_MIN,PDRYLBDAR_MIN, &
997                    PFDINFTY                                                 )
998 IF( (KDRYLBDAG/=NDRYLBDAG) .OR. (KDRYLBDAR/=NDRYLBDAR) .OR. (KND/=IND) .OR. &
999     (PALPHAG/=XALPHAG) .OR. (PNUG/=XNUG)                               .OR. &
1000     (PALPHAR/=XALPHAR) .OR. (PNUR/=XNUR)                               .OR. &
1001     (PEGR/=ZEGR) .OR. (PBR/=XBR)                                       .OR. &
1002     (PCG/=XCG) .OR. (PDG/=XDG) .OR. (PCR/=XCR) .OR. (PDR/=XDR)         .OR. &
1003     (PDRYLBDAG_MAX/=XDRYLBDAG_MAX) .OR. (PDRYLBDAR_MAX/=XDRYLBDAR_MAX) .OR. &
1004     (PDRYLBDAG_MIN/=XDRYLBDAG_MIN) .OR. (PDRYLBDAR_MIN/=XDRYLBDAR_MIN) .OR. &
1005     (PFDINFTY/=ZFDINFTY)                                               ) THEN
1006   CALL RZCOLX ( IND, XALPHAG, XNUG, XALPHAR, XNUR,                          &
1007                 ZEGR, XBR, XCG, XDG, XCR, XDR,                              &
1008                 XDRYLBDAG_MAX, XDRYLBDAR_MAX, XDRYLBDAG_MIN, XDRYLBDAR_MIN, &
1009                 ZFDINFTY, XKER_RDRYG                                        )
1010   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
1011   WRITE(UNIT=KLUOUT,FMT='("**** UPDATE NEW SET OF RDRYG KERNELS ****")')
1012   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
1013   WRITE(UNIT=KLUOUT,FMT='("!")')
1014   WRITE(UNIT=KLUOUT,FMT='("KND=",I3)') IND
1015   WRITE(UNIT=KLUOUT,FMT='("KDRYLBDAG=",I3)') NDRYLBDAG
1016   WRITE(UNIT=KLUOUT,FMT='("KDRYLBDAR=",I3)') NDRYLBDAR
1017   WRITE(UNIT=KLUOUT,FMT='("PALPHAG=",E13.6)') XALPHAG
1018   WRITE(UNIT=KLUOUT,FMT='("PNUG=",E13.6)') XNUG
1019   WRITE(UNIT=KLUOUT,FMT='("PALPHAR=",E13.6)') XALPHAR
1020   WRITE(UNIT=KLUOUT,FMT='("PNUR=",E13.6)') XNUR
1021   WRITE(UNIT=KLUOUT,FMT='("PEGR=",E13.6)') ZEGR
1022   WRITE(UNIT=KLUOUT,FMT='("PBR=",E13.6)') XBR
1023   WRITE(UNIT=KLUOUT,FMT='("PCG=",E13.6)') XCG
1024   WRITE(UNIT=KLUOUT,FMT='("PDG=",E13.6)') XDG
1025   WRITE(UNIT=KLUOUT,FMT='("PCR=",E13.6)') XCR
1026   WRITE(UNIT=KLUOUT,FMT='("PDR=",E13.6)') XDR
1027   WRITE(UNIT=KLUOUT,FMT='("PDRYLBDAG_MAX=",E13.6)') &
1028                                                     XDRYLBDAG_MAX
1029   WRITE(UNIT=KLUOUT,FMT='("PDRYLBDAR_MAX=",E13.6)') &
1030                                                     XDRYLBDAR_MAX
1031   WRITE(UNIT=KLUOUT,FMT='("PDRYLBDAG_MIN=",E13.6)') &
1032                                                     XDRYLBDAG_MIN
1033   WRITE(UNIT=KLUOUT,FMT='("PDRYLBDAR_MIN=",E13.6)') &
1034                                                     XDRYLBDAR_MIN
1035   WRITE(UNIT=KLUOUT,FMT='("PFDINFTY=",E13.6)') ZFDINFTY
1036   WRITE(UNIT=KLUOUT,FMT='("!")')
1037   WRITE(UNIT=KLUOUT,FMT='("IF( PRESENT(PKER_RDRYG) ) THEN")')
1038   DO J1 = 1 , NDRYLBDAG
1039     DO J2 = 1 , NDRYLBDAR
1040     WRITE(UNIT=KLUOUT,FMT='("PKER_RDRYG(",I3,",",I3,") = ",E13.6)') &
1041                         J1,J2,XKER_RDRYG(J1,J2)
1042     END DO
1043   END DO
1044   WRITE(UNIT=KLUOUT,FMT='("END IF")')
1045   ELSE
1046   CALL READ_XKER_RDRYG (KDRYLBDAG,KDRYLBDAR,KND,                              &
1047                      PALPHAG,PNUG,PALPHAR,PNUR,PEGR,PBR,PCG,PDG,PCR,PDR,      &
1048                      PDRYLBDAG_MAX,PDRYLBDAR_MAX,PDRYLBDAG_MIN,PDRYLBDAR_MIN, &
1049                      PFDINFTY,XKER_RDRYG                                      )
1050   WRITE(UNIT=KLUOUT,FMT='(" Read XKER_RDRYG")')
1051 END IF
1052 !
1053 !
1054 !-------------------------------------------------------------------------------
1055 !
1056 !*       9.     CONSTANTS FOR THE FAST COLD PROCESSES FOR THE HAILSTONES
1057 !               --------------------------------------------------------
1058 !
1059 !*       9.2    Constants for the wet growth of the hailstones
1060 !
1061 !
1062 !*       9.2.1  Constant for the cloud droplet and cloud ice collection
1063 !               by the hailstones
1064 !
1065 XFWETH = (XPI/4.0)*XCCH*XCH*(ZRHO00**XCEXVT)*MOMG(XALPHAH,XNUH,XDH+2.0)
1066 !
1067 !*       9.2.2  Constants for the aggregate collection by the hailstones
1068 !
1069 XFSWETH = (XPI/4.0)*XCCH*XCCS*XAS*(ZRHO00**XCEXVT)
1070 !
1071 XLBSWETH1   =    MOMG(XALPHAH,XNUH,2.)*MOMG(XALPHAS,XNUS,XBS)
1072 XLBSWETH2   = 2.*MOMG(XALPHAH,XNUH,1.)*MOMG(XALPHAS,XNUS,XBS+1.)
1073 XLBSWETH3   =                          MOMG(XALPHAS,XNUS,XBS+2.)
1074 !
1075 !*       9.2.3  Constants for the graupel collection by the hailstones
1076 !
1077 XFGWETH = (XPI/4.0)*XCCH*XCCG*XAG*(ZRHO00**XCEXVT)
1078 !
1079 XLBGWETH1   =    MOMG(XALPHAH,XNUH,2.)*MOMG(XALPHAG,XNUG,XBG)
1080 XLBGWETH2   = 2.*MOMG(XALPHAH,XNUH,1.)*MOMG(XALPHAG,XNUG,XBG+1.)
1081 XLBGWETH3   =                          MOMG(XALPHAG,XNUG,XBG+2.)
1082 !
1083 ! Notice: One magnitude of lambda discretized over 10 points
1084 !
1085 NWETLBDAS = 80
1086 XWETLBDAS_MIN = 2.5E1 ! Minimal value of Lbda_s to tabulate XKER_SWETH
1087 XWETLBDAS_MAX = 2.5E9 ! Maximal value of Lbda_s to tabulate XKER_SWETH
1088 ZRATE = LOG(XWETLBDAS_MAX/XWETLBDAS_MIN)/FLOAT(NWETLBDAS-1)
1089 XWETINTP1S = 1.0 / ZRATE
1090 XWETINTP2S = 1.0 - LOG( XWETLBDAS_MIN ) / ZRATE
1091 NWETLBDAG = 40
1092 XWETLBDAG_MIN = 1.0E3 ! Min value of Lbda_g to tabulate XKER_GWETH
1093 XWETLBDAG_MAX = 1.0E7 ! Max value of Lbda_g to tabulate XKER_GWETH
1094 ZRATE = LOG(XWETLBDAG_MAX/XWETLBDAG_MIN)/FLOAT(NWETLBDAG-1)
1095 XWETINTP1G = 1.0 / ZRATE
1096 XWETINTP2G = 1.0 - LOG( XWETLBDAG_MIN ) / ZRATE
1097 NWETLBDAH = 40
1098 XWETLBDAH_MIN = 1.0E3 ! Min value of Lbda_h to tabulate XKER_SWETH,XKER_GWETH
1099 XWETLBDAH_MAX = 1.0E7 ! Max value of Lbda_h to tabulate XKER_SWETH,XKER_GWETH
1100 ZRATE = LOG(XWETLBDAH_MAX/XWETLBDAH_MIN)/FLOAT(NWETLBDAH-1)
1101 XWETINTP1H = 1.0 / ZRATE
1102 XWETINTP2H = 1.0 - LOG( XWETLBDAH_MIN ) / ZRATE
1103 !
1104 !*       9.2.4  Computations of the tabulated normalized kernels
1105 !
1106 IND      = 50    ! Interval number, collection efficiency and infinite diameter
1107 ZEHS     = 1.0   ! factor used to integrate the dimensional distributions when
1108 ZFDINFTY = 20.0  ! computing the kernels XKER_SWETH
1109 !
1110 IF( .NOT.ALLOCATED(XKER_SWETH) ) ALLOCATE( XKER_SWETH(NWETLBDAH,NWETLBDAS) )
1111 !
1112 CALL READ_XKER_SWETH (KWETLBDAH,KWETLBDAS,KND,                              &
1113                    PALPHAH,PNUH,PALPHAS,PNUS,PEHS,PBS,PCH,PDH,PCS,PDS,      &
1114                    PWETLBDAH_MAX,PWETLBDAS_MAX,PWETLBDAH_MIN,PWETLBDAS_MIN, &
1115                    PFDINFTY                                                 )
1116 IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAS/=NWETLBDAS) .OR. (KND/=IND) .OR. &
1117     (PALPHAH/=XALPHAH) .OR. (PNUH/=XNUH)                               .OR. &
1118     (PALPHAS/=XALPHAS) .OR. (PNUS/=XNUS)                               .OR. &
1119     (PEHS/=ZEHS) .OR. (PBS/=XBS)                                       .OR. &
1120     (PCH/=XCH) .OR. (PDH/=XDH) .OR. (PCS/=XCS) .OR. (PDS/=XDS)         .OR. &
1121     (PWETLBDAH_MAX/=XWETLBDAH_MAX) .OR. (PWETLBDAS_MAX/=XWETLBDAS_MAX) .OR. &
1122     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAS_MIN/=XWETLBDAS_MIN) .OR. &
1123     (PFDINFTY/=ZFDINFTY)                                               ) THEN
1124   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAS, XNUS,                          &
1125                 ZEHS, XBS, XCH, XDH, XCS, XDS,                              &
1126                 XWETLBDAH_MAX, XWETLBDAS_MAX, XWETLBDAH_MIN, XWETLBDAS_MIN, &
1127                 ZFDINFTY, XKER_SWETH                                        )
1128   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
1129   WRITE(UNIT=KLUOUT,FMT='("**** UPDATE NEW SET OF SWETH KERNELS ****")')
1130   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
1131   WRITE(UNIT=KLUOUT,FMT='("!")')
1132   WRITE(UNIT=KLUOUT,FMT='("KND=",I3)') IND
1133   WRITE(UNIT=KLUOUT,FMT='("KWETLBDAH=",I3)') NWETLBDAH
1134   WRITE(UNIT=KLUOUT,FMT='("KWETLBDAS=",I3)') NWETLBDAS
1135   WRITE(UNIT=KLUOUT,FMT='("PALPHAH=",E13.6)') XALPHAH
1136   WRITE(UNIT=KLUOUT,FMT='("PNUH=",E13.6)') XNUH
1137   WRITE(UNIT=KLUOUT,FMT='("PALPHAS=",E13.6)') XALPHAS
1138   WRITE(UNIT=KLUOUT,FMT='("PNUS=",E13.6)') XNUS
1139   WRITE(UNIT=KLUOUT,FMT='("PEHS=",E13.6)') ZEHS
1140   WRITE(UNIT=KLUOUT,FMT='("PBS=",E13.6)') XBS
1141   WRITE(UNIT=KLUOUT,FMT='("PCH=",E13.6)') XCH
1142   WRITE(UNIT=KLUOUT,FMT='("PDH=",E13.6)') XDH
1143   WRITE(UNIT=KLUOUT,FMT='("PCS=",E13.6)') XCS
1144   WRITE(UNIT=KLUOUT,FMT='("PDS=",E13.6)') XDS
1145   WRITE(UNIT=KLUOUT,FMT='("PWETLBDAH_MAX=",E13.6)') &
1146                                                     XWETLBDAH_MAX
1147   WRITE(UNIT=KLUOUT,FMT='("PWETLBDAS_MAX=",E13.6)') &
1148                                                     XWETLBDAS_MAX
1149   WRITE(UNIT=KLUOUT,FMT='("PWETLBDAH_MIN=",E13.6)') &
1150                                                     XWETLBDAH_MIN
1151   WRITE(UNIT=KLUOUT,FMT='("PWETLBDAS_MIN=",E13.6)') &
1152                                                     XWETLBDAS_MIN
1153   WRITE(UNIT=KLUOUT,FMT='("PFDINFTY=",E13.6)') ZFDINFTY
1154   WRITE(UNIT=KLUOUT,FMT='("!")')
1155   WRITE(UNIT=KLUOUT,FMT='("IF( PRESENT(PKER_SWETH) ) THEN")')
1156   DO J1 = 1 , NWETLBDAH
1157     DO J2 = 1 , NWETLBDAS
1158     WRITE(UNIT=KLUOUT,FMT='("PKER_SWETH(",I3,",",I3,") = ",E13.6)') &
1159                         J1,J2,XKER_SWETH(J1,J2)
1160     END DO
1161   END DO
1162   WRITE(UNIT=KLUOUT,FMT='("END IF")')
1163   ELSE
1164   CALL READ_XKER_SWETH (KWETLBDAH,KWETLBDAS,KND,                              &
1165                      PALPHAH,PNUH,PALPHAS,PNUS,PEHS,PBS,PCH,PDH,PCS,PDS,      &
1166                      PWETLBDAH_MAX,PWETLBDAS_MAX,PWETLBDAH_MIN,PWETLBDAS_MIN, &
1167                      PFDINFTY,XKER_SWETH                                      )
1168   WRITE(UNIT=KLUOUT,FMT='(" Read XKER_SWETH")')
1169 END IF
1170 !
1171 !
1172 IND      = 50    ! Number of interval used to integrate the dimensional
1173 ZEHG     = 1.0   ! distributions when computing the kernel XKER_GWETH
1174 ZFDINFTY = 20.0
1175 !
1176 IF( .NOT.ALLOCATED(XKER_GWETH) ) ALLOCATE( XKER_GWETH(NWETLBDAH,NWETLBDAG) )
1177 !
1178 CALL READ_XKER_GWETH (KWETLBDAH,KWETLBDAG,KND,                              &
1179                    PALPHAH,PNUH,PALPHAG,PNUG,PEHG,PBG,PCH,PDH,PCG,PDG,      &
1180                    PWETLBDAH_MAX,PWETLBDAG_MAX,PWETLBDAH_MIN,PWETLBDAG_MIN, &
1181                    PFDINFTY                                                 )
1182 IF( (KWETLBDAH/=NWETLBDAH) .OR. (KWETLBDAG/=NWETLBDAG) .OR. (KND/=IND) .OR. &
1183     (PALPHAH/=XALPHAH) .OR. (PNUH/=XNUH)                               .OR. &
1184     (PALPHAG/=XALPHAG) .OR. (PNUG/=XNUG)                               .OR. &
1185     (PEHG/=ZEHG) .OR. (PBG/=XBG)                                       .OR. &
1186     (PCH/=XCH) .OR. (PDH/=XDH) .OR. (PCG/=XCG) .OR. (PDG/=XDG)         .OR. &
1187     (PWETLBDAH_MAX/=XWETLBDAH_MAX) .OR. (PWETLBDAG_MAX/=XWETLBDAG_MAX) .OR. &
1188     (PWETLBDAH_MIN/=XWETLBDAH_MIN) .OR. (PWETLBDAG_MIN/=XWETLBDAG_MIN) .OR. &
1189     (PFDINFTY/=ZFDINFTY)                                               ) THEN
1190   CALL RZCOLX ( IND, XALPHAH, XNUH, XALPHAG, XNUG,                          &
1191                 ZEHG, XBG, XCH, XDH, XCG, XDG,                              &
1192                 XWETLBDAH_MAX, XWETLBDAG_MAX, XWETLBDAH_MIN, XWETLBDAG_MIN, &
1193                 ZFDINFTY, XKER_GWETH                                        )
1194   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
1195   WRITE(UNIT=KLUOUT,FMT='("**** UPDATE NEW SET OF GWETH KERNELS ****")')
1196   WRITE(UNIT=KLUOUT,FMT='("*****************************************")')
1197   WRITE(UNIT=KLUOUT,FMT='("!")')
1198   WRITE(UNIT=KLUOUT,FMT='("KND=",I3)') IND
1199   WRITE(UNIT=KLUOUT,FMT='("KWETLBDAH=",I3)') NWETLBDAH
1200   WRITE(UNIT=KLUOUT,FMT='("KWETLBDAG=",I3)') NWETLBDAG
1201   WRITE(UNIT=KLUOUT,FMT='("PALPHAH=",E13.6)') XALPHAH
1202   WRITE(UNIT=KLUOUT,FMT='("PNUH=",E13.6)') XNUH
1203   WRITE(UNIT=KLUOUT,FMT='("PALPHAG=",E13.6)') XALPHAG
1204   WRITE(UNIT=KLUOUT,FMT='("PNUG=",E13.6)') XNUG
1205   WRITE(UNIT=KLUOUT,FMT='("PEHG=",E13.6)') ZEHG
1206   WRITE(UNIT=KLUOUT,FMT='("PBG=",E13.6)') XBG
1207   WRITE(UNIT=KLUOUT,FMT='("PCH=",E13.6)') XCH
1208   WRITE(UNIT=KLUOUT,FMT='("PDH=",E13.6)') XDH
1209   WRITE(UNIT=KLUOUT,FMT='("PCG=",E13.6)') XCG
1210   WRITE(UNIT=KLUOUT,FMT='("PDG=",E13.6)') XDG
1211   WRITE(UNIT=KLUOUT,FMT='("PWETLBDAH_MAX=",E13.6)') &
1212                                                     XWETLBDAH_MAX
1213   WRITE(UNIT=KLUOUT,FMT='("PWETLBDAG_MAX=",E13.6)') &
1214                                                     XWETLBDAG_MAX
1215   WRITE(UNIT=KLUOUT,FMT='("PWETLBDAH_MIN=",E13.6)') &
1216                                                     XWETLBDAH_MIN
1217   WRITE(UNIT=KLUOUT,FMT='("PWETLBDAG_MIN=",E13.6)') &
1218                                                     XWETLBDAG_MIN
1219   WRITE(UNIT=KLUOUT,FMT='("PFDINFTY=",E13.6)') ZFDINFTY
1220   WRITE(UNIT=KLUOUT,FMT='("!")')
1221   WRITE(UNIT=KLUOUT,FMT='("IF( PRESENT(PKER_GWETH) ) THEN")')
1222   DO J1 = 1 , NWETLBDAH
1223     DO J2 = 1 , NWETLBDAG
1224     WRITE(UNIT=KLUOUT,FMT='("PKER_GWETH(",I3,",",I3,") = ",E13.6)') &
1225                         J1,J2,XKER_GWETH(J1,J2)
1226     END DO
1227   END DO
1228   WRITE(UNIT=KLUOUT,FMT='("END IF")')
1229   ELSE
1230   CALL READ_XKER_GWETH (KWETLBDAH,KWETLBDAG,KND,                              &
1231                      PALPHAH,PNUH,PALPHAG,PNUG,PEHG,PBG,PCH,PDH,PCG,PDG,      &
1232                      PWETLBDAH_MAX,PWETLBDAG_MAX,PWETLBDAH_MIN,PWETLBDAG_MIN, &
1233                      PFDINFTY,XKER_GWETH                                      )
1234   WRITE(UNIT=KLUOUT,FMT='(" Read XKER_GWETH")')
1235 END IF
1236 !
1237 !
1238 !-------------------------------------------------------------------------------
1239 !
1240 !*      10.     SOME PRINTS FOR CONTROL
1241 !               -----------------------
1242 !
1243 !
1244 GFLAG = .TRUE.
1245 IF (GFLAG) THEN
1246   WRITE(UNIT=KLUOUT,FMT='(" Summary of the ice particule characteristics")')
1247   WRITE(UNIT=KLUOUT,FMT='("      PRISTINE ICE")')
1248   WRITE(UNIT=KLUOUT,FMT='("                   masse: A=",E13.6," B=",E13.6)') &
1249                                                       XAI,XBI
1250   WRITE(UNIT=KLUOUT,FMT='("                 vitesse: C=",E13.6," D=",E13.6)') &
1251                                                       XC_I,XDI
1252   WRITE(UNIT=KLUOUT,FMT='("            distribution:AL=",E13.6,"NU=",E13.6)') &
1253                                                       XALPHAI,XNUI
1254   WRITE(UNIT=KLUOUT,FMT='("              SNOW")')
1255   WRITE(UNIT=KLUOUT,FMT='("                   masse: A=",E13.6," B=",E13.6)') &
1256                                                       XAS,XBS
1257   WRITE(UNIT=KLUOUT,FMT='("                 vitesse: C=",E13.6," D=",E13.6)') &
1258                                                       XCS,XDS
1259   WRITE(UNIT=KLUOUT,FMT='("           concentration:CC=",E13.6," x=",E13.6)') &
1260                                                       XCCS,XCXS
1261   WRITE(UNIT=KLUOUT,FMT='("            distribution:AL=",E13.6,"NU=",E13.6)') &
1262                                                       XALPHAS,XNUS
1263   WRITE(UNIT=KLUOUT,FMT='("            GRAUPEL")')
1264   WRITE(UNIT=KLUOUT,FMT='("                   masse: A=",E13.6," B=",E13.6)') &
1265                                                       XAG,XBG
1266   WRITE(UNIT=KLUOUT,FMT='("                 vitesse: C=",E13.6," D=",E13.6)') &
1267                                                       XCG,XDG
1268   WRITE(UNIT=KLUOUT,FMT='("           concentration:CC=",E13.6," x=",E13.6)') &
1269                                                       XCCG,XCXG
1270   WRITE(UNIT=KLUOUT,FMT='("            distribution:AL=",E13.6,"NU=",E13.6)') &
1271                                                       XALPHAG,XNUG
1272   WRITE(UNIT=KLUOUT,FMT='("               HAIL")')
1273   WRITE(UNIT=KLUOUT,FMT='("                   masse: A=",E13.6," B=",E13.6)') &
1274                                                       XAH,XBH
1275   WRITE(UNIT=KLUOUT,FMT='("                 vitesse: C=",E13.6," D=",E13.6)') &
1276                                                       XCH,XDH
1277   WRITE(UNIT=KLUOUT,FMT='("           concentration:CC=",E13.6," x=",E13.6)') &
1278                                                       XCCH,XCXH
1279   WRITE(UNIT=KLUOUT,FMT='("            distribution:AL=",E13.6,"NU=",E13.6)') &
1280                                                       XALPHAH,XNUH
1281 END IF
1282 CONTAINS
1283 !
1284 !------------------------------------------------------------------------------
1285 !
1286   FUNCTION MOMG(PALPHA,PNU,PP) RESULT (PMOMG)
1287 !
1288 ! auxiliary routine used to compute the Pth moment order of the generalized
1289 ! gamma law
1290 !
1291   USE MODI_GAMMA
1292 !
1293   IMPLICIT NONE
1294 !
1295   REAL     :: PALPHA ! first shape parameter of the dimensionnal distribution
1296   REAL     :: PNU    ! second shape parameter of the dimensionnal distribution
1297   REAL     :: PP     ! order of the moment
1298   REAL     :: PMOMG  ! result: moment of order ZP
1299 !
1300 !------------------------------------------------------------------------------
1301 !
1302 !
1303   PMOMG = GAMMA(PNU+PP/PALPHA)/GAMMA(PNU)
1304 !
1305   END FUNCTION MOMG
1306 !
1307 !-------------------------------------------------------------------------------
1308 !
1309 !
1310 END SUBROUTINE INI_RAIN_ICE