source: LMDZ6/branches/cirrus/libf/phylmd/lmdz_lscp.F90 @ 5185

Last change on this file since 5185 was 5163, checked in by aborella, 4 months ago

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