source: LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_main.f90

Last change on this file was 5797, checked in by aborella, 2 days ago

Bugfix for saturation adjustment in cirrus mixing + bugfix for contrails sedimentation + new diagnostics + support for unadjusted contrails
File size: 74.5 KB
Line 
1!$gpum horizontal klon ngrid
2MODULE lmdz_lscp_main
3
4IMPLICIT NONE
5
6CONTAINS
7
8!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9SUBROUTINE lscp(klon, klev, dtime, missing_val,         &
10     paprs, pplay, omega, temp, qt, ql_seri, qi_seri,   &
11     ratqs, sigma_qtherm, ptconv, cfcon_old, qvcon_old, &
12     qccon_old, cfcon, qvcon, qccon,                    &
13     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
14     pfraclr, pfracld,                                  &
15     cldfraliq, cldfraliqth,                            &
16     sigma2_icefracturb,sigma2_icefracturbth,           &
17     mean_icefracturb,mean_icefracturbth,               &
18     radocond, radicefrac, rain, snow,                  &
19     frac_impa, frac_nucl, beta,                        &
20     prfl, psfl, rhcl, qta, fraca,                      &
21     tv, pspsk, tla, thl, wth, iflag_cld_th,            &
22     iflag_ice_thermo, distcltop, temp_cltop,           &
23     tke, tke_dissip,                                   &
24     entr_therm, detr_therm,                            &
25     cell_area, stratomask,                             &
26     cf_seri, qvc_seri, u_seri, v_seri,                 &
27     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
28     dcf_sub, dcf_con, dcf_mix,                         &
29     dqi_sed, dcf_sed, dqvc_sed,                        &
30     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
31     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
32     cfc_seri, qtc_seri, qic_seri, nic_seri,            &
33     qice_cont, flight_dist, flight_fuel,               &
34     contfrarad, qradice_cont,                          &
35     Tcritcont, qcritcont, potcontfraP, potcontfraNP,   &
36     cloudth_sth,                                       &
37     cloudth_senv, cloudth_sigmath, cloudth_sigmaenv,   &
38     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
39     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
40     dqsmelt, dqsfreez)
41
42!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
43! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
44!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
45!--------------------------------------------------------------------------------
46! Date: 01/2021
47!--------------------------------------------------------------------------------
48! Aim: Large Scale Clouds and Precipitation (LSCP)
49!
50! This code is a new version of the fisrtilp.F90 routine, which is itself a
51! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
52! and 'ilp' (il pleut, L. Li)
53!
54! Compared to the original fisrtilp code, lscp
55! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
56! -> consider always precipitation thermalisation (fl_cor_ebil>0)
57! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
58! -> option "oldbug" by JYG has been removed
59! -> iflag_t_glace >0 always
60! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
61! -> rectangular distribution from L. Li no longer available
62! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
63!--------------------------------------------------------------------------------
64! References:
65!
66! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
67! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
68! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
69! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
70! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
71! - Touzze-Peifert Ludo, PhD thesis, p117-124
72! -------------------------------------------------------------------------------
73! Code structure:
74!
75! P0>     Thermalization of the precipitation coming from the overlying layer
76! P1>     Evaporation of the precipitation (falling from the k+1 level)
77! P2>     Cloud formation (at the k level)
78! P2.A.1> With the PDFs, calculation of cloud properties using the inital
79!         values of T and Q
80! P2.A.2> Coupling between condensed water and temperature
81! P2.A.3> Calculation of final quantities associated with cloud formation
82! P2.B>   Release of Latent heat after cloud formation
83! P3>     Autoconversion to precipitation (k-level)
84! P4>     Wet scavenging
85!------------------------------------------------------------------------------
86! Some preliminary comments (JBM) :
87!
88! The cloud water that the radiation scheme sees is not the same that the cloud
89! water used in the physics and the dynamics
90!
91! During the autoconversion to precipitation (P3 step), radocond (cloud water used
92! by the radiation scheme) is calculated as an average of the water that remains
93! in the cloud during the precipitation and not the water remaining at the end
94! of the time step. The latter is used in the rest of the physics and advected
95! by the dynamics.
96!
97! In summary:
98!
99! Radiation:
100! xflwc(newmicro)+xfiwc(newmicro) =
101!   radocond=lwcon(nc)+iwcon(nc)
102!
103! Notetheless, be aware of:
104!
105! radocond .NE. ocond(nc)
106! i.e.:
107! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
108! but oliq+(ocond-oliq) .EQ. ocond
109! (which is not trivial)
110!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111!
112
113! USE de modules contenant des fonctions.
114USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
115USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
116USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
117USE lmdz_lscp_condensation, ONLY : cloudth, cloudth_v3, cloudth_v6, condensation_cloudth
118USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
119USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
120
121! Use du module lmdz_lscp_ini contenant les constantes
122USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
123USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
124USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
125USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
126USE lmdz_lscp_ini, ONLY : min_frac_th_cld, temp_nowater, rho_ice
127USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG, RPI
128USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
129USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
130USE lmdz_lscp_ini, ONLY : ok_weibull_warm_clouds, ok_no_issr_strato, ok_ice_sedim
131USE lmdz_lscp_ini, ONLY : ok_plane_contrail, ok_unadjusted_contrails
132USE lmdz_lscp_ini, ONLY : ok_nodeep_lscp, ok_higher_cirrus_cover
133USE lmdz_lscp_ini, ONLY : ok_lscp_mergecond, gamma_mixth
134USE lmdz_lscp_ini, ONLY : eff2vol_radius_contrails
135
136! Temporary call for Lamquin et al (2012) diagnostics
137USE phys_local_var_mod, ONLY : issrfra100to150, issrfra150to200, issrfra200to250
138USE phys_local_var_mod, ONLY : issrfra250to300, issrfra300to400, issrfra400to500
139USE phys_local_var_mod, ONLY : AEI_contrails, AEI_surv_contrails
140USE phys_local_var_mod, ONLY : fsurv_contrails, section_contrails
141USE phys_local_var_mod, ONLY : dcfc_ini, dqic_ini, dqtc_ini, dnic_ini
142USE phys_local_var_mod, ONLY : dcfc_sub, dqic_sub, dqtc_sub, dnic_sub
143USE phys_local_var_mod, ONLY : dcfc_mix, dqic_mix, dqtc_mix, dnic_mix
144USE phys_local_var_mod, ONLY : dnic_agg, dqic_adj, dqtc_adj
145USE phys_local_var_mod, ONLY : dcfc_sed, dqic_sed, dqtc_sed, dnic_sed
146USE phys_local_var_mod, ONLY : dcfc_auto, dqic_auto, dqtc_auto, dnic_auto
147USE phys_local_var_mod, ONLY : dcf_auto, dqi_auto, dqvc_auto
148USE phys_local_var_mod, ONLY : nice_ygcont, iwc_ygcont, rvol_ygcont, tau_ygcont
149USE phys_local_var_mod, ONLY : nice_vscont, iwc_vscont, rvol_vscont, tau_vscont
150USE phys_local_var_mod, ONLY : nice_cont, iwc_cont, rvol_cont, tau_cont
151
152IMPLICIT NONE
153
154!===============================================================================
155! VARIABLES DECLARATION
156!===============================================================================
157
158! INPUT VARIABLES:
159!-----------------
160
161  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
162  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
163 
164
165  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
166  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
167
168  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
169  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
170  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
171  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
172  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qt              ! total specific humidity (in vapor phase in input) [kg/kg]
173  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: ql_seri         ! liquid specific content [kg/kg]
174  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qi_seri         ! ice specific content [kg/kg]
175  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke             ! turbulent kinetic energy [m2/s2]
176  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke_dissip      ! TKE dissipation [m2/s3]
177  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: entr_therm      ! thermal plume entrainment rate * dz [kg/s/m2]
178  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: detr_therm      ! thermal plume detrainment rate * dz [kg/s/m2]
179
180
181 
182  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
183  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon_old       ! cloud fraction from deep convection from previous timestep [-]
184  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qvcon_old       ! in-cloud vapor specific humidity from deep convection from previous timestep [kg/kg]
185  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qccon_old       ! in-cloud condensed specific humidity from deep convection from previous timestep [kg/kg]
186  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: cfcon           ! cloud fraction from deep convection [-]
187  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qvcon           ! in-cloud vapor specific humidity from deep convection [kg/kg]
188  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qccon           ! in-cloud condensed specific humidity from deep convection [kg/kg]
189
190  !Inputs associated with thermal plumes
191
192  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
193  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
194  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
195  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
196  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
197  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid potential temperature within thermals [K]
198  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: wth                 ! vertical velocity within thermals [m/s]
199  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: sigma_qtherm        ! controls the saturation deficit distribution width in thermals [-]
200
201
202
203  ! INPUT/OUTPUT variables
204  !------------------------
205 
206  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
207  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs        ! function that sets the large-scale water distribution width
208
209
210  ! INPUT/OUTPUT condensation and ice supersaturation
211  !--------------------------------------------------
212  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
213  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qvc_seri         ! cloudy water vapor [kg/kg]
214  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
215  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
216  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
217  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: stratomask       ! fraction of stratosphere (0 or 1)
218
219  ! INPUT/OUTPUT aviation
220  !--------------------------------------------------
221  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cfc_seri         ! contrail fraction [-]
222  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qtc_seri         ! contrail total specific humidity [kg/kg]
223  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: qic_seri         ! contrail ice specific humidity [kg/kg]
224  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: nic_seri         ! contrail ice crystals concentration [#/kg]
225  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! aviation distance flown concentration [m/s/m3]
226  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_fuel      ! aviation fuel consumption concentration [kg/s/m3]
227
228  ! OUTPUT variables
229  !-----------------
230
231  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
238  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
239  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq        ! liquid fraction of cloud fraction [-]
240  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliqth      ! liquid fraction of cloud fraction in thermals [-]
241  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
242  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturbth ! Variance of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturbth   ! Mean of the diagnostic supersaturation distribution in thermals (icefrac_turb) [-]
245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
246  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
247  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
248  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
249  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
250  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
251  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
252  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
255
256  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
257 
258  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
259  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
260 
261  ! for condensation and ice supersaturation
262
263  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
264  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
265  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
266  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
267  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
268  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
269  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
270  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
271  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
272  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
273  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
274  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
275  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
276  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
277  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
278  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
279  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
280  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sed        !--ice water content tendency due to sedmentation of ice crystals [kg/kg/s]
281  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sed        !--cloud fraction tendency due to sedimentation of ice crystals [kg/kg/s]
282  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sed       !--cloud water vapor tendency due to sedimentation of ice crystals [kg/kg/s]
283  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
284  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
285
286  ! for contrails and aviation
287
288  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qice_cont      !--condensed water in contrails [kg/kg]
289  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: contfrarad     !--contrail fraction for radiation [-]
290  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qradice_cont   !--condensed water in contrails used in the radiation scheme [kg/kg]
291  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcritcont      !--critical temperature for contrail formation [K]
292  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcritcont      !--critical specific humidity for contrail formation [kg/kg]
293  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: potcontfraP    !--potential persistent contrail fraction [-]
294  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: potcontfraNP   !--potential non-persistent contrail fraction [-]
295
296
297  ! for POPRECIP
298
299  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
300  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
301  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
302  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
303  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
304  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
305  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
306  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
307  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
308  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
309  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
310  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
311  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
312
313  ! for thermals
314
315  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
316  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
317  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
318  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
319
320
321  ! LOCAL VARIABLES:
322  !----------------
323  REAL, DIMENSION(klon) :: qliq_in, qice_in, qvc_in, cldfra_in
324  REAL, DIMENSION(klon,klev) :: ctot, rnebth, ctot_vol
325  REAL, DIMENSION(klon,klev) :: wls                                 !-- large scalce vertical velocity [m/s]
326  REAL, DIMENSION(klon) :: zqs, zdqs, zqsl, zdqsl, zqsi, zdqsi
327  REAL, DIMENSION(klon) :: zqsth, zqslth, zqsith, zdqsth, zdqslth, zdqsith
328  REAL :: zdelta
329  REAL, DIMENSION(klon) :: zdqsdT_raw
330  REAL, DIMENSION(klon) :: gammasat,dgammasatdt                   ! coefficient to make cold condensation at the correct RH and derivative wrt T
331  REAL, DIMENSION(klon) :: Tbef,Tbefth,Tbefthm1,qlibef,DT                  ! temperature, humidity and temp. variation during condensation iteration
332  REAL :: num,denom
333  REAL :: cste
334  REAL, DIMENSION(klon) :: qincloud
335  REAL, DIMENSION(klon) :: zrfl, zifl
336  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqth, zqn, zqnl
337  REAL, DIMENSION(klon) :: zoliql, zoliqi
338  REAL, DIMENSION(klon) :: zt, zp
339  REAL, DIMENSION(klon) :: zfice, zficeth, zficeenv, zneb, zcf, zsnow
340  REAL, DIMENSION(klon) :: dzfice, dzficeth, dzficeenv
341  REAL, DIMENSION(klon) :: qtot, zeroklon
342  ! Variables precipitation energy conservation
343  REAL, DIMENSION(klon) :: zmqc
344  REAL :: zalpha_tr
345  REAL :: zfrac_lessi
346  REAL, DIMENSION(klon) :: zprec_cond
347  REAL, DIMENSION(klon) :: zlh_solid
348  REAL, DIMENSION(klon) :: ztupnew
349  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
350  REAL, DIMENSION(klon) :: zrflclr, zrflcld
351  REAL, DIMENSION(klon) :: ziflclr, ziflcld
352  REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
353  REAL, DIMENSION(klon) :: tot_zneb
354  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop, zdeltaz
355  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl, zqliqth, zqiceth, zqvapclth, sursat_e, invtau_e ! for icefrac_lscp_turb
356  ! for ice sedimentation
357  REAL, DIMENSION(klon) :: dzsed, flsed, cfsed
358  REAL, DIMENSION(klon) :: dzsed_abv, flsed_abv, cfsed_abv
359  REAL, DIMENSION(klon) :: flauto
360  REAL :: qice_sedim
361
362  ! for quantity of condensates seen by radiation
363  REAL, DIMENSION(klon) :: zradocond, zradoice
364  REAL, DIMENSION(klon) :: zrho_up, zvelo_up
365  REAL :: dz, coef_cover
366 
367  ! for condensation and ice supersaturation
368  REAL, DIMENSION(klon) :: qvc, qvcl, shear
369  REAL, DIMENSION(klon) :: zrneb_deep, zcond_deep
370  REAL :: delta_z, deepconv_coef
371  ! for contrails
372  REAL, DIMENSION(klon) :: contfra, qcont, qvcont, Ncont
373  REAL, DIMENSION(klon) :: totfra_in, qtot_in
374  LOGICAL, DIMENSION(klon) :: pt_pron_clds
375  REAL, DIMENSION(klon) :: dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont
376  REAL, DIMENSION(klon) :: dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv
377  REAL :: rho, rhodz, iwp_cont, rei_cont
378  !--for Lamquin et al 2012 diagnostics
379  REAL, DIMENSION(klon) :: issrfra100to150UP, issrfra150to200UP, issrfra200to250UP
380  REAL, DIMENSION(klon) :: issrfra250to300UP, issrfra300to400UP, issrfra400to500UP
381
382  INTEGER i, k, kk, iter
383  INTEGER, DIMENSION(klon) :: n_i
384  INTEGER ncoreczq
385  LOGICAL iftop
386
387  LOGICAL, DIMENSION(klon) :: stratiform_or_all_distrib,pticefracturb
388  LOGICAL, DIMENSION(klon) :: keepgoing
389
390  CHARACTER (len = 20) :: modname = 'lscp'
391  CHARACTER (len = 80) :: abort_message
392 
393
394!===============================================================================
395! INITIALISATION
396!===============================================================================
397
398! Few initial checks
399
400
401IF (iflag_fisrtilp_qsat .LT. 0) THEN
402    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
403    CALL abort_physic(modname,abort_message,1)
404ENDIF
405
406! AA for 'safety' reasons
407zalpha_tr   = 0.
408zfrac_lessi = 0.
409beta(:,:)= 0.
410
411! Initialisation of variables:
412
413prfl(:,:) = 0.0
414psfl(:,:) = 0.0
415d_t(:,:) = 0.0
416d_q(:,:) = 0.0
417d_ql(:,:) = 0.0
418d_qi(:,:) = 0.0
419rneb(:,:) = 0.0
420rnebth(:,:)=0.0
421pfraclr(:,:)=0.0
422pfracld(:,:)=0.0
423cldfraliq(:,:)=0.
424sigma2_icefracturb(:,:)=0.
425mean_icefracturb(:,:)=0.
426cldfraliqth(:,:)=0.
427sigma2_icefracturbth(:,:)=0.
428mean_icefracturbth(:,:)=0.
429radocond(:,:) = 0.0
430radicefrac(:,:) = 0.0
431frac_nucl(:,:) = 1.0
432frac_impa(:,:) = 1.0
433rain(:) = 0.0
434snow(:) = 0.0
435zrfl(:) = 0.0
436zifl(:) = 0.0
437zneb(:) = seuil_neb
438zrflclr(:) = 0.0
439ziflclr(:) = 0.0
440zrflcld(:) = 0.0
441ziflcld(:) = 0.0
442tot_zneb(:) = 0.0
443zeroklon(:) = 0.0
444zdistcltop(:)=0.0
445ztemp_cltop(:) = 0.0
446ztupnew(:)=0.0
447ctot_vol(:,:)=0.0
448rneblsvol(:,:)=0.0
449znebprecip(:)=0.0
450znebprecipclr(:)=0.0
451znebprecipcld(:)=0.0
452distcltop(:,:)=0.
453temp_cltop(:,:)=0.
454
455
456!--Ice supersaturation
457gamma_cond(:,:) = 1.
458qissr(:,:)      = 0.
459issrfra(:,:)    = 0.
460dcf_sub(:,:)    = 0.
461dcf_con(:,:)    = 0.
462dcf_mix(:,:)    = 0.
463dcf_sed(:,:)    = 0.
464dcf_auto(:,:)   = 0.
465dqi_adj(:,:)    = 0.
466dqi_sub(:,:)    = 0.
467dqi_con(:,:)    = 0.
468dqi_mix(:,:)    = 0.
469dqi_sed(:,:)    = 0.
470dqi_auto(:,:)   = 0.
471dqvc_adj(:,:)   = 0.
472dqvc_sub(:,:)   = 0.
473dqvc_con(:,:)   = 0.
474dqvc_mix(:,:)   = 0.
475dqvc_sed(:,:)   = 0.
476dqvc_auto(:,:)  = 0.
477qvc(:)          = 0.
478shear(:)        = 0.
479flsed(:)        = 0.
480flauto(:)       = 0.
481flsed_cont(:)   = 0.
482Nflsed_cont(:)  = 0.
483pt_pron_clds(:) = .FALSE.
484
485!--Contrails
486dcfc_ini(:,:)  = 0.
487dqic_ini(:,:)  = 0.
488dqtc_ini(:,:)  = 0.
489dnic_ini(:,:)  = 0.
490dcfc_sub(:,:)  = 0.
491dqic_sub(:,:)  = 0.
492dqtc_sub(:,:)  = 0.
493dnic_sub(:,:)  = 0.
494dqic_adj(:,:)  = 0.
495dqtc_adj(:,:)  = 0.
496dcfc_mix(:,:)  = 0.
497dqic_mix(:,:)  = 0.
498dqtc_mix(:,:)  = 0.
499dnic_mix(:,:)  = 0.
500dnic_agg(:,:)  = 0.
501dcfc_sed(:,:)  = 0.
502dqic_sed(:,:)  = 0.
503dqtc_sed(:,:)  = 0.
504dnic_sed(:,:)  = 0.
505dcfc_auto(:,:) = 0.
506dqic_auto(:,:) = 0.
507dqtc_auto(:,:) = 0.
508dnic_auto(:,:) = 0.
509Tcritcont(:,:) = missing_val
510qcritcont(:,:) = missing_val
511potcontfraP(:,:)  = missing_val
512potcontfraNP(:,:) = missing_val
513AEI_contrails(:,:) = missing_val
514AEI_surv_contrails(:,:) = missing_val
515fsurv_contrails(:,:) = missing_val
516section_contrails(:,:) = missing_val
517
518
519!--for Lamquin et al (2012) diagnostics
520issrfra100to150(:)   = 0.
521issrfra100to150UP(:) = 0.
522issrfra150to200(:)   = 0.
523issrfra150to200UP(:) = 0.
524issrfra200to250(:)   = 0.
525issrfra200to250UP(:) = 0.
526issrfra250to300(:)   = 0.
527issrfra250to300UP(:) = 0.
528issrfra300to400(:)   = 0.
529issrfra300to400UP(:) = 0.
530issrfra400to500(:)   = 0.
531issrfra400to500UP(:) = 0.
532
533!-- poprecip
534qraindiag(:,:)= 0.
535qsnowdiag(:,:)= 0.
536dqreva(:,:)   = 0.
537dqrauto(:,:)  = 0.
538dqrmelt(:,:)  = 0.
539dqrfreez(:,:) = 0.
540dqrcol(:,:)   = 0.
541dqssub(:,:)   = 0.
542dqsauto(:,:)  = 0.
543dqsrim(:,:)   = 0.
544dqsagg(:,:)   = 0.
545dqsfreez(:,:) = 0.
546dqsmelt(:,:)  = 0.
547zqupnew(:)    = 0.
548zqvapclr(:)   = 0.
549
550!-- cloud phase useful variables
551wls(:,:)=-omega(:,:) / RG / (pplay(:,:)/RD/temp(:,:))
552zqliq(:)=0.
553zqice(:)=0.
554zqvapcl(:)=0.
555zqliqth(:)=0.
556zqiceth(:)=0.
557zqvapclth(:)=0.
558sursat_e(:)=0.
559invtau_e(:)=0.
560pticefracturb(:)=.FALSE.
561
562
563!===============================================================================
564!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
565!===============================================================================
566
567ncoreczq=0
568
569DO k = klev, 1, -1
570
571    IF (k.LE.klev-1) THEN
572       iftop=.false.
573    ELSE
574       iftop=.true.
575    ENDIF
576
577    ! Initialisation temperature and specific humidity
578    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
579    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
580    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
581    ! d_t = temperature tendency due to lscp
582    ! The temperature of the overlying layer is updated here because needed for thermalization
583    DO i = 1, klon
584        zt(i)=temp(i,k)
585        zq(i)=qt(i,k)
586        zp(i)=pplay(i,k)
587        qliq_in(i) = ql_seri(i,k)
588        qice_in(i) = qi_seri(i,k)
589        zcf(i) = 0.
590        zfice(i) = 1.0   ! initialized at 1 as by default we assume mpc to be at ice saturation
591        dzfice(i) = 0.0
592        zficeth(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
593        dzficeth(i) = 0.0
594        zficeenv(i) = 1.0 ! initialized at 1 as by default we assume mpc to be at ice saturation
595        dzficeenv(i) = 0.0
596       
597
598        IF (.not. iftop) THEN
599           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
600           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
601           !--zqs(i) is the saturation specific humidity in the layer above
602           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
603        ENDIF
604        !c_iso init of iso
605    ENDDO
606    IF ( ok_ice_supersat ) THEN
607      cldfra_in(:) = cf_seri(:,k)
608      qvc_in(:) = qvc_seri(:,k)
609    ENDIF
610
611    ! --------------------------------------------------------------------
612    ! P1> Precipitation processes, before cloud formation:
613    !     Thermalization of precipitation falling from the overlying layer AND
614    !     Precipitation evaporation/sublimation/melting
615    !---------------------------------------------------------------------
616   
617    !================================================================
618    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
619    IF ( ok_poprecip ) THEN
620
621      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
622                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
623                        zqvapclr, zqupnew, flsed, flsed_cont, &
624                        cldfra_in, qvc_in, qliq_in, qice_in, &
625                        zrfl, zrflclr, zrflcld, &
626                        zifl, ziflclr, ziflcld, &
627                        dqreva(:,k), dqssub(:,k) &
628                        )
629
630    !================================================================
631    ELSE
632
633      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
634                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
635                        flsed, flsed_cont, &
636                        zrfl, zrflclr, zrflcld, &
637                        zifl, ziflclr, ziflcld, &
638                        dqreva(:,k), dqssub(:,k) &
639                        )
640
641    ENDIF ! (ok_poprecip)
642   
643    ! Calculation of qsat,L/cp*dqsat/dT and ncoreczq counter
644    !-------------------------------------------------------
645
646     qtot(:)=zq(:)+zmqc(:)
647     CALL calc_qsat_ecmwf(klon,zt,qtot,zp,RTT,0,.false.,zqs,zdqs)
648     DO i = 1, klon
649            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
650            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
651            IF (zq(i) .LT. 1.e-15) THEN
652                ncoreczq=ncoreczq+1
653                zq(i)=1.e-15
654            ENDIF
655            ! c_iso: do something similar for isotopes
656
657     ENDDO
658   
659    ! -------------------------------------------------------------------------
660    ! P2> Cloud formation including condensation and cloud phase determination
661    !--------------------------------------------------------------------------
662    !
663    ! We always assume a 'fractional cloud' approach
664    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
665    ! is prescribed and depends on large scale variables and boundary layer
666    ! properties)
667    ! The decrease in condensed part due tu latent heating is taken into
668    ! account
669    ! -------------------------------------------------------------------
670
671        ! P2.1> With the PDFs (log-normal, bigaussian)
672        ! cloud properties calculation with the initial values of t and q
673        ! ----------------------------------------------------------------
674
675        ! initialise gammasat and stratiform_or_all_distrib
676        stratiform_or_all_distrib(:)=.TRUE.
677        gammasat(:)=1.
678
679        ! part of the code that is supposed to become obsolete soon
680        !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
681        IF (.NOT. ok_lscp_mergecond) THEN
682          IF (iflag_cld_th.GE.5) THEN
683               ! Cloud cover and content in meshes affected by shallow convection,
684               ! are retrieved from a bi-gaussian distribution of the saturation deficit
685               ! following Jam et al. 2013
686
687               IF (iflag_cloudth_vert.LE.2) THEN
688                  ! Old version of Arnaud Jam
689
690                    CALL cloudth(klon,klev,k,tv,                  &
691                         zq,qta,fraca,                            &
692                         qincloud,ctot,pspsk,paprs,pplay,tla,thl, &
693                         ratqs,zqs,temp,                              &
694                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
695
696
697                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
698                   ! Default version of Arnaud Jam
699                       
700                    CALL cloudth_v3(klon,klev,k,tv,                        &
701                         zq,qta,fraca,                                     &
702                         qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
703                         ratqs,sigma_qtherm,zqs,temp, &
704                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
705
706
707                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
708                   ! Jean Jouhaud's version, with specific separation between surface and volume
709                   ! cloud fraction Decembre 2018
710
711                    CALL cloudth_v6(klon,klev,k,tv,                        &
712                         zq,qta,fraca,                                     &
713                         qincloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
714                         ratqs,zqs,temp, &
715                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
716
717                ENDIF
718
719
720                DO i=1,klon
721                    rneb(i,k)=ctot(i,k)
722                    ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
723                    rneblsvol(i,k)=ctot_vol(i,k)
724                    zqn(i)=qincloud(i)
725                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
726                    qvc(i) = rneb(i,k) * zqs(i)
727                ENDDO
728
729               ! Cloud phase final determination for clouds after "old" cloudth calls
730               ! for those clouds, only temperature dependent phase partitioning (eventually modulated by
731               ! distance to cloud top) is available
732               ! compute distance to cloud top when cloud phase is computed depending on temperature
733               ! and distance to cloud top
734               IF (iflag_t_glace .GE. 4) THEN
735                    CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
736               ENDIF
737               CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
738
739          ENDIF
740
741          IF (iflag_cld_th .EQ. 5) THEN
742               ! When iflag_cld_th=5, we always assume
743                ! bi-gaussian distribution
744            stratiform_or_all_distrib(:) = .FALSE.
745
746          ELSEIF (iflag_cld_th .GE. 6) THEN
747                ! stratiform distribution only when no thermals
748            stratiform_or_all_distrib(:) = fraca(:,k) < min_frac_th_cld
749
750          ENDIF
751
752        ENDIF ! .not. ok_lscp_mergecond
753        !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
754
755
756        DO i = 1, klon
757          !--Calculate the shear value (input for condensation and ice supersat)
758          !--Cell thickness [m]
759          delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(i) * RD
760          IF ( iftop ) THEN
761            ! top
762            shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
763                           + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
764          ELSEIF ( k .EQ. 1 ) THEN
765            ! surface
766            shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
767                           + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
768          ELSE
769            ! other layers
770            shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
771                               - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
772                           + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
773                               - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
774          ENDIF
775        ENDDO
776
777        IF ( ok_ice_supersat ) THEN
778
779          !--Initialisation
780          IF ( ok_plane_contrail ) THEN
781            IF ( iftop ) THEN
782              dzsed_cont_abv(:) = 0.
783              flsed_cont_abv(:) = 0.
784              Nflsed_cont_abv(:)= 0.
785              cfsed_cont_abv(:) = 0.
786            ELSE
787              dzsed_cont_abv(:) = dzsed_cont(:)
788              flsed_cont_abv(:) = flsed_cont(:)
789              Nflsed_cont_abv(:)= Nflsed_cont(:)
790              cfsed_cont_abv(:) = cfsed_cont(:)
791            ENDIF
792            contfra(:) = 0.
793            qcont(:)   = 0.
794            qvcont(:)  = 0.
795            Ncont(:)   = 0.
796            dzsed_cont(:) = 0.
797            flsed_cont(:) = 0.
798            Nflsed_cont(:)= 0.
799            cfsed_cont(:) = 0.
800          ENDIF
801
802          IF ( iftop ) THEN
803            dzsed_abv(:) = 0.
804            flsed_abv(:) = 0.
805            cfsed_abv(:) = 0.
806          ELSE
807            dzsed_abv(:) = dzsed(:)
808            flsed_abv(:) = flsed(:)
809            cfsed_abv(:) = cfsed(:)
810          ENDIF
811          dzsed(:) = 0.
812          flsed(:) = 0.
813          cfsed(:) = 0.
814          flauto(:) = 0.
815
816          DO i = 1, klon
817            pt_pron_clds(i) = .TRUE.
818          ENDDO
819          IF ( .NOT. ok_weibull_warm_clouds ) THEN
820            DO i = 1, klon
821              pt_pron_clds(i) = pt_pron_clds(i) .AND. ( zt(i) .LE. temp_nowater )
822            ENDDO
823          ENDIF
824          IF ( ok_no_issr_strato ) THEN
825            DO i = 1, klon
826              pt_pron_clds(i) = pt_pron_clds(i) .AND. ( stratomask(i,k) .EQ. 0. )
827            ENDDO
828          ENDIF
829
830          totfra_in(:) = 1.
831          qtot_in(:) = zq(:)
832
833          IF ( ok_nodeep_lscp ) THEN
834            DO i = 1, klon
835              !--If deep convection is activated, the condensation scheme activates
836              !--only in the environment. NB. the clear sky fraction will the be
837              !--maximised by 1. - cfcon(i,k)
838              IF ( pt_pron_clds(i) .AND. ptconv(i,k) ) THEN
839                totfra_in(i) = 1. - cfcon(i,k)
840                qtot_in(i) = zq(i) - ( qvcon(i,k) + qccon(i,k) ) * cfcon(i,k)
841              ENDIF
842            ENDDO
843          ENDIF
844
845          DO i = 1, klon
846            IF ( pt_pron_clds(i) ) THEN
847              IF ( cfcon(i,k) .LT. cfcon_old(i,k) ) THEN
848                !--If deep convection is weakening, we add the clouds that are not anymore
849                !--'in' deep convection to the advected clouds
850                cldfra_in(i) = cldfra_in(i) + ( cfcon_old(i,k) - cfcon(i,k) )
851                qvc_in(i) = qvc_in(i) + qvcon_old(i,k) * ( cfcon_old(i,k) - cfcon(i,k) )
852                qice_in(i) = qice_in(i) + qccon_old(i,k) * ( cfcon_old(i,k) - cfcon(i,k) )
853              ELSEIF ( cldfra_in(i) .GT. eps ) THEN
854                !--Else if deep convection is strengthening, it consumes the existing cloud
855                !--fraction (which does not at this moment represent deep convection)
856                !deepconv_coef = 1. - ( cfcon(i,k) - cfcon_old(i,k) ) / ( 1. - cfcon_old(i,k) )
857                deepconv_coef = MAX(0., 1. - ( cfcon(i,k) - cfcon_old(i,k) ) / cldfra_in(i) )
858                cldfra_in(i) = cldfra_in(i) * deepconv_coef
859                qvc_in(i)    = qvc_in(i)    * deepconv_coef
860                qice_in(i)   = qice_in(i)   * deepconv_coef
861                IF ( ok_plane_contrail ) THEN
862                  !--If contrails are activated, their fraction is also reduced when deep
863                  !--convection is active
864                  cfc_seri(i,k) = cfc_seri(i,k) * deepconv_coef
865                  qtc_seri(i,k) = qtc_seri(i,k) * deepconv_coef
866                  qic_seri(i,k) = qic_seri(i,k) * deepconv_coef
867                  nic_seri(i,k) = nic_seri(i,k) * deepconv_coef
868                ENDIF
869              ENDIF
870            ENDIF
871          ENDDO
872        ENDIF
873
874
875        DT(:) = 0.
876        n_i(:)=0
877        Tbef(:)=zt(:)
878        qlibef(:)=0.
879        Tbefth(:)=tla(:,k)*pspsk(:,k)
880        IF (k .GT. 1) THEN
881         Tbefthm1(:)=tla(:,k-1)*pspsk(:,k-1)
882        ELSE
883         Tbefthm1(:)=Tbefth(:)       
884        ENDIF
885        zqth(:)=qta(:,k)
886        zdeltaz(:)=(paprs(:,k)-paprs(:,k+1))/RG/zp(:)*RD*zt(:)
887
888        ! Treatment of stratiform clouds (lognormale or ice-sursat) or all clouds (including cloudth
889        ! in case of ok_lscp_mergecond)
890        ! We iterate here to take into account the change in qsat(T) and ice fraction
891        ! during the condensation process
892        ! the increment in temperature is calculated using the first principle of
893        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
894        ! and a clear sky fraction)
895        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
896
897        DO iter=1,iflag_fisrtilp_qsat+1
898
899                keepgoing(:) = .FALSE.
900
901                DO i=1,klon
902
903                    ! keepgoing = .true. while convergence is not satisfied
904
905                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. stratiform_or_all_distrib(i)) THEN
906
907                        ! if not convergence:
908                        ! we calculate a new iteration
909                        keepgoing(i) = .TRUE.
910
911                        ! P2.2.1> cloud fraction and condensed water mass calculation
912                        ! Calculated variables:
913                        ! rneb : cloud fraction
914                        ! zqn : total water within the cloud
915                        ! zcond: mean condensed water within the mesh
916                        ! rhcl: clear-sky relative humidity
917                        !---------------------------------------------------------------
918
919                        ! new temperature that only serves in the iteration process:
920                        Tbef(i)=Tbef(i)+DT(i)
921
922                        ! total water including mass of precip
923                        qtot(i)=zq(i)+zmqc(i)
924
925                      ENDIF
926
927                  ENDDO
928
929                  ! Calculation of saturation specific humidity
930                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,0,.false.,zqs,zdqs)
931                  ! also in thermals
932                  CALL calc_qsat_ecmwf(klon,Tbefth,zqth,zp,RTT,0,.false.,zqsth,zdqsth)
933
934                  IF (iflag_icefrac .GE. 1) THEN
935                  ! consider a ice weighted qs to ensure that liquid clouds at T<0 have a consistent cloud fraction
936                  ! and cloud condensed water content. Principle from Dietlitcher et al. 2018, GMD
937                  ! we make this option works only for the physically-based and tke-dependent param from Raillard et al.
938                  ! (iflag_icefrac>=1)
939                      CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,1,.false.,zqsl,zdqsl)
940                      CALL calc_qsat_ecmwf(klon,Tbef,qtot,zp,RTT,2,.false.,zqsi,zdqsi)
941                      DO i=1,klon
942                           zqs(i)=zfice(i)*zqsi(i)+(1.-zfice(i))*zqsl(i)
943                           zdqs(i)=zfice(i)*zdqsi(i)+zqsi(i)*dzfice(i)+(1.-zfice(i))*zdqsl(i)-zqsl(i)*dzfice(i)
944                      ENDDO
945                  ENDIF
946                  IF (iflag_icefrac .GE. 2) THEN
947                  ! same adjustment for thermals
948                      CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,1,.false.,zqslth,zdqslth)
949                      CALL calc_qsat_ecmwf(klon,Tbefth,qtot,zp,RTT,2,.false.,zqsith,zdqsith)
950                      DO i=1,klon
951                           zqsth(i)=zficeth(i)*zqsith(i)+(1.-zficeth(i))*zqslth(i)
952                           zdqsth(i)=zficeth(i)*zdqsith(i)+zqsith(i)*dzficeth(i)+(1.-zficeth(i))*zdqslth(i)-zqslth(i)*dzficeth(i)
953                      ENDDO
954                  ENDIF
955
956                  CALL calc_gammasat(klon,Tbef,qtot,zp,gammasat,dgammasatdt)
957                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
958                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
959
960                  ! Cloud condensation based on subgrid pdf
961                  !--AB Activates a condensation scheme that allows for
962                  !--ice supersaturation and contrails evolution from aviation
963                  IF (ok_ice_supersat) THEN
964
965                    !---------------------------------------------
966                    !--   CONDENSATION AND ICE SUPERSATURATION  --
967                    !---------------------------------------------
968
969                    CALL condensation_ice_supersat( &
970                        klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), &
971                        totfra_in, cldfra_in, qvc_in, qliq_in, qice_in, &
972                        shear, tke_dissip(:,k), cell_area, Tbef, qtot_in, zqs, &
973                        gammasat, ratqs(:,k), keepgoing, pt_pron_clds, &
974                        dzsed_abv, flsed_abv, cfsed_abv, dzsed, flsed, cfsed, flauto, &
975                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
976                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
977                        dcf_sed(:,k), dcf_auto(:,k), &
978                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
979                        dqi_sed(:,k), dqi_auto(:,k), &
980                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
981                        dqvc_sed(:,k), dqvc_auto(:,k), &
982                        cfc_seri(:,k), qtc_seri(:,k), qic_seri(:,k), nic_seri(:,k), &
983                        flight_dist(:,k), flight_fuel(:,k), &
984                        contfra, qcont, qvcont, Ncont, &
985                        Tcritcont(:,k), qcritcont(:,k), potcontfraP(:,k), potcontfraNP(:,k), &
986                        AEI_contrails(:,k), AEI_surv_contrails(:,k), &
987                        fsurv_contrails(:,k), section_contrails(:,k), &
988                        dzsed_cont_abv, flsed_cont_abv, Nflsed_cont_abv, cfsed_cont_abv, &
989                        dzsed_cont, flsed_cont, Nflsed_cont, cfsed_cont, &
990                        dcfc_ini(:,k), dqic_ini(:,k), dqtc_ini(:,k), dnic_ini(:,k), &
991                        dcfc_sub(:,k), dqic_sub(:,k), dqtc_sub(:,k), dnic_sub(:,k), &
992                        dcfc_mix(:,k), dqic_mix(:,k), dqtc_mix(:,k), dnic_mix(:,k), &
993                        dnic_agg(:,k), dqic_adj(:,k), dqtc_adj(:,k), &
994                        dcfc_sed(:,k), dqic_sed(:,k), dqtc_sed(:,k), dnic_sed(:,k), &
995                        dcfc_auto(:,k), dqic_auto(:,k), dqtc_auto(:,k), dnic_auto(:,k))
996
997                  ELSE
998                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
999
1000                  CALL condensation_lognormal( &
1001                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
1002                       keepgoing, rneb(:,k), zqn, qvc)
1003
1004
1005                  ENDIF ! .NOT. ok_ice_supersat
1006
1007                  ! Volume cloud fraction
1008                  rneblsvol(:,k)=rneb(:,k)
1009
1010
1011                  IF  (ok_lscp_mergecond) THEN
1012                      ! in that case we overwrite rneb, rneblsvol and zqn for shallow convective clouds only
1013                      CALL condensation_cloudth(klon,Tbef,zq,zqth,fraca(:,k), &
1014                                     pspsk(:,k),zp,tla(:,k), &
1015                                     ratqs(:,k),sigma_qtherm(:,k),zqsth,zqs,zqn,rneb(:,k),rnebth(:,k),rneblsvol(:,k), &
1016                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
1017                  ENDIF
1018
1019
1020
1021                  ! Cloud phase determination
1022
1023
1024                  IF (iflag_icefrac .LE. 1) THEN
1025                     ! "old" phase partitioning depending on temperature and eventually distance to cloud top everywhere
1026                     IF (iflag_t_glace .GE. 4) THEN
1027                        CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
1028                     ENDIF
1029                     CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
1030                  ENDIF
1031
1032                  IF (iflag_icefrac .EQ. 1) THEN
1033                     ! phase partitioning depending on turbulence, vertical velocity and ice crystal microphysics
1034                     ! only in stratiform clouds in the mixed phase regime (Raillard et al. 2025)
1035                     ! it overwrites temperature-dependent phase partitioning except for convective boundary layer clouds
1036                     pticefracturb(:) = (fraca(:,k) < min_frac_th_cld) .AND. (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT)
1037                     DO i=1,klon
1038                        qincloud(i)=zqn(i)
1039                        zcf(i)=MIN(MAX(rneb(i,k), 0.),1.)
1040                        sursat_e(i) = 0.
1041                        invtau_e(i) = 0.
1042                     ENDDO
1043                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), qice_in,    &
1044                     ziflcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, &
1045                     zficeenv, dzficeenv, cldfraliq(:,k),sigma2_icefracturb(:,k),mean_icefracturb(:,k))                     
1046                     DO i=1,klon
1047                        IF (pticefracturb(i)) THEN
1048                          zfice(i)=zficeenv(i)
1049                          dzfice(i)=dzficeenv(i)
1050                        ENDIF
1051                     ENDDO
1052                   
1053
1054                  ELSEIF (iflag_icefrac .GE. 2) THEN
1055                     ! compute also phase partitioning also in thermal clouds (neglecting tke in thermals as first assumption)
1056                     ! moreover, given the upward velocity of thermals, we assume that precipitation falls in the environment
1057                     ! hence we assume that no snow falls in thermals.
1058                     ! we activate only in the mixed phase regime not to distrub the cirrus param at very cold T
1059                     pticefracturb(:) = (Tbef(:) .GT. temp_nowater) .AND. (Tbef(:) .LT. RTT)
1060                     !Thermals
1061                     DO i=1,klon
1062                        IF (fraca(i,k) .GT. min_frac_th_cld) THEN
1063                           zcf(i)=MIN(MAX(rnebth(i,k),0.), 1.)/fraca(i,k)
1064                           qincloud(i)=zqn(i)*fraca(i,k)
1065                        ELSE
1066                           zcf(i) = 0.
1067                           qincloud(i) = 0.
1068                        ENDIF
1069                        sursat_e(i)=cloudth_senv(i,k)/zqsi(i)
1070                        invtau_e(i)=gamma_mixth*MAX(entr_therm(i,k)-detr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i)
1071                     ENDDO
1072                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbefth, zp, paprs(:,k), paprs(:,k+1), wth(:,k), qice_in,  &
1073                     zeroklon, qincloud, zcf, zeroklon, zeroklon, sursat_e, invtau_e, zqliqth, zqvapclth, zqiceth, &
1074                     zficeth, dzficeth,cldfraliqth(:,k), sigma2_icefracturbth(:,k), mean_icefracturbth(:,k))
1075                     !Environment
1076                     DO i=1,klon
1077                        qincloud(i)=zqn(i)*(1.-fraca(i,k))
1078                        zcf(i)=MIN(MAX(rneb(i,k)-rnebth(i,k), 0.),1.)/(1.-fraca(i,k))
1079                        IF (k .GT. 1) THEN
1080                           ! evaluate the mixing sursaturation using saturation deficit at level below
1081                           ! as air pacels detraining into clouds have not (less) seen yet entrainement from above
1082                           sursat_e(i)=cloudth_sth(i,k-1)/(zqsith(i)+zdqsith(i)*RCPD/RLSTT*(Tbefthm1(i)-Tbefth(i)))
1083                           ! mixing is assumed to scales with intensity of net detrainment/entrainment rate (D/dz-E/dz) / rho
1084                           invtau_e(i)=gamma_mixth*MAX(detr_therm(i,k)-entr_therm(i,k),0.)*RD*Tbef(i)/zp(i)/zdeltaz(i)
1085                        ELSE
1086                           sursat_e(i)=0.
1087                           invtau_e(i)=0.
1088                        ENDIF
1089                     ENDDO
1090                     CALL icefrac_lscp_turb(klon, dtime, pticefracturb, Tbef, zp, paprs(:,k), paprs(:,k+1), wls(:,k), qice_in,    &
1091                     ziflcld, qincloud, zcf, tke(:,k), tke_dissip(:,k), sursat_e, invtau_e, zqliq, zqvapcl, zqice, &
1092                     zfice, dzfice, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
1093 
1094                    ! adjust zfice to account for condensates in thermals'fraction
1095                     DO i=1,klon
1096                        IF ( zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i) .GT. 0.) THEN
1097                          zfice(i)=MIN(1., (zqiceth(i)+zqice(i))/(zqliqth(i)+zqliq(i)+zqiceth(i)+zqice(i)))
1098                          dzfice(i)=0. ! dxice/dT=0. always when using icefrac_lscp_turb
1099                        ENDIF
1100                     ENDDO
1101
1102                  ENDIF ! iflag_icefrac
1103
1104
1105
1106
1107                  DO i=1,klon
1108                      IF (keepgoing(i)) THEN
1109
1110                       ! P2.2.2> Approximative calculation of temperature variation DT
1111                       !           due to condensation.
1112                       ! Calculated variables:
1113                       ! dT : temperature change due to condensation
1114                       !---------------------------------------------------------------
1115
1116               
1117                        IF (zfice(i).LT.1) THEN
1118                            cste=RLVTT
1119                        ELSE
1120                            cste=RLSTT
1121                        ENDIF
1122                       
1123                        IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN
1124                          !--AB We relax the saturation adjustment assumption
1125                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1126                          IF ( rneb(i,k) .GT. eps ) THEN
1127                            qlibef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
1128                          ELSE
1129                            qlibef(i) = 0.
1130                          ENDIF
1131                        ELSE
1132                          qlibef(i)=max(0.,zqn(i)-zqs(i))
1133                        ENDIF
1134
1135                        IF ( ok_ice_sedim ) THEN
1136                          qice_sedim = (flauto(i) + flsed(i) + flsed_cont(i)) &
1137                               / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
1138                          ! Add the ice that was sedimented, as it is not included in zqn
1139                          qlibef(i) = qlibef(i) + qice_sedim
1140                        ENDIF
1141
1142                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
1143                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlibef(i)
1144                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
1145                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
1146                              *qlibef(i)*dzfice(i)
1147                        ! here we update a provisory temperature variable that only serves in the iteration
1148                        ! process
1149                        DT(i)=num/denom
1150                        n_i(i)=n_i(i)+1
1151
1152                    ENDIF ! end keepgoing
1153
1154                ENDDO     ! end loop on i
1155
1156        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
1157
1158        ! P2.2> Final quantities calculation
1159        ! Calculated variables:
1160        !   rneb : cloud fraction
1161        !   zcond: mean condensed water in the mesh
1162        !   zqn  : mean water vapor in the mesh
1163        !   zfice: ice fraction in clouds
1164        !   zt   : temperature
1165        !   rhcl : clear-sky relative humidity
1166        ! ----------------------------------------------------------------
1167
1168
1169        ! Water vapor and condensed water update, subsequent latent heat exchange for each cloud type
1170        !---------------------------------------------------------------------------------------------
1171        DO i=1, klon
1172                ! Checks on rneb, rhcl and zqn
1173                IF (rneb(i,k) .LE. 0.0) THEN
1174                    zqn(i) = 0.0
1175                    rneb(i,k) = 0.0
1176                    zcond(i) = 0.0
1177                    rhcl(i,k)=zq(i)/zqs(i)
1178                ELSE IF (rneb(i,k) .GE. 1.0) THEN
1179                    zqn(i) = zq(i)
1180                    rneb(i,k) = 1.0
1181                    IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN
1182                      !--AB We relax the saturation adjustment assumption
1183                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1184                      zcond(i) = MAX(0., zqn(i) - qvc(i))
1185                    ELSE
1186                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
1187                    ENDIF
1188                    rhcl(i,k)=1.0
1189                ELSE
1190                    IF ( ok_unadjusted_clouds .OR. ok_unadjusted_contrails ) THEN
1191                      !--AB We relax the saturation adjustment assumption
1192     
1193                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
1194                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
1195                    ELSE
1196                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
1197                    ENDIF
1198
1199                    ! following line is very strange and probably wrong
1200                    rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
1201                    ! Correct calculation of clear-sky relative humidity. To activate once
1202                    ! issues related to zqn>zq in bi-gaussian clouds are corrected
1203                    !--Relative humidity (no unit) in clear sky, calculated as rh = q / qsat
1204                    !--This is slighly approximated, the actual formula is
1205                    !-- rh = q / qsat * (eps + (1-eps)*qsat) / (eps + (1-eps)*q)
1206                    !--Here, rh = (qtot - qincld * cldfra) / clrfra / qsat
1207                    !--where (qtot - qincld * cldfra) is the grid-mean clear sky water content
1208                    ! rhcl(i,k) = ( zq(i) - zqn(i) * rneb(i,k) ) / ( 1. - rneb(i,k) ) / zqs(i)
1209                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
1210
1211                ENDIF
1212
1213                ! water vapor update
1214                zq(i) = zq(i) - zcond(i)
1215
1216                IF ( ok_ice_sedim ) THEN
1217                  qice_sedim = (flauto(i) + flsed(i) + flsed_cont(i)) &
1218                      / ( paprs(i,k) - paprs(i,k+1) ) * RG * dtime
1219                  ! Remove the ice that was sedimented. As it is not included in zqn,
1220                  ! we only remove it from the total water
1221                  zq(i) = zq(i) - qice_sedim
1222                  ! Temperature update due to phase change (sedimented ice was condensed)
1223                  zt(i) = zt(i) + qice_sedim &
1224                        * RLSTT / RCPD / ( 1. + RVTMP2 * ( zq(i) + zmqc(i) + zcond(i) ) )
1225                ENDIF
1226                       
1227                ! temperature update due to phase change
1228                zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
1229                &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
1230                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
1231        ENDDO
1232
1233    ! ----------------------------------------------------------------
1234    ! P3> Precipitation processes, after cloud formation
1235    !     - precipitation formation, melting/freezing
1236    ! ----------------------------------------------------------------
1237
1238    ! Initiate the quantity of liquid and solid condensates
1239    ! Note that in the following, zcond is the total amount of condensates
1240    ! including newly formed precipitations (i.e., condensates formed by the
1241    ! condensation process), while zoliq is the total amount of condensates in
1242    ! the cloud (i.e., on which precipitation processes have applied)
1243    DO i = 1, klon
1244      zoliq(i)  = zcond(i)
1245      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
1246      zoliqi(i) = zcond(i) * zfice(i)
1247    ENDDO
1248
1249    !================================================================
1250    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
1251    IF (ok_poprecip) THEN
1252
1253      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), zp, &
1254                            ctot_vol(:,k), ptconv(:,k), &
1255                            zt, zq, zoliql, zoliqi, zfice, &
1256                            rneb(:,k), znebprecipclr, znebprecipcld, &
1257                            zrfl, zrflclr, zrflcld, &
1258                            zifl, ziflclr, ziflcld, &
1259                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
1260                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
1261                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
1262                            dqsmelt(:,k), dqsfreez(:,k) &
1263                            )
1264      DO i = 1, klon
1265          zoliq(i) = zoliql(i) + zoliqi(i)
1266      ENDDO
1267
1268    !================================================================
1269    ELSE
1270
1271      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), zp, &
1272                            ctot_vol(:,k), ptconv(:,k), pt_pron_clds, zdqsdT_raw, &
1273                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
1274                            rneb(:,k), znebprecipclr, znebprecipcld, &
1275                            zneb, tot_zneb, zrho_up, zvelo_up, &
1276                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
1277                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
1278                            )
1279
1280    ENDIF ! ok_poprecip
1281
1282    IF ( ok_ice_sedim ) THEN
1283      zifl(:) = zifl(:) + flauto(:)
1284      ziflcld(:) = ziflcld(:) + flauto(:)
1285    ENDIF
1286
1287    ! End of precipitation processes after cloud formation
1288    ! ----------------------------------------------------
1289
1290    !----------------------------------------------------------------------
1291    ! P4> Calculation of cloud condensates amount seen by radiative scheme
1292    !----------------------------------------------------------------------
1293
1294    DO i=1,klon
1295
1296      IF (ok_poprecip) THEN
1297        IF (ok_radocond_snow) THEN
1298          radocond(i,k) = zoliq(i)
1299          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
1300        ELSE
1301          radocond(i,k) = zoliq(i)
1302          zradoice(i) = zoliqi(i)
1303        ENDIF
1304      ELSE
1305        radocond(i,k) = zradocond(i)
1306      ENDIF
1307
1308      ! calculate the percentage of ice in "radocond" so cloud+precip seen
1309      ! by the radiation scheme
1310      IF (radocond(i,k) .GT. 0.) THEN
1311        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
1312      ENDIF
1313    ENDDO
1314
1315    ! ----------------------------------------------------------------
1316    ! P5> Wet scavenging
1317    ! ----------------------------------------------------------------
1318
1319    !Scavenging through nucleation in the layer
1320
1321    DO i = 1,klon
1322       
1323        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
1324            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
1325        ELSE
1326            beta(i,k) = 0.
1327        ENDIF
1328
1329        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
1330
1331        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1332
1333            IF (temp(i,k) .GE. temp_nowater) THEN
1334                zalpha_tr = a_tr_sca(3)
1335            ELSE
1336                zalpha_tr = a_tr_sca(4)
1337            ENDIF
1338
1339            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1340            frac_nucl(i,k)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1341
1342            ! Nucleation with a factor of -1 instead of -0.5
1343            zfrac_lessi = 1. - EXP(-zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1344
1345        ENDIF
1346
1347    ENDDO
1348
1349    ! Scavenging through impaction in the underlying layer
1350
1351    DO kk = k-1, 1, -1
1352
1353        DO i = 1, klon
1354
1355            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1356
1357                IF (temp(i,kk) .GE. temp_nowater) THEN
1358                    zalpha_tr = a_tr_sca(1)
1359                ELSE
1360                    zalpha_tr = a_tr_sca(2)
1361                ENDIF
1362
1363              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/MAX(rneb(i,k),seuil_neb))
1364              frac_impa(i,kk)= 1.-MAX(rneb(i,k),seuil_neb)*zfrac_lessi
1365
1366            ENDIF
1367
1368        ENDDO
1369
1370    ENDDO
1371   
1372    !------------------------------------------------------------
1373    ! P6 > write diagnostics and outputs
1374    !------------------------------------------------------------
1375
1376    CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,1,.false.,qsatl(:,k),zdqs)
1377    CALL calc_qsat_ecmwf(klon,zt,zeroklon,zp,RTT,2,.false.,qsati(:,k),zdqs)
1378   
1379    !--AB Write diagnostics and tracers for ice supersaturation
1380    IF ( ok_plane_contrail ) THEN
1381      DO i = 1, klon
1382        IF ( zoliq(i) .LE. 0. ) THEN
1383          contfra(i) = 0.
1384          qcont(i) = 0.
1385          qvcont(i) = 0.
1386          Ncont(i) = 0.
1387        ENDIF
1388      ENDDO
1389      cfc_seri(:,k) = contfra(:)
1390      nic_seri(:,k) = Ncont(:)
1391      !--Ice water content of contrails
1392      qice_cont(:,k) = qcont(:) - qvcont(:)
1393      IF ( ok_unadjusted_contrails ) THEN
1394        qtc_seri(:,k) = qvcont(:)
1395        qic_seri(:,k) = qice_cont(:,k)
1396      ELSE
1397        qtc_seri(:,k) = qcont(:)
1398      ENDIF
1399      !--Radiative properties
1400      contfrarad(:,k) = contfra(:)
1401      qradice_cont(:,k) = qice_cont(:,k)
1402    ENDIF
1403
1404    IF ( ok_ice_supersat ) THEN
1405      DO i = 1, klon
1406
1407        !--We save the cloud properties that will be advected
1408        cf_seri(i,k) = rneb(i,k)
1409        qvc_seri(i,k) = qvc(i)
1410
1411        !--We keep convective clouds properties in memory, and account for
1412        !--the sink of condensed water from precipitation
1413        IF ( ptconv(i,k) ) THEN
1414          qvcon_old(i,k) = qvcon(i,k)
1415          qccon_old(i,k) = qccon(i,k)
1416        ELSE
1417          qvcon_old(i,k) = 0.
1418          qccon_old(i,k) = 0.
1419        ENDIF
1420
1421        !--If everything was precipitated, the remaining empty cloud is dissipated
1422        !--and everything is transfered to the subsaturated clear sky region
1423        !--NB. we do not change rneb, as it is a diagnostic only
1424        IF ( zoliq(i) .LE. 0. ) THEN
1425          cf_seri(i,k) = 0.
1426          qvc_seri(i,k) = 0.
1427        ENDIF
1428
1429        !--Diagnostics
1430        gamma_cond(i,k) = gammasat(i)
1431        subfra(i,k) = totfra_in(i) - cf_seri(i,k) - issrfra(i,k)
1432        qsub(i,k) = qtot_in(i) - qvc_seri(i,k) - qissr(i,k)
1433        qcld(i,k) = qvc_seri(i,k) + zoliq(i)
1434
1435        IF ( ok_higher_cirrus_cover .AND. pt_pron_clds(i) .AND. .NOT. ptconv(i,k) ) THEN
1436          !--Following Brooks et al. (2005)
1437          !--This is only valid for cirrus clouds
1438          !--We do not apply it do convective clouds
1439          IF ( ( rneb(i,k) .GT. eps ) .AND. ( rneb(i,k) .LT. (1. - eps) ) ) THEN
1440            dz = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * zt(i) * RD
1441            coef_cover = (0.0706 + 0.1274 * shear(i)**0.3015) &
1442                * dz**0.7679 / SQRT(cell_area(i))**0.2254
1443            rneb(i,k) = 1. / (1. + EXP(-coef_cover) * (1. / rneb(i,k) - 1.))
1444            IF ( contfrarad(i,k) .GT. eps ) THEN
1445              contfrarad(i,k) = cfc_seri(i,k) / cf_seri(i,k) * rneb(i,k)
1446            ENDIF
1447          ENDIF
1448        ENDIF
1449
1450        !--Calculation of the ice supersaturated fraction following Lamquin et al (2012)
1451        !--methodology: in each layer, we make a maximum random overlap assumption for
1452        !--ice supersaturation
1453        IF ( ( paprs(i,k) .GT. 10000. ) .AND. ( paprs(i,k) .LE. 15000. ) ) THEN
1454                IF ( issrfra100to150UP(i) .GT. ( 1. - eps ) ) THEN
1455                        issrfra100to150(i) = 1.
1456                ELSE
1457                        issrfra100to150(i) = 1. - ( 1. - issrfra100to150(i) ) * &
1458                                ( 1. - MAX( issrfra(i,k), issrfra100to150UP(i) ) ) &
1459                              / ( 1. - issrfra100to150UP(i) )
1460                        issrfra100to150UP(i) = issrfra(i,k)
1461                ENDIF
1462        ELSEIF ( ( paprs(i,k) .GT. 15000. ) .AND. ( paprs(i,k) .LE. 20000. ) ) THEN
1463                IF ( issrfra150to200UP(i) .GT. ( 1. - eps ) ) THEN
1464                        issrfra150to200(i) = 1.
1465                ELSE
1466                        issrfra150to200(i) = 1. - ( 1. - issrfra150to200(i) ) * &
1467                                ( 1. - MAX( issrfra(i,k), issrfra150to200UP(i) ) ) &
1468                              / ( 1. - issrfra150to200UP(i) )
1469                        issrfra150to200UP(i) = issrfra(i,k)
1470                ENDIF
1471        ELSEIF ( ( paprs(i,k) .GT. 20000. ) .AND. ( paprs(i,k) .LE. 25000. ) ) THEN
1472                IF ( issrfra200to250UP(i) .GT. ( 1. - eps ) ) THEN
1473                        issrfra200to250(i) = 1.
1474                ELSE
1475                        issrfra200to250(i) = 1. - ( 1. - issrfra200to250(i) ) * &
1476                                ( 1. - MAX( issrfra(i,k), issrfra200to250UP(i) ) ) &
1477                              / ( 1. - issrfra200to250UP(i) )
1478                        issrfra200to250UP(i) = issrfra(i,k)
1479                ENDIF
1480        ELSEIF ( ( paprs(i,k) .GT. 25000. ) .AND. ( paprs(i,k) .LE. 30000. ) ) THEN
1481                IF ( issrfra250to300UP(i) .GT. ( 1. - eps ) ) THEN
1482                        issrfra250to300(i) = 1.
1483                ELSE
1484                        issrfra250to300(i) = 1. - ( 1. - issrfra250to300(i) ) * &
1485                                ( 1. - MAX( issrfra(i,k), issrfra250to300UP(i) ) ) &
1486                              / ( 1. - issrfra250to300UP(i) )
1487                        issrfra250to300UP(i) = issrfra(i,k)
1488                ENDIF
1489        ELSEIF ( ( paprs(i,k) .GT. 30000. ) .AND. ( paprs(i,k) .LE. 40000. ) ) THEN
1490                IF ( issrfra300to400UP(i) .GT. ( 1. - eps ) ) THEN
1491                        issrfra300to400(i) = 1.
1492                ELSE
1493                        issrfra300to400(i) = 1. - ( 1. - issrfra300to400(i) ) * &
1494                                ( 1. - MAX( issrfra(i,k), issrfra300to400UP(i) ) ) &
1495                              / ( 1. - issrfra300to400UP(i) )
1496                        issrfra300to400UP(i) = issrfra(i,k)
1497                ENDIF
1498        ELSEIF ( ( paprs(i,k) .GT. 40000. ) .AND. ( paprs(i,k) .LE. 50000. ) ) THEN
1499                IF ( issrfra400to500UP(i) .GT. ( 1. - eps ) ) THEN
1500                        issrfra400to500(i) = 1.
1501                ELSE
1502                        issrfra400to500(i) = 1. - ( 1. - issrfra400to500(i) ) * &
1503                                ( 1. - MAX( issrfra(i,k), issrfra400to500UP(i) ) ) &
1504                              / ( 1. - issrfra400to500UP(i) )
1505                        issrfra400to500UP(i) = issrfra(i,k)
1506                ENDIF
1507        ENDIF
1508
1509      ENDDO
1510    ENDIF
1511
1512    IF ( ok_plane_contrail ) THEN
1513      !--Other useful diagnostics
1514      DO i = 1, klon
1515        !--Very young countrails
1516        IF ( dcfc_ini(i,k) .GT. eps ) THEN
1517          rho = pplay(i,k) / zt(i) / RD
1518          nice_ygcont(i,k) = dnic_ini(i,k) / dcfc_ini(i,k) / 1e6 * rho
1519          iwc_ygcont(i,k) = dqic_ini(i,k) / dcfc_ini(i,k) * 1e3 * rho
1520          rvol_ygcont(i,k) = (dqic_ini(i,k) / MAX(eps, dnic_ini(i,k)) / rho_ice * 3. / 4. / RPI)**(1./3.) * 1e6
1521
1522          rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
1523          iwp_cont = 1e3 * dqic_ini(i,k) / dcfc_ini(i,k) * rhodz
1524          rei_cont = MIN(200., MAX(10., rvol_ygcont(i,k) / eff2vol_radius_contrails))
1525          tau_ygcont(i,k) = iwp_cont*(3.448e-3+2.431/rei_cont)
1526        ELSE
1527          nice_ygcont(i,k) = missing_val
1528          iwc_ygcont(i,k) = missing_val
1529          rvol_ygcont(i,k) = missing_val
1530          tau_ygcont(i,k) = missing_val
1531        ENDIF
1532        !--All contrails
1533        IF ( cfc_seri(i,k) .GT. 1e-8 ) THEN
1534          rho = pplay(i,k) / zt(i) / RD
1535          nice_cont(i,k) = nic_seri(i,k) / cfc_seri(i,k) / 1e6 * rho
1536          iwc_cont(i,k) = qice_cont(i,k) / cfc_seri(i,k) * 1e3 * rho
1537          rvol_cont(i,k) = (qice_cont(i,k) / MAX(eps, nic_seri(i,k)) / rho_ice * 3. / 4. / RPI)**(1./3.) * 1e6
1538
1539          rhodz = ( paprs(i,k) - paprs(i,k+1) ) / RG
1540          iwp_cont = 1e3 * qice_cont(i,k) / contfrarad(i,k) * rhodz
1541          rei_cont = MIN(200., MAX(10., rvol_cont(i,k) / eff2vol_radius_contrails))
1542          tau_cont(i,k) = iwp_cont*(3.448e-3+2.431/rei_cont)
1543          IF ( tau_cont(i,k) .GT. 0.05 ) THEN
1544            nice_vscont(i,k) = nice_cont(i,k)
1545            iwc_vscont(i,k) = iwc_cont(i,k)
1546            rvol_vscont(i,k) = rvol_cont(i,k)
1547            tau_vscont(i,k) = tau_cont(i,k)
1548          ELSE
1549            nice_vscont(i,k) = missing_val
1550            iwc_vscont(i,k) = missing_val
1551            rvol_vscont(i,k) = missing_val
1552            tau_vscont(i,k) = missing_val
1553          ENDIF
1554        ELSE
1555          nice_cont(i,k) = missing_val
1556          iwc_cont(i,k) = missing_val
1557          rvol_cont(i,k) = missing_val
1558          tau_cont(i,k) = missing_val
1559          nice_vscont(i,k) = missing_val
1560          iwc_vscont(i,k) = missing_val
1561          rvol_vscont(i,k) = missing_val
1562          tau_vscont(i,k) = missing_val
1563        ENDIF
1564      ENDDO
1565    ENDIF
1566
1567    ! Outputs:
1568    !-------------------------------
1569    ! Precipitation fluxes at layer interfaces
1570    ! + precipitation fractions +
1571    ! temperature and water species tendencies
1572    DO i = 1, klon
1573        psfl(i,k)=zifl(i)
1574        prfl(i,k)=zrfl(i)
1575        pfraclr(i,k)=znebprecipclr(i)
1576        pfracld(i,k)=znebprecipcld(i)
1577        distcltop(i,k)=zdistcltop(i)
1578        temp_cltop(i,k)=ztemp_cltop(i)
1579        d_q(i,k) = zq(i) - qt(i,k)
1580        d_t(i,k) = zt(i) - temp(i,k)
1581
1582        IF (ok_bug_phase_lscp) THEN
1583           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1584           d_qi(i,k) = zfice(i)*zoliq(i)
1585        ELSE
1586           d_ql(i,k) = zoliql(i)
1587           d_qi(i,k) = zoliqi(i)   
1588        ENDIF
1589
1590    ENDDO
1591
1592
1593ENDDO ! loop on k from top to bottom
1594
1595
1596  ! Rain or snow at the surface (depending on the first layer temperature)
1597  DO i = 1, klon
1598      snow(i) = zifl(i)
1599      rain(i) = zrfl(i)
1600      ! c_iso final output
1601  ENDDO
1602
1603  IF ( ok_ice_sedim ) THEN
1604    DO i = 1, klon
1605      snow(i) = snow(i) + flsed(i) + flsed_cont(i)
1606    ENDDO
1607  ENDIF
1608
1609  IF (ncoreczq>0) THEN
1610      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1611  ENDIF
1612
1613
1614RETURN
1615
1616END SUBROUTINE lscp
1617!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1618
1619END MODULE lmdz_lscp_main
Note: See TracBrowser for help on using the repository browser.