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

Last change on this file since 4553 was 4535, checked in by evignon, 14 months ago

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