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

Last change on this file since 5422 was 5203, checked in by Laurent Fairhead, 3 months ago

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