source: LMDZ6/trunk/libf/phylmd/lscp_mod.F90 @ 4393

Last change on this file since 4393 was 4392, checked in by evignon, 2 years ago

super session de commentaires pour preparer linclusion des isotopes
dans lscp. Camille, Cecile, Niels, Sebastien et Etienne
File size: 49.6 KB
Line 
1MODULE lscp_mod
2
3IMPLICIT NONE
4
5CONTAINS
6
7!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
9     paprs,pplay,t,q,ptconv,ratqs,                      &
10     d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri,  &
11     radliq, radicefrac, rain, snow,                    &
12     frac_impa, frac_nucl, beta,                        &
13     prfl, psfl, rhcl, zqta, fraca,                     &
14     ztv, zpspsk, ztla, zthl, iflag_cld_th,             &
15     iflag_ice_thermo, ok_ice_sursat, qsatl, qsats,     &
16     qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss,   &
17     Tcontr, qcontr, qcontr2, fcontrN, fcontrP)
18
19!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
21!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
22!--------------------------------------------------------------------------------
23! Date: 01/2021
24!--------------------------------------------------------------------------------
25! Aim: Large Scale Clouds and Precipitation (LSCP)
26!
27! This code is a new version of the fisrtilp.F90 routine, which is itself a
28! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
29! and 'ilp' (il pleut, L. Li)
30!
31! Compared to the original fisrtilp code, lscp
32! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
33! -> consider always precipitation thermalisation (fl_cor_ebil>0)
34! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
35! -> option "oldbug" by JYG has been removed
36! -> iflag_t_glace >0 always
37! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
38! -> rectangular distribution from L. Li no longer available
39! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
40!--------------------------------------------------------------------------------
41! References:
42!
43! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
44! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
45! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
46! -------------------------------------------------------------------------------
47! Code structure:
48!
49! P0>     Thermalization of the precipitation coming from the overlying layer
50! P1>     Evaporation of the precipitation (falling from the k+1 level)
51! P2>     Cloud formation (at the k level)
52! P2.A.0> Cloud properties calculation from a rectangular pdf
53! P2.A.1> With the new PDFs, calculation of cloud properties using the inital
54!         values of T and Q
55! P2.A.2> Coupling between condensed water and temperature
56! P2.A.3> Calculation of final quantities associated with cloud formation
57! P2.B>   'All or nothing' cloud
58! P2.C>   Release of Latent heat after cloud formation
59! P3>     Autoconversion to precipitation (k-level)
60! P4>     Wet scavenging
61!------------------------------------------------------------------------------
62! Some preliminary comments (JBM) :
63!
64! The cloud water that the radiation scheme sees is not the same that the cloud
65! water used in the physics and the dynamics
66!
67! During the autoconversion to precipitation (P3 step), radliq (cloud water used
68! by the radiation scheme) is calculated as an average of the water that remains
69! in the cloud during the precipitation and not the water remaining at the end
70! of the time step. The latter is used in the rest of the physics and advected
71! by the dynamics.
72!
73! In summary:
74!
75! Radiation:
76! xflwc(newmicro)+xfiwc(newmicro) =
77!   cldliq(physiq)=radliq(fisrt)=lwcon(nc)+iwcon(nc)
78!
79! Notetheless, be aware of:
80!
81! radliq(fisrt) .NE. ocond(nc)
82! i.e.:
83! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
84! (which is not trivial)
85!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
86!
87USE print_control_mod, ONLY: prt_level, lunout
88USE cloudth_mod, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
89USE lscp_tools_mod, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat, fallice_velocity
90USE ice_sursat_mod, ONLY : ice_sursat
91
92USE lscp_ini_mod, ONLY : seuil_neb, ninter, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
93USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radliq_snow, a_tr_sca
94
95
96IMPLICIT NONE
97
98!===============================================================================
99! VARIABLES DECLARATION
100!===============================================================================
101
102include "YOMCST.h"
103include "YOETHF.h"
104include "fisrtilp.h"
105include "nuage.h"
106
107! INPUT VARIABLES:
108!-----------------
109
110  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
111  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
112  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
113
114  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
115  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
116  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: t               ! temperature (K)
117  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: q               ! specific humidity [kg/kg]
118  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
119  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
120                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
121  LOGICAL,                         INTENT(IN)   :: ok_ice_sursat   ! flag to determine if ice sursaturation is activated
122
123  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
124
125
126  !Inputs associated with thermal plumes
127
128  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztv             ! virtual potential temperature [K]
129  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zqta            ! specific humidity within thermals [kg/kg]
130  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca           ! fraction of thermals within the mesh [-]
131  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: zpspsk          ! exner potential (p/100000)**(R/cp)
132  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ztla            ! liquid temperature within thermals [K]
133  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: zthl         ! liquid potential temperature [K]
134
135  ! INPUT/OUTPUT variables
136  !------------------------
137
138  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: ratqs            ! function of pressure that sets the large-scale
139  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: beta             ! conversion rate of condensed water
140
141
142  ! Input sursaturation en glace
143  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rneb_seri        ! fraction nuageuse en memoire
144 
145  ! OUTPUT variables
146  !-----------------
147
148  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
149  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
150  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
151  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
152  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
153  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
154  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radliq           ! condensed water used in the radiation scheme [kg/kg]
155  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
156  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
157  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! large-scale rainfall [kg/s/m2]
158  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! large-scale snowfall [kg/s/m2]
159  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl            ! saturation specific humidity wrt liquid [kg/kg]
160  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsats            ! saturation specific humidity wrt ice [kg/kg] 
161  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
162  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
163
164  ! fraction of aerosol scavenging through impaction and nucleation (for on-line)
165 
166  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa        ! scavenging fraction due tu impaction [-]
167  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl        ! scavenging fraction due tu nucleation [-]
168 
169  ! for supersaturation and contrails parameterisation
170 
171  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qclr             ! specific total water content in clear sky region [kg/kg]
172  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld             ! specific total water content in cloudy region [kg/kg]
173  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qss              ! specific total water content in supersat region [kg/kg]
174  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qvc              ! specific vapor content in clouds [kg/kg]
175  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebclr          ! mesh fraction of clear sky [-]   
176  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rnebss           ! mesh fraction of ISSR [-] 
177  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_ss         ! coefficient governing the ice nucleation RHi threshold [-]     
178  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr           ! ?
179  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr           ! ?
180  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2          ! ?         
181  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN          ! ?
182  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP          ! ?
183 
184
185  ! LOCAL VARIABLES:
186  !----------------
187
188  REAL qsl(klon), qsi(klon)
189  REAL zct, zcl
190  REAL ctot(klon,klev)
191  REAL ctot_vol(klon,klev)
192  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
193  REAL zdqsdT_raw(klon)
194  REAL gammasat(klon),dgammasatdt(klon)                ! coefficient to make cold condensation at the correct RH and derivative wrt T
195  REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
196  REAL cste
197  REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
198  REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
199  REAL erf
200  REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon)
201  REAL zrfl(klon), zrfln(klon), zqev, zqevt
202  REAL zifl(klon), zifln(klon), ziflprev(klon),zqev0,zqevi, zqevti
203  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon)
204  REAL zoliql(klon), zoliqi(klon)
205  REAL zt(klon)
206  REAL zdz(klon),zrho(klon,klev),iwc(klon)
207  REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon)
208  REAL zmelt,zrain,zsnow,zprecip
209  REAL dzfice(klon)
210  REAL zsolid
211  REAL qtot(klon), qzero(klon)
212  REAL dqsl(klon),dqsi(klon)
213  REAL smallestreal
214  !  Variables for Bergeron process
215  REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl
216  REAL zqpreci(klon), zqprecl(klon)
217  ! Variables precipitation energy conservation
218  REAL zmqc(klon)
219  REAL zalpha_tr
220  REAL zfrac_lessi
221  REAL zprec_cond(klon)
222  REAL zmair(klon), zcpair, zcpeau
223  REAL zlh_solid(klon), zm_solid         ! for liquid -> solid conversion
224  REAL zrflclr(klon), zrflcld(klon)
225  REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon)
226  REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon)
227  REAL ziflclr(klon), ziflcld(klon)
228  REAL znebprecipclr(klon), znebprecipcld(klon)
229  REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon)
230  REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon)
231  REAL velo(klon,klev), vr
232  REAL qlmpc, qimpc, rnebmpc
233  REAL radliqi(klon,klev), radliql(klon,klev)
234
235  INTEGER i, k, n, kk
236  INTEGER n_i(klon), iter
237  INTEGER ncoreczq
238  INTEGER mpc_bl_points(klon,klev)
239
240  LOGICAL lognormale(klon)
241  LOGICAL convergence(klon)
242  LOGICAL,SAVE :: appel1er
243  !$OMP THREADPRIVATE(appel1er)
244
245  CHARACTER (len = 20) :: modname = 'lscp'
246  CHARACTER (len = 80) :: abort_message
247 
248  DATA appel1er /.TRUE./
249
250!===============================================================================
251! INITIALISATION
252!===============================================================================
253
254! Few initial checks
255
256IF (iflag_t_glace.EQ.0) THEN
257    abort_message = 'lscp cannot be used if iflag_t_glace=0'
258    CALL abort_physic(modname,abort_message,1)
259ENDIF
260
261IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN
262    abort_message = 'lscp cannot be used without ice thermodynamics'
263    CALL abort_physic(modname,abort_message,1)
264ENDIF
265
266IF (iflag_fisrtilp_qsat .LT. 0) THEN
267    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
268    CALL abort_physic(modname,abort_message,1)
269ENDIF
270
271! Few initialisations
272
273znebprecip(:)=0.0
274ctot_vol(1:klon,1:klev)=0.0
275rneblsvol(1:klon,1:klev)=0.0
276smallestreal=1.e-9
277znebprecipclr(:)=0.0
278znebprecipcld(:)=0.0
279mpc_bl_points(:,:)=0
280
281IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
282
283! AA for 'safety' reasons
284zalpha_tr   = 0.
285zfrac_lessi = 0.
286beta(:,:)= 0.
287
288! Initialisation of variables:
289
290prfl(:,:) = 0.0
291psfl(:,:) = 0.0
292d_t(:,:) = 0.0
293d_q(:,:) = 0.0
294d_ql(:,:) = 0.0
295d_qi(:,:) = 0.0
296rneb(:,:) = 0.0
297radliq(:,:) = 0.0
298radicefrac(:,:) = 0.0
299frac_nucl(:,:) = 1.0
300frac_impa(:,:) = 1.0
301rain(:) = 0.0
302snow(:) = 0.0
303zoliq(:)=0.0
304zfice(:)=0.0
305dzfice(:)=0.0
306zqprecl(:)=0.0
307zqpreci(:)=0.0
308zrfl(:) = 0.0
309zifl(:) = 0.0
310ziflprev(:)=0.0
311zneb(:) = seuil_neb
312zrflclr(:) = 0.0
313ziflclr(:) = 0.0
314zrflcld(:) = 0.0
315ziflcld(:) = 0.0
316tot_zneb(:) = 0.0
317tot_znebn(:) = 0.0
318d_tot_zneb(:) = 0.0
319qzero(:) = 0.0
320
321!--ice supersaturation
322gamma_ss(:,:) = 1.0
323qss(:,:) = 0.0
324rnebss(:,:) = 0.0
325Tcontr(:,:) = missing_val
326qcontr(:,:) = missing_val
327qcontr2(:,:) = missing_val
328fcontrN(:,:) = 0.0
329fcontrP(:,:) = 0.0
330
331!c_iso: variable initialisation for iso
332
333
334!===============================================================================
335!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
336!===============================================================================
337
338ncoreczq=0
339
340DO k = klev, 1, -1
341
342    ! Initialisation temperature and specific humidity
343    DO i = 1, klon
344        zt(i)=t(i,k)
345        zq(i)=q(i,k)
346        !c_iso init of iso
347    ENDDO
348   
349    ! --------------------------------------------------------------------
350    ! P0> Thermalization of precipitation falling from the overlying layer
351    ! --------------------------------------------------------------------
352    ! Computes air temperature variation due to latent heat transported by
353    ! precipitation. Precipitation is then thermalized with the air in the
354    ! layer.
355    ! The precipitation should remain thermalized throughout the different
356    ! thermodynamical transformations. The corresponding water mass should
357    ! be added when calculating the layer's enthalpy change with
358    ! temperature
359    ! ---------------------------------------------------------------------
360   
361    IF (k.LE.klev-1) THEN
362
363        DO i = 1, klon
364 
365            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
366            ! no condensed water so cp=cp(vapor+dry air)
367            ! RVTMP2=rcpv/rcpd-1
368            zcpair=RCPD*(1.0+RVTMP2*zq(i))
369            zcpeau=RCPD*RVTMP2
370
371            ! zmqc: precipitation mass that has to be thermalized with
372            ! layer's air so that precipitation at the ground has the
373            ! same temperature as the lowermost layer
374            zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i)
375            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
376            zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
377             / (zcpair + zmqc(i)*zcpeau)
378
379        ENDDO
380 
381    ELSE
382 
383        DO i = 1, klon
384            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
385            zmqc(i) = 0.
386        ENDDO
387 
388    ENDIF
389
390    ! --------------------------------------------------------------------
391    ! P1> Precipitation evaporation/sublimation
392    ! --------------------------------------------------------------------
393    ! A part of the precipitation coming from above is evaporated/sublimated.
394    ! --------------------------------------------------------------------
395
396    IF (iflag_evap_prec.GE.1) THEN ! useless -> to be removed
397
398        ! Calculation of saturation specific humidity
399        ! depending on temperature:
400        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
401        ! wrt liquid water
402        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
403        ! wrt ice
404        CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
405
406        DO i = 1, klon
407
408            ! if precipitation
409            IF (zrfl(i)+zifl(i).GT.0.) THEN
410
411            ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4).
412            ! c_iso: likely important to distinguish cs from neb isotope precipitation
413
414                IF (iflag_evap_prec.EQ.4) THEN
415                    zrfl(i) = zrflclr(i)
416                    zifl(i) = ziflclr(i)
417                ENDIF
418
419                IF (iflag_evap_prec.EQ.1) THEN
420                    znebprecip(i)=zneb(i)
421                ELSE
422                    znebprecip(i)=MAX(zneb(i),znebprecip(i))
423                ENDIF
424
425                IF (iflag_evap_prec.EQ.4) THEN
426                ! Max evaporation not to saturate the whole mesh
427                    zqev0 = MAX(0.0, zqs(i)-zq(i))
428                ELSE
429                ! Max evap not to saturate the fraction below the cloud
430                    zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
431                ENDIF
432
433                ! Evaporation of liquid precipitation coming from above
434                ! dP/dz=beta*(1-q/qsat)*sqrt(P)
435                ! formula from Sundquist 1988, Klemp & Wilhemson 1978
436                ! LTP: evaporation only in the clear sky part (iflag_evap_prec=4)
437
438                IF (iflag_evap_prec.EQ.3) THEN
439                    zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
440                    *SQRT(zrfl(i)/max(1.e-4,znebprecip(i)))        &
441                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
442                ELSE IF (iflag_evap_prec.EQ.4) THEN
443                     zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) &
444                    *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
445                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
446                ELSE
447                    zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
448                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
449                ENDIF
450
451                zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
452                    *RG*dtime/(paprs(i,k)-paprs(i,k+1))
453
454                ! sublimation of the solid precipitation coming from above
455                IF (iflag_evap_prec.EQ.3) THEN
456                    zqevti = znebprecip(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) &
457                    *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
458                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
459                ELSE IF (iflag_evap_prec.EQ.4) THEN
460                     zqevti = znebprecipclr(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) &
461                    *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
462                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
463                ELSE
464                    zqevti = 1.*coef_eva_i*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
465                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
466                ENDIF
467
468                zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
469                *RG*dtime/(paprs(i,k)-paprs(i,k+1))
470
471                ! A. JAM
472                ! Evaporation limit: we ensure that the layer's fraction below
473                ! the cloud or the whole mesh (depending on iflag_evap_prec)
474                ! does not reach saturation. In this case, we
475                ! redistribute zqev0 conserving the ratio liquid/ice
476
477                IF (zqevt+zqevti.GT.zqev0) THEN
478                    zqev=zqev0*zqevt/(zqevt+zqevti)
479                    zqevi=zqev0*zqevti/(zqevt+zqevti)
480                ELSE
481                    zqev=zqevt
482                    zqevi=zqevti
483                ENDIF
484
485
486                ! New solid and liquid precipitation fluxes after evap and sublimation
487                zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))  &
488                         /RG/dtime)
489                zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
490                         /RG/dtime)
491
492
493                ! vapor, temperature, precip fluxes update
494                zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
495                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
496                ! precip thermalization
497                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
498                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
499                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
500                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
501                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
502                + (zifln(i)-zifl(i))                      &
503                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
504                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
505
506                ! New values of liquid and solid precipitation
507                zrfl(i) = zrfln(i)
508                zifl(i) = zifln(i)
509
510                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
511                !                         due to evap + sublim                                     
512                                           
513
514                IF (iflag_evap_prec.EQ.4) THEN
515                    zrflclr(i) = zrfl(i)
516                    ziflclr(i) = zifl(i)
517                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN
518                        znebprecipclr(i) = 0.0
519                    ENDIF
520                    zrfl(i) = zrflclr(i) + zrflcld(i)
521                    zifl(i) = ziflclr(i) + ziflcld(i)
522                ENDIF
523
524                ! c_iso duplicate for isotopes or loop on isotopes
525
526                ! Melting:
527                zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! JYG
528                ! precip fraction that is melted
529                zmelt = MIN(MAX(zmelt,0.),1.)
530
531                ! Ice melting
532                IF (iflag_evap_prec.EQ.4) THEN
533                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
534                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
535                    zrfl(i)=zrflclr(i)+zrflcld(i)
536                ELSE
537                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
538                ENDIF
539                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
540
541                ! Latent heat of melting with precipitation thermalization
542                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
543                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
544               
545                IF (iflag_evap_prec.EQ.4) THEN
546                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
547                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
548                    zifl(i)=ziflclr(i)+ziflcld(i)
549                ELSE
550                    zifl(i)=zifl(i)*(1.-zmelt)
551                ENDIF
552                !c_iso: same for isotopic precip.
553
554            ELSE
555                ! if no precip, we reinitialize the cloud fraction used for the precip to 0
556                znebprecip(i)=0.
557
558            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
559
560        ENDDO
561
562    ENDIF ! (iflag_evap_prec>=1)
563
564    ! --------------------------------------------------------------------
565    ! End precip evaporation
566    ! --------------------------------------------------------------------
567   
568    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
569    !-------------------------------------------------------
570
571     qtot(:)=zq(:)+zmqc(:)
572     CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
573     DO i = 1, klon
574            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
575            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
576            IF (zq(i) .LT. 1.e-15) THEN
577                ncoreczq=ncoreczq+1
578                zq(i)=1.e-15
579            ENDIF
580            ! c_iso: do something similar for isotopes
581
582     ENDDO
583   
584    ! --------------------------------------------------------------------
585    ! P2> Cloud formation
586    !---------------------------------------------------------------------
587    !
588    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
589    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
590    ! is prescribed and depends on large scale variables and boundary layer
591    ! properties)
592    ! The decrease in condensed part due tu latent heating is taken into
593    ! account
594    ! -------------------------------------------------------------------
595
596        ! P2.1> With the PDFs (log-normal, bigaussian)
597        ! cloud propertues calculation with the initial values of t and q
598        ! ----------------------------------------------------------------
599
600        ! initialise gammasat and qincloud_mpc
601        gammasat(:)=1.
602        qincloud_mpc(:)=0.
603
604        IF (iflag_cld_th.GE.5) THEN
605
606             IF (iflag_mpc_bl .LT. 1) THEN
607
608                IF (iflag_cloudth_vert.LE.2) THEN
609
610                    CALL cloudth(klon,klev,k,ztv,                  &
611                         zq,zqta,fraca,                            &
612                         qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
613                         ratqs,zqs,t)
614
615                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
616
617                    CALL cloudth_v3(klon,klev,k,ztv,                        &
618                         zq,zqta,fraca,                                     &
619                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
620                         ratqs,zqs,t)
621
622                !Jean Jouhaud's version, Decembre 2018
623                !-------------------------------------
624
625                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
626
627                    CALL cloudth_v6(klon,klev,k,ztv,                        &
628                         zq,zqta,fraca,                                     &
629                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
630                         ratqs,zqs,t)
631
632                ENDIF
633
634         ELSE
635                ! cloudth_v3 + cold microphysical considerations
636                ! + stationary mixed-phase cloud model
637
638                    CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, &
639                         t,ztv,zq,zqta,fraca, zpspsk,                      &
640                         paprs,pplay,ztla,zthl,ratqs,zqs,psfl,             &
641                         qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol)
642
643         ENDIF ! iflag_mpc_bl
644
645                DO i=1,klon
646                    rneb(i,k)=ctot(i,k)
647                    rneblsvol(i,k)=ctot_vol(i,k)
648                    zqn(i)=qcloud(i)
649                ENDDO
650
651        ENDIF
652
653        IF (iflag_cld_th .LE. 4) THEN
654           
655                ! lognormal
656            lognormale = .TRUE.
657
658        ELSEIF (iflag_cld_th .GE. 6) THEN
659           
660                ! lognormal distribution when no thermals
661            lognormale = fraca(:,k) < 1e-10
662
663        ELSE
664                ! When iflag_cld_th=5, we always assume
665                ! bi-gaussian distribution
666            lognormale = .FALSE.
667       
668        ENDIF
669
670        DT(:) = 0.
671        n_i(:)=0
672        Tbef(:)=zt(:)
673        qlbef(:)=0.
674
675        ! Treatment of non-boundary layer clouds (lognormale)
676        ! condensation with qsat(T) variation (adaptation)
677        ! Iterative Loop to converge towards qsat
678
679        DO iter=1,iflag_fisrtilp_qsat+1
680
681                DO i=1,klon
682
683                    ! convergence = .true. until when convergence is not satisfied
684                    convergence(i)=ABS(DT(i)).GT.DDT0
685
686                    IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
687
688                        ! if not convergence:
689
690                        ! P2.2.1> cloud fraction and condensed water mass calculation
691                        ! Calculated variables:
692                        ! rneb : cloud fraction
693                        ! zqn : total water within the cloud
694                        ! zcond: mean condensed water within the mesh
695                        ! rhcl: clear-sky relative humidity
696                        !---------------------------------------------------------------
697
698                        ! new temperature:
699                        Tbef(i)=Tbef(i)+DT(i)
700
701                        ! Rneb, qzn and zcond for lognormal PDFs
702                        qtot(i)=zq(i)+zmqc(i)
703
704                      ENDIF
705
706                  ENDDO
707
708                  ! Calculation of saturation specific humidity and ice fraction
709                  CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
710                  CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
711                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
712                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
713                  CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
714
715                  DO i=1,klon
716
717                      IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
718
719                        zpdf_sig(i)=ratqs(i,k)*zq(i)
720                        zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
721                        zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
722                        zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
723                        zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
724                        zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
725                        zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
726                        zpdf_e1(i)=1.-erf(zpdf_e1(i))
727                        zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
728                        zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
729                        zpdf_e2(i)=1.-erf(zpdf_e2(i))
730
731                        !--ice sursaturation by Audran
732                        IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN
733
734                          IF (zpdf_e1(i).LT.1.e-10) THEN
735                              rneb(i,k)=0.
736                              zqn(i)=gammasat(i)*zqs(i)
737                          ELSE
738                              rneb(i,k)=0.5*zpdf_e1(i)
739                              zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
740                          ENDIF
741
742                          rnebss(i,k)=0.0   !--added by OB (needed because of the convergence loop on time)
743                          fcontrN(i,k)=0.0  !--idem
744                          fcontrP(i,k)=0.0  !--idem
745                          qss(i,k)=0.0      !--idem
746
747                        ELSE
748
749                        !------------------------------------
750                        ! ICE SUPERSATURATION
751                        !------------------------------------
752
753                        CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), &
754                             gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k),                 &
755                             rneb(i,k), zqn(i), rnebss(i,k), qss(i,k),                                   &
756                             Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
757
758                        ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
759
760                        ! If vertical heterogeneity, change fraction by volume as well
761                        IF (iflag_cloudth_vert.GE.3) THEN
762                            ctot_vol(i,k)=rneb(i,k)
763                            rneblsvol(i,k)=ctot_vol(i,k)
764                        ENDIF
765
766
767                       ! P2.2.2> Approximative calculation of temperature variation DT
768                       !           due to condensation.
769                       ! Calculated variables:
770                       ! dT : temperature change due to condensation
771                       !---------------------------------------------------------------
772
773               
774                        IF (zfice(i).LT.1) THEN
775                            cste=RLVTT
776                        ELSE
777                            cste=RLSTT
778                        ENDIF
779
780                        qlbef(i)=max(0.,zqn(i)-zqs(i))
781                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
782                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
783                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
784                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
785                              *qlbef(i)*dzfice(i)
786                        DT(i)=num/denom
787                        n_i(i)=n_i(i)+1
788
789                    ENDIF ! end convergence
790
791                ENDDO     ! end loop on i
792
793        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
794
795        ! P2.3> Final quantities calculation
796        ! Calculated variables:
797        !   rneb : cloud fraction
798        !   zcond: mean condensed water in the mesh
799        !   zqn  : mean water vapor in the mesh
800        !   zt   : temperature
801        !   rhcl : clear-sky relative humidity
802        ! ----------------------------------------------------------------
803
804        ! Water vapor update, Phase determination and subsequent latent heat exchange
805
806        ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
807        CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
808
809        DO i=1, klon
810
811
812            IF (mpc_bl_points(i,k) .GT. 0) THEN
813                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
814                ! following line is very strange and probably wrong
815                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
816                ! water vapor update and partition function if thermals
817                zq(i) = zq(i) - zcond(i)       
818                zfice(i)=icefrac_mpc(i,k)
819
820            ELSE
821
822                ! Checks on rneb, rhcl and zqn
823                IF (rneb(i,k) .LE. 0.0) THEN
824                    zqn(i) = 0.0
825                    rneb(i,k) = 0.0
826                    zcond(i) = 0.0
827                    rhcl(i,k)=zq(i)/zqs(i)
828                ELSE IF (rneb(i,k) .GE. 1.0) THEN
829                    zqn(i) = zq(i)
830                    rneb(i,k) = 1.0
831                    zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))
832                    rhcl(i,k)=1.0
833                ELSE
834                    zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))*rneb(i,k)
835                    ! following line is very strange and probably wrong:
836                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
837                ENDIF
838
839                ! water vapor update
840                zq(i) = zq(i) - zcond(i)
841
842            ENDIF
843
844            ! c_iso : routine that computes in-cloud supersaturation
845            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
846
847            ! temperature update due to phase change
848            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
849            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
850                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
851        ENDDO
852
853        ! If vertical heterogeneity, change volume fraction
854        IF (iflag_cloudth_vert .GE. 3) THEN
855          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
856          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
857        ENDIF
858
859    ! ----------------------------------------------------------------
860    ! P3> Precipitation formation
861    ! ----------------------------------------------------------------
862
863    ! LTP:
864    IF (iflag_evap_prec==4) THEN
865
866        !Partitionning between precipitation coming from clouds and that coming from CS
867
868        !0) Calculate tot_zneb, total cloud fraction above the cloud
869        !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
870       
871        DO i=1, klon
872                tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
873                        /(1-min(zneb(i),1-smallestreal))
874                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
875                tot_zneb(i) = tot_znebn(i)
876
877
878                !1) Cloudy to clear air
879                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
880                IF (znebprecipcld(i) .GT. 0) THEN
881                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
882                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
883                ELSE
884                        d_zrfl_cld_clr(i) = 0.
885                        d_zifl_cld_clr(i) = 0.
886                ENDIF
887
888                !2) Clear to cloudy air
889                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
890                IF (znebprecipclr(i) .GT. 0) THEN
891                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
892                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
893                ELSE
894                        d_zrfl_clr_cld(i) = 0.
895                        d_zifl_clr_cld(i) = 0.
896                ENDIF
897
898                !Update variables
899                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
900                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
901                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
902                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
903                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
904                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
905
906                ! c_iso  do the same thing for isotopes precip
907        ENDDO
908    ENDIF
909
910    ! Initialisation of zoliq and radliq variables
911
912    DO i = 1, klon
913            zoliq(i) = zcond(i)
914            zoliqi(i)= zoliq(i)*zfice(i)
915            zoliql(i)= zoliq(i)*(1.-zfice(i))
916            ! c_iso : initialisation of zoliq* also for isotopes
917            zrho(i,k)  = pplay(i,k) / zt(i) / RD
918            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
919            iwc(i)   = 0.
920            zneb(i)  = MAX(rneb(i,k),seuil_neb)
921            radliq(i,k) = zoliq(i)/REAL(ninter+1)
922            radicefrac(i,k) = zfice(i)
923            radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
924            radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
925    ENDDO
926
927    ! Autoconversion
928    ! -------------------------------------------------------------------------------
929     
930    DO n = 1, ninter
931
932        DO i=1,klon
933            IF (rneb(i,k).GT.0.0) THEN
934                iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
935            ENDIF
936        ENDDO
937
938        CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
939
940        DO i = 1, klon
941
942            IF (rneb(i,k).GT.0.0) THEN
943
944                ! Initialization of zrain, zsnow and zprecip:
945                zrain=0.
946                zsnow=0.
947                zprecip=0.
948                ! c_iso same init for isotopes. Externalisation?
949
950                IF (zneb(i).EQ.seuil_neb) THEN
951                    zprecip = 0.0
952                    zsnow = 0.0
953                    zrain= 0.0
954                ELSE
955                    ! water quantity to remove: zchau (Sundqvist, 1978)
956                    ! same thing for the ice: zfroi (Zender & Kiehl, 1997)
957                    IF (ptconv(i,k)) THEN ! if convective point
958                        zcl=cld_lc_con
959                        zct=1./cld_tau_con
960                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*zfice(i)
961                    ELSE
962                        zcl=cld_lc_lsc
963                        zct=1./cld_tau_lsc
964                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
965                            *velo(i,k) * zfice(i)
966                    ENDIF
967
968                    ! if vertical heterogeneity is taken into account, we use
969                    ! the "true" volume fraction instead of a modified
970                    ! surface fraction (which is larger and artificially
971                    ! reduces the in-cloud water).
972
973
974                    ! warning: here should be zoliql?
975                    IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
976                        zchau = zct   *dtime/REAL(ninter) * zoliq(i) &
977                        *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl)**2)) *(1.-zfice(i))
978                    ELSE
979                        zchau = zct   *dtime/REAL(ninter) * zoliq(i) &
980                        *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl)**2)) *(1.-zfice(i)) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
981                    ENDIF
982
983                    zrain   = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
984                    zsnow   = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
985                    zprecip = MAX(zrain + zsnow,0.0)
986
987                ENDIF
988
989                zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
990                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
991                zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
992
993                ! c_iso: call isotope_conversion (for readibility) or duplicate
994
995                radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
996                radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1)
997                radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1)
998
999            ENDIF ! rneb >0
1000
1001        ENDDO  ! i = 1,klon
1002
1003    ENDDO      ! n = 1,niter
1004
1005    ! Precipitation flux calculation
1006
1007    DO i = 1, klon
1008           
1009            IF (iflag_evap_prec.EQ.4) THEN
1010                ziflprev(i)=ziflcld(i)
1011            ELSE
1012                ziflprev(i)=zifl(i)*zneb(i)
1013            ENDIF
1014
1015            IF (rneb(i,k) .GT. 0.0) THEN
1016
1017               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
1018               ! If T<0C, liquid precip are converted into ice, which leads to an increase in
1019               ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
1020               ! taken into account through a linearization of the equations and by approximating
1021               ! the liquid precipitation process with a "threshold" process. We assume that
1022               ! condensates are not modified during this operation. Liquid precipitation is
1023               ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
1024               ! Water vapor increases as well
1025               ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
1026
1027                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1028                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1029                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1030                    coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1031                    ! Computation of DT if all the liquid precip freezes
1032                    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1033                    ! T should not exceed the freezing point
1034                    ! that is Delta > RTT-zt(i)
1035                    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1036                    zt(i)  = zt(i) + DeltaT
1037                    ! water vaporization due to temp. increase
1038                    Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1039                    ! we add this vaporized water to the total vapor and we remove it from the precip
1040                    zq(i) = zq(i) +  Deltaq
1041                    ! The three "max" lines herebelow protect from rounding errors
1042                    zcond(i) = max( zcond(i)- Deltaq, 0. )
1043                    ! liquid precipitation converted to ice precip
1044                    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1045                    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1046                    ! iced water budget
1047                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1048
1049                    ! c_iso : mv here condensation of isotopes + redispatchage en precip
1050
1051                IF (iflag_evap_prec.EQ.4) THEN
1052                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
1053                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1054                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1055                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1056                    znebprecipcld(i) = rneb(i,k)
1057                    zrfl(i) = zrflcld(i) + zrflclr(i)
1058                    zifl(i) = ziflcld(i) + ziflclr(i)
1059                ELSE
1060                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
1061                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1062                    zifl(i) = zifl(i)+ zqpreci(i)        &
1063                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1064                ENDIF
1065                ! c_iso : same for isotopes
1066 
1067           ENDIF ! rneb>0
1068
1069   ENDDO
1070
1071    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
1072    ! if iflag_evap_pre=4
1073    IF (iflag_evap_prec.EQ.4) THEN
1074
1075        DO i=1,klon
1076
1077            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
1078                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
1079                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
1080            ELSE
1081                znebprecipclr(i)=0.0
1082            ENDIF
1083
1084            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
1085                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
1086                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
1087            ELSE
1088                znebprecipcld(i)=0.0
1089            ENDIF
1090
1091        ENDDO
1092
1093    ENDIF
1094
1095    ! End of precipitation formation
1096    ! --------------------------------
1097
1098    ! Outputs:
1099    ! Precipitation fluxes at layer interfaces
1100    ! and temperature and water species tendencies
1101    DO i = 1, klon
1102        psfl(i,k)=zifl(i)
1103        prfl(i,k)=zrfl(i)
1104        d_ql(i,k) = (1-zfice(i))*zoliq(i)
1105        d_qi(i,k) = zfice(i)*zoliq(i)
1106        d_q(i,k) = zq(i) - q(i,k)
1107        ! c_iso: same for isotopes
1108        d_t(i,k) = zt(i) - t(i,k)
1109    ENDDO
1110
1111    ! Calculation of the concentration of condensates seen by the radiation scheme
1112    ! for liquid phase, we take radliql
1113    ! for ice phase, we take radliqi if we neglect snowfall, otherwise (ok_radliq_snow=true)
1114    ! we recaulate radliqi to account for contributions from the precipitation flux
1115
1116    IF ((ok_radliq_snow) .AND. (k .LT. klev)) THEN
1117        ! for the solid phase (crystals + snowflakes)
1118        ! we recalculate radliqi by summing
1119        ! the ice content calculated in the mesh
1120        ! + the contribution of the non-evaporated snowfall
1121        ! from the overlying layer
1122        DO i=1,klon
1123           IF (ziflprev(i) .NE. 0.0) THEN
1124              radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
1125           ELSE
1126              radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
1127           ENDIF
1128           radliq(i,k)=radliql(i,k)+radliqi(i,k)
1129        ENDDO
1130    ENDIF
1131
1132    ! caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme
1133    DO i=1,klon
1134        IF (radliq(i,k) .GT. 0.) THEN
1135            radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.)
1136        ENDIF
1137    ENDDO
1138
1139    ! ----------------------------------------------------------------
1140    ! P4> Wet scavenging
1141    ! ----------------------------------------------------------------
1142
1143    !Scavenging through nucleation in the layer
1144
1145    DO i = 1,klon
1146       
1147        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1148            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1149        ELSE
1150            beta(i,k) = 0.
1151        ENDIF
1152
1153        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1154
1155        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1156
1157            IF (t(i,k) .GE. t_glace_min) THEN
1158                zalpha_tr = a_tr_sca(3)
1159            ELSE
1160                zalpha_tr = a_tr_sca(4)
1161            ENDIF
1162
1163            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1164            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1165
1166            ! Nucleation with a factor of -1 instead of -0.5
1167            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1168
1169        ENDIF
1170
1171    ENDDO
1172
1173    ! Scavenging through impaction in the underlying layer
1174
1175    DO kk = k-1, 1, -1
1176
1177        DO i = 1, klon
1178
1179            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1180
1181                IF (t(i,kk) .GE. t_glace_min) THEN
1182                    zalpha_tr = a_tr_sca(1)
1183                ELSE
1184                    zalpha_tr = a_tr_sca(2)
1185                ENDIF
1186
1187              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1188              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1189
1190            ENDIF
1191
1192        ENDDO
1193
1194    ENDDO
1195
1196    !--save some variables for ice supersaturation
1197    !
1198    DO i = 1, klon
1199        ! for memory
1200        rneb_seri(i,k) = rneb(i,k)
1201
1202        ! for diagnostics
1203        rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
1204
1205        qvc(i,k) = zqs(i) * rneb(i,k)
1206        qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k))  !--added by OB because of pathologic cases with the lognormal
1207        qcld(i,k) = qvc(i,k) + zcond(i)
1208   ENDDO
1209   !q_sat
1210   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:))
1211   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
1212
1213ENDDO
1214
1215!======================================================================
1216!                      END OF VERTICAL LOOP
1217!======================================================================
1218
1219  ! Rain or snow at the surface (depending on the first layer temperature)
1220  DO i = 1, klon
1221      snow(i) = zifl(i)
1222      rain(i) = zrfl(i)
1223      ! c_iso final output
1224  ENDDO
1225
1226  IF (ncoreczq>0) THEN
1227      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1228  ENDIF
1229
1230END SUBROUTINE lscp
1231!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1232
1233END MODULE lscp_mod
Note: See TracBrowser for help on using the repository browser.