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

Last change on this file since 4412 was 4412, checked in by evignon, 17 months ago

travail de l'atelier nuages du 30/01/23: on renomme cldliq (ou radliq) en radocond
car le nom est tres trompeur + on ajoute des commentaires dans lscp_mod
File size: 49.8 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
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           ! ?
181  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr           ! ?
182  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2          ! ?         
183  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN          ! ?
184  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP          ! ?
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 convergence(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
504                ! precip thermalization
505                zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
506                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
507                zt(i) = zt(i) + (zrfln(i)-zrfl(i))        &
508                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
509                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
510                + (zifln(i)-zifl(i))                      &
511                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
512                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
513
514                ! New values of liquid and solid precipitation
515                zrfl(i) = zrfln(i)
516                zifl(i) = zifln(i)
517
518                ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
519                !                         due to evap + sublim                                     
520                                           
521
522                IF (iflag_evap_prec.EQ.4) THEN
523                    zrflclr(i) = zrfl(i)
524                    ziflclr(i) = zifl(i)
525                    IF(zrflclr(i) + ziflclr(i).LE.0) THEN
526                        znebprecipclr(i) = 0.0
527                    ENDIF
528                    zrfl(i) = zrflclr(i) + zrflcld(i)
529                    zifl(i) = ziflclr(i) + ziflcld(i)
530                ENDIF
531
532                ! c_iso duplicate for isotopes or loop on isotopes
533
534                ! Melting:
535                zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
536                ! precip fraction that is melted
537                zmelt = MIN(MAX(zmelt,0.),1.)
538
539                ! update of rainfall and snowfall due to melting
540                IF (iflag_evap_prec.EQ.4) THEN
541                    zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
542                    zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
543                    zrfl(i)=zrflclr(i)+zrflcld(i)
544
545                    ziflclr(i)=ziflclr(i)*(1.-zmelt)
546                    ziflcld(i)=ziflcld(i)*(1.-zmelt)
547                    zifl(i)=ziflclr(i)+ziflcld(i)
548
549                ELSE
550                    zrfl(i)=zrfl(i)+zmelt*zifl(i)
551
552                    zifl(i)=zifl(i)*(1.-zmelt)
553                ENDIF
554
555
556                ! c_iso: melting of isotopic precipi with zmelt (no fractionation)
557
558                ! Latent heat of melting with precipitation thermalization
559                zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
560                *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
561               
562
563            ELSE
564                ! if no precip, we reinitialize the cloud fraction used for the precip to 0
565                znebprecip(i)=0.
566
567            ENDIF ! (zrfl(i)+zifl(i).GT.0.)
568
569        ENDDO
570
571    ENDIF ! (iflag_evap_prec>=1)
572
573    ! --------------------------------------------------------------------
574    ! End precip evaporation
575    ! --------------------------------------------------------------------
576   
577    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
578    !-------------------------------------------------------
579
580     qtot(:)=zq(:)+zmqc(:)
581     CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
582     DO i = 1, klon
583            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
584            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
585            IF (zq(i) .LT. 1.e-15) THEN
586                ncoreczq=ncoreczq+1
587                zq(i)=1.e-15
588            ENDIF
589            ! c_iso: do something similar for isotopes
590
591     ENDDO
592   
593    ! --------------------------------------------------------------------
594    ! P2> Cloud formation
595    !---------------------------------------------------------------------
596    !
597    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
598    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
599    ! is prescribed and depends on large scale variables and boundary layer
600    ! properties)
601    ! The decrease in condensed part due tu latent heating is taken into
602    ! account
603    ! -------------------------------------------------------------------
604
605        ! P2.1> With the PDFs (log-normal, bigaussian)
606        ! cloud propertues calculation with the initial values of t and q
607        ! ----------------------------------------------------------------
608
609        ! initialise gammasat and qincloud_mpc
610        gammasat(:)=1.
611        qincloud_mpc(:)=0.
612
613        IF (iflag_cld_th.GE.5) THEN
614
615             IF (iflag_mpc_bl .LT. 1) THEN
616
617                IF (iflag_cloudth_vert.LE.2) THEN
618
619                    CALL cloudth(klon,klev,k,ztv,                  &
620                         zq,zqta,fraca,                            &
621                         qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, &
622                         ratqs,zqs,t)
623
624                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
625
626                    CALL cloudth_v3(klon,klev,k,ztv,                        &
627                         zq,zqta,fraca,                                     &
628                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
629                         ratqs,zqs,t)
630
631                !Jean Jouhaud's version, Decembre 2018
632                !-------------------------------------
633
634                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
635
636                    CALL cloudth_v6(klon,klev,k,ztv,                        &
637                         zq,zqta,fraca,                                     &
638                         qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, &
639                         ratqs,zqs,t)
640
641                ENDIF
642
643         ELSE
644                ! cloudth_v3 + cold microphysical considerations
645                ! + stationary mixed-phase cloud model
646
647                    CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, &
648                         t,ztv,zq,zqta,fraca, zpspsk,                      &
649                         paprs,pplay,ztla,zthl,ratqs,zqs,psfl,             &
650                         qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol)
651
652         ENDIF ! iflag_mpc_bl
653
654                DO i=1,klon
655                    rneb(i,k)=ctot(i,k)
656                    rneblsvol(i,k)=ctot_vol(i,k)
657                    zqn(i)=qcloud(i)
658                ENDDO
659
660        ENDIF
661
662        IF (iflag_cld_th .LE. 4) THEN
663           
664                ! lognormal
665            lognormale = .TRUE.
666
667        ELSEIF (iflag_cld_th .GE. 6) THEN
668           
669                ! lognormal distribution when no thermals
670            lognormale = fraca(:,k) < 1e-10
671
672        ELSE
673                ! When iflag_cld_th=5, we always assume
674                ! bi-gaussian distribution
675            lognormale = .FALSE.
676       
677        ENDIF
678
679        DT(:) = 0.
680        n_i(:)=0
681        Tbef(:)=zt(:)
682        qlbef(:)=0.
683
684        ! Treatment of non-boundary layer clouds (lognormale)
685        ! condensation with qsat(T) variation (adaptation)
686        ! Iterative Loop to converge towards qsat
687
688        DO iter=1,iflag_fisrtilp_qsat+1
689
690                DO i=1,klon
691
692                    ! convergence = .true. until when convergence is not satisfied
693                    convergence(i)=ABS(DT(i)).GT.DDT0
694
695                    IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
696
697                        ! if not convergence:
698
699                        ! P2.2.1> cloud fraction and condensed water mass calculation
700                        ! Calculated variables:
701                        ! rneb : cloud fraction
702                        ! zqn : total water within the cloud
703                        ! zcond: mean condensed water within the mesh
704                        ! rhcl: clear-sky relative humidity
705                        !---------------------------------------------------------------
706
707                        ! new temperature:
708                        Tbef(i)=Tbef(i)+DT(i)
709
710                        ! Rneb, qzn and zcond for lognormal PDFs
711                        qtot(i)=zq(i)+zmqc(i)
712
713                      ENDIF
714
715                  ENDDO
716
717                  ! Calculation of saturation specific humidity and ice fraction
718                  CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
719                  CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
720                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
721                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
722                  CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:))
723
724                  DO i=1,klon
725
726                      IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
727
728                        zpdf_sig(i)=ratqs(i,k)*zq(i)
729                        zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
730                        zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
731                        zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
732                        zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
733                        zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
734                        zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
735                        zpdf_e1(i)=1.-erf(zpdf_e1(i))
736                        zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
737                        zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
738                        zpdf_e2(i)=1.-erf(zpdf_e2(i))
739
740                        !--ice sursaturation by Audran
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
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                        DT(i)=num/denom
796                        n_i(i)=n_i(i)+1
797
798                    ENDIF ! end convergence
799
800                ENDDO     ! end loop on i
801
802        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
803
804        ! P2.3> Final quantities calculation
805        ! Calculated variables:
806        !   rneb : cloud fraction
807        !   zcond: mean condensed water in the mesh
808        !   zqn  : mean water vapor in the mesh
809        !   zt   : temperature
810        !   rhcl : clear-sky relative humidity
811        ! ----------------------------------------------------------------
812
813        ! Water vapor update, Phase determination and subsequent latent heat exchange
814
815        ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs)
816        CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))
817
818        DO i=1, klon
819
820
821            IF (mpc_bl_points(i,k) .GT. 0) THEN
822                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
823                ! following line is very strange and probably wrong
824                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
825                ! water vapor update and partition function if thermals
826                zq(i) = zq(i) - zcond(i)       
827                zfice(i)=icefrac_mpc(i,k)
828
829            ELSE
830
831                ! Checks on rneb, rhcl and zqn
832                IF (rneb(i,k) .LE. 0.0) THEN
833                    zqn(i) = 0.0
834                    rneb(i,k) = 0.0
835                    zcond(i) = 0.0
836                    rhcl(i,k)=zq(i)/zqs(i)
837                ELSE IF (rneb(i,k) .GE. 1.0) THEN
838                    zqn(i) = zq(i)
839                    rneb(i,k) = 1.0
840                    zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))
841                    rhcl(i,k)=1.0
842                ELSE
843                    zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))*rneb(i,k)
844                    ! following line is very strange and probably wrong:
845                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
846                ENDIF
847
848                ! water vapor update
849                zq(i) = zq(i) - zcond(i)
850
851            ENDIF
852
853            ! c_iso : routine that computes in-cloud supersaturation
854            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
855
856            ! temperature update due to phase change
857            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
858            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
859                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
860        ENDDO
861
862        ! If vertical heterogeneity, change volume fraction
863        IF (iflag_cloudth_vert .GE. 3) THEN
864          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
865          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
866        ENDIF
867
868    ! ----------------------------------------------------------------
869    ! P3> Precipitation formation
870    ! ----------------------------------------------------------------
871
872    ! LTP:
873    IF (iflag_evap_prec==4) THEN
874
875        !Partitionning between precipitation coming from clouds and that coming from CS
876
877        !0) Calculate tot_zneb, total cloud fraction above the cloud
878        !assuming a maximum-random overlap (voir Jakob and Klein, 2000)
879       
880        DO i=1, klon
881                tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
882                        /(1-min(zneb(i),1-smallestreal))
883                d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
884                tot_zneb(i) = tot_znebn(i)
885
886
887                !1) Cloudy to clear air
888                d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
889                IF (znebprecipcld(i) .GT. 0) THEN
890                        d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
891                        d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
892                ELSE
893                        d_zrfl_cld_clr(i) = 0.
894                        d_zifl_cld_clr(i) = 0.
895                ENDIF
896
897                !2) Clear to cloudy air
898                d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
899                IF (znebprecipclr(i) .GT. 0) THEN
900                        d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
901                        d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
902                ELSE
903                        d_zrfl_clr_cld(i) = 0.
904                        d_zifl_clr_cld(i) = 0.
905                ENDIF
906
907                !Update variables
908                znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
909                znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
910                zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
911                ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
912                zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
913                ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
914
915                ! c_iso  do the same thing for isotopes precip
916        ENDDO
917    ENDIF
918
919    ! Initialisation of zoliq and radocond variables
920
921    DO i = 1, klon
922            zoliq(i) = zcond(i)
923            zoliqi(i)= zoliq(i)*zfice(i)
924            zoliql(i)= zoliq(i)*(1.-zfice(i))
925            ! c_iso : initialisation of zoliq* also for isotopes
926            zrho(i,k)  = pplay(i,k) / zt(i) / RD
927            zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
928            iwc(i)   = 0.
929            zneb(i)  = MAX(rneb(i,k),seuil_neb)
930            radocond(i,k) = zoliq(i)/REAL(ninter+1)
931            radicefrac(i,k) = zfice(i)
932            radocondi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1)
933            radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1)
934    ENDDO
935
936    ! Autoconversion
937    ! -------------------------------------------------------------------------------
938     
939    DO n = 1, ninter
940
941        DO i=1,klon
942            IF (rneb(i,k).GT.0.0) THEN
943                iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
944            ENDIF
945        ENDDO
946
947        CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
948
949        DO i = 1, klon
950
951            IF (rneb(i,k).GT.0.0) THEN
952
953                ! Initialization of zrain, zsnow and zprecip:
954                zrain=0.
955                zsnow=0.
956                zprecip=0.
957                ! c_iso same init for isotopes. Externalisation?
958
959                IF (zneb(i).EQ.seuil_neb) THEN
960                    zprecip = 0.0
961                    zsnow = 0.0
962                    zrain= 0.0
963                ELSE
964                    ! water quantity to remove: zchau (Sundqvist, 1978)
965                    ! same thing for the ice: zfroi (Zender & Kiehl, 1997)
966                    IF (ptconv(i,k)) THEN ! if convective point
967                        zcl=cld_lc_con
968                        zct=1./cld_tau_con
969                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i)*velo(i,k)
970                    ELSE
971                        zcl=cld_lc_lsc
972                        zct=1./cld_tau_lsc
973                        zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i) &   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
974                            *velo(i,k)
975                    ENDIF
976
977                    ! if vertical heterogeneity is taken into account, we use
978                    ! the "true" volume fraction instead of a modified
979                    ! surface fraction (which is larger and artificially
980                    ! reduces the in-cloud water).
981
982
983                    IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN
984                        zchau = zct   *dtime/REAL(ninter) * zoliql(i) &
985                        *(1.0-EXP(-(zoliql(i)/ctot_vol(i,k)/zcl)**2))
986                    ELSE
987                        zchau = zct   *dtime/REAL(ninter) * zoliql(i) &
988                        *(1.0-EXP(-(zoliql(i)/zneb(i)/zcl)**2))        ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2)
989                    ENDIF
990
991                    zrain   = MIN(MAX(zchau,0.0),zoliql(i))
992                    zsnow   = MIN(MAX(zfroi,0.0),zoliqi(i))
993                    zprecip = MAX(zrain + zsnow,0.0)
994
995                ENDIF
996
997                zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain  , 0.0)
998                zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow  , 0.0)
999                zoliq(i)  = MAX(zoliq(i)-zprecip , 0.0)
1000
1001                ! c_iso: call isotope_conversion (for readibility) or duplicate
1002
1003                radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(ninter+1)
1004                radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(ninter+1)
1005                radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(ninter+1)
1006
1007            ENDIF ! rneb >0
1008
1009        ENDDO  ! i = 1,klon
1010
1011    ENDDO      ! n = 1,niter
1012
1013    ! Precipitation flux calculation
1014
1015    DO i = 1, klon
1016           
1017            IF (iflag_evap_prec.EQ.4) THEN
1018                ziflprev(i)=ziflcld(i)
1019            ELSE
1020                ziflprev(i)=zifl(i)*zneb(i)
1021            ENDIF
1022
1023            IF (rneb(i,k) .GT. 0.0) THEN
1024
1025               ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
1026               ! If T<0C, liquid precip are converted into ice, which leads to an increase in
1027               ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
1028               ! taken into account through a linearization of the equations and by approximating
1029               ! the liquid precipitation process with a "threshold" process. We assume that
1030               ! condensates are not modified during this operation. Liquid precipitation is
1031               ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
1032               ! Water vapor increases as well
1033               ! Note that compared to fisrtilp, we always assume iflag_bergeron=2
1034
1035                    zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
1036                    zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
1037                    zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1038                    coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i)
1039                    ! Computation of DT if all the liquid precip freezes
1040                    DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1))
1041                    ! T should not exceed the freezing point
1042                    ! that is Delta > RTT-zt(i)
1043                    DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
1044                    zt(i)  = zt(i) + DeltaT
1045                    ! water vaporization due to temp. increase
1046                    Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT
1047                    ! we add this vaporized water to the total vapor and we remove it from the precip
1048                    zq(i) = zq(i) +  Deltaq
1049                    ! The three "max" lines herebelow protect from rounding errors
1050                    zcond(i) = max( zcond(i)- Deltaq, 0. )
1051                    ! liquid precipitation converted to ice precip
1052                    Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT
1053                    zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
1054                    ! iced water budget
1055                    zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
1056
1057                    ! c_iso : mv here condensation of isotopes + redispatchage en precip
1058
1059                IF (iflag_evap_prec.EQ.4) THEN
1060                    zrflcld(i) = zrflcld(i)+zqprecl(i) &
1061                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1062                    ziflcld(i) = ziflcld(i)+ zqpreci(i) &
1063                        *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1064                    znebprecipcld(i) = rneb(i,k)
1065                    zrfl(i) = zrflcld(i) + zrflclr(i)
1066                    zifl(i) = ziflcld(i) + ziflclr(i)
1067                ELSE
1068                    zrfl(i) = zrfl(i)+ zqprecl(i)        &
1069                    *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1070                    zifl(i) = zifl(i)+ zqpreci(i)        &
1071                   *(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
1072                ENDIF
1073                ! c_iso : same for isotopes
1074 
1075           ENDIF ! rneb>0
1076
1077   ENDDO
1078
1079    ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
1080    ! if iflag_evap_pre=4
1081    IF (iflag_evap_prec.EQ.4) THEN
1082
1083        DO i=1,klon
1084
1085            IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN
1086                znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
1087                (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min)))
1088            ELSE
1089                znebprecipclr(i)=0.0
1090            ENDIF
1091
1092            IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN
1093                znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
1094                (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min)))
1095            ELSE
1096                znebprecipcld(i)=0.0
1097            ENDIF
1098
1099        ENDDO
1100
1101    ENDIF
1102
1103    ! End of precipitation formation
1104    ! --------------------------------
1105
1106    ! Outputs:
1107    ! Precipitation fluxes at layer interfaces
1108    ! and temperature and water species tendencies
1109    DO i = 1, klon
1110        psfl(i,k)=zifl(i)
1111        prfl(i,k)=zrfl(i)
1112        d_ql(i,k) = (1-zfice(i))*zoliq(i)
1113        d_qi(i,k) = zfice(i)*zoliq(i)
1114        d_q(i,k) = zq(i) - q(i,k)
1115        ! c_iso: same for isotopes
1116        d_t(i,k) = zt(i) - t(i,k)
1117    ENDDO
1118
1119    ! Calculation of the concentration of condensates seen by the radiation scheme
1120    ! for liquid phase, we take radocondl
1121    ! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
1122    ! we recaulate radocondi to account for contributions from the precipitation flux
1123
1124    IF ((ok_radocond_snow) .AND. (k .LT. klev)) THEN
1125        ! for the solid phase (crystals + snowflakes)
1126        ! we recalculate radocondi by summing
1127        ! the ice content calculated in the mesh
1128        ! + the contribution of the non-evaporated snowfall
1129        ! from the overlying layer
1130        DO i=1,klon
1131           IF (ziflprev(i) .NE. 0.0) THEN
1132              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
1133           ELSE
1134              radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
1135           ENDIF
1136           radocond(i,k)=radocondl(i,k)+radocondi(i,k)
1137        ENDDO
1138    ENDIF
1139
1140    ! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
1141    DO i=1,klon
1142        IF (radocond(i,k) .GT. 0.) THEN
1143            radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
1144        ENDIF
1145    ENDDO
1146
1147    ! ----------------------------------------------------------------
1148    ! P4> Wet scavenging
1149    ! ----------------------------------------------------------------
1150
1151    !Scavenging through nucleation in the layer
1152
1153    DO i = 1,klon
1154       
1155        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1156            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1157        ELSE
1158            beta(i,k) = 0.
1159        ENDIF
1160
1161        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1162
1163        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1164
1165            IF (t(i,k) .GE. t_glace_min) THEN
1166                zalpha_tr = a_tr_sca(3)
1167            ELSE
1168                zalpha_tr = a_tr_sca(4)
1169            ENDIF
1170
1171            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1172            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
1173
1174            ! Nucleation with a factor of -1 instead of -0.5
1175            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
1176
1177        ENDIF
1178
1179    ENDDO
1180
1181    ! Scavenging through impaction in the underlying layer
1182
1183    DO kk = k-1, 1, -1
1184
1185        DO i = 1, klon
1186
1187            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1188
1189                IF (t(i,kk) .GE. t_glace_min) THEN
1190                    zalpha_tr = a_tr_sca(1)
1191                ELSE
1192                    zalpha_tr = a_tr_sca(2)
1193                ENDIF
1194
1195              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1196              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1197
1198            ENDIF
1199
1200        ENDDO
1201
1202    ENDDO
1203
1204    !--save some variables for ice supersaturation
1205    !
1206    DO i = 1, klon
1207        ! for memory
1208        rneb_seri(i,k) = rneb(i,k)
1209
1210        ! for diagnostics
1211        rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
1212
1213        qvc(i,k) = zqs(i) * rneb(i,k)
1214        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
1215        qcld(i,k) = qvc(i,k) + zcond(i)
1216   ENDDO
1217   !q_sat
1218   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:))
1219   CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
1220
1221ENDDO
1222
1223!======================================================================
1224!                      END OF VERTICAL LOOP
1225!======================================================================
1226
1227  ! Rain or snow at the surface (depending on the first layer temperature)
1228  DO i = 1, klon
1229      snow(i) = zifl(i)
1230      rain(i) = zrfl(i)
1231      ! c_iso final output
1232  ENDDO
1233
1234  IF (ncoreczq>0) THEN
1235      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1236  ENDIF
1237
1238END SUBROUTINE lscp
1239!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1240
1241END MODULE lscp_mod
Note: See TracBrowser for help on using the repository browser.