source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90 @ 4910

Last change on this file since 4910 was 4910, checked in by evignon, 4 weeks ago

debut du chantier de refonte de cloudth_mpc

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