source: LMDZ6/branches/Test_modipsl/libf/phylmd/lscp_mod.F90 @ 5440

Last change on this file since 5440 was 4563, checked in by evignon, 19 months ago

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