source: LMDZ6/branches/blowing_snow/libf/phylmd/lscp_mod.F90 @ 4484

Last change on this file since 4484 was 4425, checked in by evignon, 21 months ago

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