source: LMDZ6/trunk/libf/phylmd/lmdz_lscp.f90 @ 5433

Last change on this file since 5433 was 5430, checked in by evignon, 10 days ago

ajout de la deposition de vapeur sur les precipitations
glacees (croissance) dans poprecip, A borella
File size: 49.1 KB
Line 
1MODULE lmdz_lscp
2
3IMPLICIT NONE
4
5CONTAINS
6
7!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8SUBROUTINE lscp(klon,klev,dtime,missing_val,            &
9     paprs, pplay, omega, temp, qt, ql_seri, qi_seri,   &
10     ptconv, ratqs, sigma_qtherm,                       &
11     d_t, d_q, d_ql, d_qi, rneb, rneblsvol,             &
12     pfraclr, pfracld,                                  &
13     cldfraliq, sigma2_icefracturb,mean_icefracturb,    &
14     radocond, radicefrac, rain, snow,                  &
15     frac_impa, frac_nucl, beta,                        &
16     prfl, psfl, rhcl, qta, fraca,                      &
17     tv, pspsk, tla, thl, iflag_cld_th,                 &
18     iflag_ice_thermo, distcltop, temp_cltop,           &
19     tke, tke_dissip,                                   &
20     cell_area,                                         &
21     cf_seri, rvc_seri, u_seri, v_seri,                 &
22     qsub, qissr, qcld, subfra, issrfra, gamma_cond,    &
23     dcf_sub, dcf_con, dcf_mix,          &
24     dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj,      &
25     dqvc_sub, dqvc_con, dqvc_mix, qsatl, qsati,        &
26     Tcontr, qcontr, qcontr2, fcontrN, fcontrP, dcf_avi,&
27     dqi_avi, dqvc_avi, flight_dist, flight_h2o,        &
28     cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
29     qraindiag, qsnowdiag, dqreva, dqssub, dqrauto,     &
30     dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
31     dqsmelt, dqsfreez)
32
33!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
35!          A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
36!--------------------------------------------------------------------------------
37! Date: 01/2021
38!--------------------------------------------------------------------------------
39! Aim: Large Scale Clouds and Precipitation (LSCP)
40!
41! This code is a new version of the fisrtilp.F90 routine, which is itself a
42! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
43! and 'ilp' (il pleut, L. Li)
44!
45! Compared to the original fisrtilp code, lscp
46! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
47! -> consider always precipitation thermalisation (fl_cor_ebil>0)
48! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
49! -> option "oldbug" by JYG has been removed
50! -> iflag_t_glace >0 always
51! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
52! -> rectangular distribution from L. Li no longer available
53! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
54!--------------------------------------------------------------------------------
55! References:
56!
57! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
58! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
59! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
60! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
61! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
62! - Touzze-Peifert Ludo, PhD thesis, p117-124
63! -------------------------------------------------------------------------------
64! Code structure:
65!
66! P0>     Thermalization of the precipitation coming from the overlying layer
67! P1>     Evaporation of the precipitation (falling from the k+1 level)
68! P2>     Cloud formation (at the k level)
69! P2.A.1> With the PDFs, calculation of cloud properties using the inital
70!         values of T and Q
71! P2.A.2> Coupling between condensed water and temperature
72! P2.A.3> Calculation of final quantities associated with cloud formation
73! P2.B>   Release of Latent heat after cloud formation
74! P3>     Autoconversion to precipitation (k-level)
75! P4>     Wet scavenging
76!------------------------------------------------------------------------------
77! Some preliminary comments (JBM) :
78!
79! The cloud water that the radiation scheme sees is not the same that the cloud
80! water used in the physics and the dynamics
81!
82! During the autoconversion to precipitation (P3 step), radocond (cloud water used
83! by the radiation scheme) is calculated as an average of the water that remains
84! in the cloud during the precipitation and not the water remaining at the end
85! of the time step. The latter is used in the rest of the physics and advected
86! by the dynamics.
87!
88! In summary:
89!
90! Radiation:
91! xflwc(newmicro)+xfiwc(newmicro) =
92!   radocond=lwcon(nc)+iwcon(nc)
93!
94! Notetheless, be aware of:
95!
96! radocond .NE. ocond(nc)
97! i.e.:
98! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
99! but oliq+(ocond-oliq) .EQ. ocond
100! (which is not trivial)
101!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102!
103
104! USE de modules contenant des fonctions.
105USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
106USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
107USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top
108USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat
109USE lmdz_lscp_precip, ONLY : histprecip_precld, histprecip_postcld
110USE lmdz_lscp_precip, ONLY : poprecip_precld, poprecip_postcld
111
112! Use du module lmdz_lscp_ini contenant les constantes
113USE lmdz_lscp_ini, ONLY : prt_level, lunout, eps
114USE lmdz_lscp_ini, ONLY : seuil_neb, iflag_evap_prec, DDT0
115USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca
116USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_t_glace, iflag_fisrtilp_qsat
117USE lmdz_lscp_ini, ONLY : t_glace_min, min_frac_th_cld
118USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RVTMP2, RTT, RD, RG
119USE lmdz_lscp_ini, ONLY : ok_poprecip, ok_bug_phase_lscp
120USE lmdz_lscp_ini, ONLY : ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac
121
122IMPLICIT NONE
123
124!===============================================================================
125! VARIABLES DECLARATION
126!===============================================================================
127
128! INPUT VARIABLES:
129!-----------------
130
131  INTEGER,                         INTENT(IN)   :: klon,klev       ! number of horizontal grid points and vertical levels
132  REAL,                            INTENT(IN)   :: dtime           ! time step [s]
133  REAL,                            INTENT(IN)   :: missing_val     ! missing value for output
134
135  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: paprs           ! inter-layer pressure [Pa]
136  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pplay           ! mid-layer pressure [Pa]
137  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: temp            ! temperature (K)
138  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: omega           ! vertical velocity [Pa/s]
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)   :: ql_seri         ! liquid specific content [kg/kg]
141  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qi_seri         ! ice specific content [kg/kg]
142  INTEGER,                         INTENT(IN)   :: iflag_cld_th    ! flag that determines the distribution of convective clouds
143  INTEGER,                         INTENT(IN)   :: iflag_ice_thermo! flag to activate the ice thermodynamics
144                                                                   ! CR: if iflag_ice_thermo=2, only convection is active
145  LOGICAL, DIMENSION(klon,klev),   INTENT(IN)   :: ptconv          ! grid points where deep convection scheme is active
146
147  !Inputs associated with thermal plumes
148
149  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tv                  ! virtual potential temperature [K]
150  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: qta                 ! specific humidity within thermals [kg/kg]
151  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: fraca               ! fraction of thermals within the mesh [-]
152  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: pspsk               ! exner potential (p/100000)**(R/cp)
153  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: tla                 ! liquid temperature within thermals [K]
154  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke                 !--turbulent kinetic energy [m2/s2]
155  REAL, DIMENSION(klon,klev+1),    INTENT(IN)   :: tke_dissip          !--TKE dissipation [m2/s3]
156
157  ! INPUT/OUTPUT variables
158  !------------------------
159 
160  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: thl          ! liquid potential temperature [K]
161  REAL, DIMENSION(klon,klev),      INTENT(INOUT)   :: ratqs,sigma_qtherm        ! function of pressure that sets the large-scale
162
163
164  ! INPUT/OUTPUT condensation and ice supersaturation
165  !--------------------------------------------------
166  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: cf_seri          ! cloud fraction [-]
167  REAL, DIMENSION(klon,klev),      INTENT(INOUT):: rvc_seri         ! cloudy water vapor to total water vapor ratio [-]
168  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: u_seri           ! eastward wind [m/s]
169  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: v_seri           ! northward wind [m/s]
170  REAL, DIMENSION(klon),           INTENT(IN)   :: cell_area        ! area of each cell [m2]
171
172  ! INPUT/OUTPUT aviation
173  !--------------------------------------------------
174  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_dist      ! Aviation distance flown within the mesh [m/s/mesh]
175  REAL, DIMENSION(klon,klev),      INTENT(IN)   :: flight_h2o       ! Aviation H2O emitted within the mesh [kg H2O/s/mesh]
176 
177  ! OUTPUT variables
178  !-----------------
179
180  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_t              ! temperature increment [K]
181  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_q              ! specific humidity increment [kg/kg]
182  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_ql             ! liquid water increment [kg/kg]
183  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: d_qi             ! cloud ice mass increment [kg/kg]
184  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneb             ! cloud fraction [-]
185  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rneblsvol        ! cloud fraction per unit volume [-] 
186  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfraclr          ! precip fraction clear-sky part [-]
187  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: pfracld          ! precip fraction cloudy part [-]
188  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cldfraliq           ! liquid fraction of cloud [-]
189  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: sigma2_icefracturb  ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
190  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: mean_icefracturb    ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
191  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radocond         ! condensed water used in the radiation scheme [kg/kg]
192  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: radicefrac       ! ice fraction of condensed water for radiation scheme
193  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: rhcl             ! clear-sky relative humidity [-]
194  REAL, DIMENSION(klon),           INTENT(OUT)  :: rain             ! surface large-scale rainfall [kg/s/m2]
195  REAL, DIMENSION(klon),           INTENT(OUT)  :: snow             ! surface large-scale snowfall [kg/s/m2]
196  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: prfl             ! large-scale rainfall flux in the column [kg/s/m2]
197  REAL, DIMENSION(klon,klev+1),    INTENT(OUT)  :: psfl             ! large-scale snowfall flux in the column [kg/s/m2]
198  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: distcltop        ! distance to cloud top [m]
199  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: temp_cltop       ! temperature of cloud top [K]
200  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: beta             ! conversion rate of condensed water
201
202  ! fraction of aerosol scavenging through impaction and nucleation    (for on-line)
203 
204  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_impa           ! scavenging fraction due tu impaction [-]
205  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: frac_nucl           ! scavenging fraction due tu nucleation [-]
206 
207  ! for condensation and ice supersaturation
208
209  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsub           !--specific total water content in sub-saturated clear sky region [kg/kg]
210  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qissr          !--specific total water content in supersat region [kg/kg]
211  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcld           !--specific total water content in cloudy region [kg/kg]
212  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: subfra         !--mesh fraction of subsaturated clear sky [-]   
213  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: issrfra        !--mesh fraction of ISSR [-] 
214  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: gamma_cond     !--coefficient governing the ice nucleation RHi threshold [-]     
215  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_sub        !--cloud fraction tendency because of sublimation [s-1]
216  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_con        !--cloud fraction tendency because of condensation [s-1]
217  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_mix        !--cloud fraction tendency because of cloud mixing [s-1]
218  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_adj        !--specific ice content tendency because of temperature adjustment [kg/kg/s]
219  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_sub        !--specific ice content tendency because of sublimation [kg/kg/s]
220  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_con        !--specific ice content tendency because of condensation [kg/kg/s]
221  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_mix        !--specific ice content tendency because of cloud mixing [kg/kg/s]
222  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_adj       !--specific cloud water vapor tendency because of temperature adjustment [kg/kg/s]
223  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_sub       !--specific cloud water vapor tendency because of sublimation [kg/kg/s]
224  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_con       !--specific cloud water vapor tendency because of condensation [kg/kg/s]
225  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_mix       !--specific cloud water vapor tendency because of cloud mixing [kg/kg/s]
226  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsatl          !--saturation specific humidity wrt liquid [kg/kg]
227  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsati          !--saturation specific humidity wrt ice [kg/kg] 
228
229  ! for contrails and aviation
230
231  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: Tcontr         !--threshold temperature for contrail formation [K]
232  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr         !--threshold humidity for contrail formation [kg/kg]
233  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qcontr2        !--// (2nd expression more consistent with LMDZ expression of q)
234  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrN        !--fraction of grid favourable to non-persistent contrails
235  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: fcontrP        !--fraction of grid favourable to persistent contrails
236  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dcf_avi        !--cloud fraction tendency because of aviation [s-1]
237  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqi_avi        !--specific ice content tendency because of aviation [kg/kg/s]
238  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqvc_avi       !--specific cloud water vapor tendency because of aviation [kg/kg/s]
239
240
241  ! for POPRECIP
242
243  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qraindiag      !--DIAGNOSTIC specific rain content [kg/kg]
244  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: qsnowdiag      !--DIAGNOSTIC specific snow content [kg/kg]
245  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqreva         !--rain tendendy due to evaporation [kg/kg/s]
246  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqssub         !--snow tendency due to sublimation [kg/kg/s]
247  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrcol         !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
248  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsagg         !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
249  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrauto        !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
250  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsauto        !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
251  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsrim         !--snow tendency due to riming [kg/kg/s]
252  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsmelt        !--snow tendency due to melting [kg/kg/s]
253  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrmelt        !--rain tendency due to melting [kg/kg/s]
254  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqsfreez       !--snow tendency due to freezing [kg/kg/s]
255  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: dqrfreez       !--rain tendency due to freezing [kg/kg/s]
256
257  ! for thermals
258
259  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sth      !--mean saturation deficit in thermals
260  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_senv     !--mean saturation deficit in environment
261  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmath  !--std of saturation deficit in thermals
262  REAL, DIMENSION(klon,klev),      INTENT(OUT)  :: cloudth_sigmaenv !--std of saturation deficit in environment
263
264
265  ! LOCAL VARIABLES:
266  !----------------
267  REAL,DIMENSION(klon) :: qice_ini
268  REAL, DIMENSION(klon,klev) :: ctot
269  REAL, DIMENSION(klon,klev) :: ctot_vol
270  REAL, DIMENSION(klon) :: zqs, zdqs
271  REAL :: zdelta
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) :: zfice_th
278  REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
279  REAL, DIMENSION(klon) :: zrfl, zifl
280  REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
281  REAL, DIMENSION(klon) :: zoliql, zoliqi
282  REAL, DIMENSION(klon) :: zt
283  REAL, DIMENSION(klon) :: zfice,zneb
284  REAL, DIMENSION(klon) :: dzfice
285  REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
286  REAL, DIMENSION(klon) :: qtot, qzero
287  ! Variables precipitation energy conservation
288  REAL, DIMENSION(klon) :: zmqc
289  REAL :: zalpha_tr
290  REAL :: zfrac_lessi
291  REAL, DIMENSION(klon) :: zprec_cond
292  REAL, DIMENSION(klon) :: zlh_solid
293  REAL, DIMENSION(klon) :: ztupnew
294  REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
295  REAL, DIMENSION(klon) :: zrflclr, zrflcld
296  REAL, DIMENSION(klon) :: ziflclr, ziflcld
297  REAL, DIMENSION(klon) :: znebprecip, znebprecipclr, znebprecipcld
298  REAL, DIMENSION(klon) :: tot_zneb
299  REAL :: qlmpc, qimpc, rnebmpc
300  REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
301  REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl        ! for icefrac_lscp_turb
302
303  ! for quantity of condensates seen by radiation
304  REAL, DIMENSION(klon) :: zradocond, zradoice
305  REAL, DIMENSION(klon) :: zrho_up, zvelo_up
306 
307  ! for condensation and ice supersaturation
308  REAL, DIMENSION(klon) :: qvc, shear
309  REAL :: delta_z
310  !--Added for ice supersaturation (ok_ice_supersat) and contrails (ok_plane_contrails)
311  ! Constants used for calculating ratios that are advected (using a parent-child
312  ! formalism). This is not done in the dynamical core because at this moment,
313  ! only isotopes can use this parent-child formalism. Note that the two constants
314  ! are the same as the one use in the dynamical core, being also defined in
315  ! dyn3d_common/infotrac.F90
316  REAL :: min_qParent, min_ratio
317
318  INTEGER i, k, kk, iter
319  INTEGER, DIMENSION(klon) :: n_i
320  INTEGER ncoreczq
321  INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
322  LOGICAL iftop
323
324  LOGICAL, DIMENSION(klon) :: lognormale
325  LOGICAL, DIMENSION(klon) :: keepgoing
326
327  CHARACTER (len = 20) :: modname = 'lscp'
328  CHARACTER (len = 80) :: abort_message
329 
330
331!===============================================================================
332! INITIALISATION
333!===============================================================================
334
335! Few initial checks
336
337
338IF (iflag_fisrtilp_qsat .LT. 0) THEN
339    abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
340    CALL abort_physic(modname,abort_message,1)
341ENDIF
342
343! Few initialisations
344
345ctot_vol(1:klon,1:klev)=0.0
346rneblsvol(1:klon,1:klev)=0.0
347znebprecip(:)=0.0
348znebprecipclr(:)=0.0
349znebprecipcld(:)=0.0
350mpc_bl_points(:,:)=0
351
352IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
353
354! AA for 'safety' reasons
355zalpha_tr   = 0.
356zfrac_lessi = 0.
357beta(:,:)= 0.
358
359! Initialisation of variables:
360
361prfl(:,:) = 0.0
362psfl(:,:) = 0.0
363d_t(:,:) = 0.0
364d_q(:,:) = 0.0
365d_ql(:,:) = 0.0
366d_qi(:,:) = 0.0
367rneb(:,:) = 0.0
368pfraclr(:,:)=0.0
369pfracld(:,:)=0.0
370cldfraliq(:,:)=0.
371sigma2_icefracturb(:,:)=0.
372mean_icefracturb(:,:)=0.
373radocond(:,:) = 0.0
374radicefrac(:,:) = 0.0
375frac_nucl(:,:) = 1.0
376frac_impa(:,:) = 1.0
377rain(:) = 0.0
378snow(:) = 0.0
379zfice(:)=0.0
380dzfice(:)=0.0
381zfice_turb(:)=0.0
382dzfice_turb(:)=0.0
383zrfl(:) = 0.0
384zifl(:) = 0.0
385zneb(:) = seuil_neb
386zrflclr(:) = 0.0
387ziflclr(:) = 0.0
388zrflcld(:) = 0.0
389ziflcld(:) = 0.0
390tot_zneb(:) = 0.0
391qzero(:) = 0.0
392zdistcltop(:)=0.0
393ztemp_cltop(:) = 0.0
394ztupnew(:)=0.0
395
396distcltop(:,:)=0.
397temp_cltop(:,:)=0.
398
399!--Ice supersaturation
400gamma_cond(:,:) = 1.
401qissr(:,:)      = 0.
402issrfra(:,:)    = 0.
403dcf_sub(:,:)    = 0.
404dcf_con(:,:)    = 0.
405dcf_mix(:,:)    = 0.
406dqi_adj(:,:)    = 0.
407dqi_sub(:,:)    = 0.
408dqi_con(:,:)    = 0.
409dqi_mix(:,:)    = 0.
410dqvc_adj(:,:)   = 0.
411dqvc_sub(:,:)   = 0.
412dqvc_con(:,:)   = 0.
413dqvc_mix(:,:)   = 0.
414fcontrN(:,:)    = 0.
415fcontrP(:,:)    = 0.
416Tcontr(:,:)     = missing_val
417qcontr(:,:)     = missing_val
418qcontr2(:,:)    = missing_val
419dcf_avi(:,:)    = 0.
420dqi_avi(:,:)    = 0.
421dqvc_avi(:,:)   = 0.
422qvc(:)          = 0.
423shear(:)        = 0.
424min_qParent     = 1.e-30
425min_ratio       = 1.e-16
426
427!-- poprecip
428qraindiag(:,:)= 0.
429qsnowdiag(:,:)= 0.
430dqreva(:,:)   = 0.
431dqrauto(:,:)  = 0.
432dqrmelt(:,:)  = 0.
433dqrfreez(:,:) = 0.
434dqrcol(:,:)   = 0.
435dqssub(:,:)   = 0.
436dqsauto(:,:)  = 0.
437dqsrim(:,:)   = 0.
438dqsagg(:,:)   = 0.
439dqsfreez(:,:) = 0.
440dqsmelt(:,:)  = 0.
441zqupnew(:)    = 0.
442zqvapclr(:)   = 0.
443
444
445
446!c_iso: variable initialisation for iso
447
448
449!===============================================================================
450!              BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
451!===============================================================================
452
453ncoreczq=0
454
455DO k = klev, 1, -1
456
457    qice_ini = qi_seri(:,k)
458
459    IF (k.LE.klev-1) THEN
460       iftop=.false.
461    ELSE
462       iftop=.true.
463    ENDIF
464
465    ! Initialisation temperature and specific humidity
466    ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
467    ! at the end of the klon loop, a temperature incremtent d_t due to all processes
468    ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
469    ! d_t = temperature tendency due to lscp
470    ! The temperature of the overlying layer is updated here because needed for thermalization
471    DO i = 1, klon
472        zt(i)=temp(i,k)
473        zq(i)=qt(i,k)
474        IF (.not. iftop) THEN
475           ztupnew(i)  = temp(i,k+1) + d_t(i,k+1)
476           zqupnew(i)  = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
477           !--zqs(i) is the saturation specific humidity in the layer above
478           zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
479        ENDIF
480        !c_iso init of iso
481    ENDDO
482
483    ! --------------------------------------------------------------------
484    ! P1> Precipitation processes, before cloud formation:
485    !     Thermalization of precipitation falling from the overlying layer AND
486    !     Precipitation evaporation/sublimation/melting
487    !---------------------------------------------------------------------
488   
489    !================================================================
490    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
491    IF ( ok_poprecip ) THEN
492
493      CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
494                        zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
495                        zqvapclr, zqupnew, &
496                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
497                        zrfl, zrflclr, zrflcld, &
498                        zifl, ziflclr, ziflcld, &
499                        dqreva(:,k), dqssub(:,k) &
500                        )
501
502    !================================================================
503    ELSE
504
505      CALL histprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
506                        zt, ztupnew, zq, zmqc, zneb, znebprecip, znebprecipclr, &
507                        zrfl, zrflclr, zrflcld, &
508                        zifl, ziflclr, ziflcld, &
509                        dqreva(:,k), dqssub(:,k) &
510                        )
511
512    ENDIF ! (ok_poprecip)
513   
514    ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
515    !-------------------------------------------------------
516
517     qtot(:)=zq(:)+zmqc(:)
518     CALL calc_qsat_ecmwf(klon,zt,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
519     DO i = 1, klon
520            zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
521            zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta)
522            IF (zq(i) .LT. 1.e-15) THEN
523                ncoreczq=ncoreczq+1
524                zq(i)=1.e-15
525            ENDIF
526            ! c_iso: do something similar for isotopes
527
528     ENDDO
529   
530    ! --------------------------------------------------------------------
531    ! P2> Cloud formation
532    !---------------------------------------------------------------------
533    !
534    ! Unlike fisrtilp, we always assume a 'fractional cloud' approach
535    ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
536    ! is prescribed and depends on large scale variables and boundary layer
537    ! properties)
538    ! The decrease in condensed part due tu latent heating is taken into
539    ! account
540    ! -------------------------------------------------------------------
541
542        ! P2.1> With the PDFs (log-normal, bigaussian)
543        ! cloud properties calculation with the initial values of t and q
544        ! ----------------------------------------------------------------
545
546        ! initialise gammasat and qincloud_mpc
547        gammasat(:)=1.
548        qincloud_mpc(:)=0.
549
550        IF (iflag_cld_th.GE.5) THEN
551               ! Cloud cover and content in meshes affected by shallow convection,
552               ! are retrieved from a bi-gaussian distribution of the saturation deficit
553               ! following Jam et al. 2013
554
555               IF (iflag_cloudth_vert.LE.2) THEN
556                  ! Old version of Arnaud Jam
557
558                    CALL cloudth(klon,klev,k,tv,                  &
559                         zq,qta,fraca,                            &
560                         qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
561                         ratqs,zqs,temp,                              &
562                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
563
564
565                ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN
566                   ! Default version of Arnaud Jam
567                       
568                    CALL cloudth_v3(klon,klev,k,tv,                        &
569                         zq,qta,fraca,                                     &
570                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
571                         ratqs,sigma_qtherm,zqs,temp, &
572                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
573
574
575                ELSEIF (iflag_cloudth_vert.EQ.6) THEN
576                   ! Jean Jouhaud's version, with specific separation between surface and volume
577                   ! cloud fraction Decembre 2018
578
579                    CALL cloudth_v6(klon,klev,k,tv,                        &
580                         zq,qta,fraca,                                     &
581                         qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
582                         ratqs,zqs,temp, &
583                         cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
584
585                ELSEIF (iflag_cloudth_vert .EQ. 7) THEN
586                   ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
587                   ! for boundary-layer mixed phase clouds
588                    CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
589                                     pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
590                                     ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
591                                     cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
592
593                ENDIF
594
595
596                DO i=1,klon
597                    rneb(i,k)=ctot(i,k)
598                    rneblsvol(i,k)=ctot_vol(i,k)
599                    zqn(i)=qcloud(i)
600                    !--AB grid-mean vapor in the cloud - we assume saturation adjustment
601                    qvc(i) = rneb(i,k) * zqs(i)
602                ENDDO
603
604        ENDIF
605
606        IF (iflag_cld_th .LE. 4) THEN
607           
608                ! lognormal
609            lognormale(:) = .TRUE.
610
611        ELSEIF (iflag_cld_th .GE. 6) THEN
612           
613                ! lognormal distribution when no thermals
614            lognormale(:) = fraca(:,k) < min_frac_th_cld
615
616        ELSE
617                ! When iflag_cld_th=5, we always assume
618                ! bi-gaussian distribution
619            lognormale(:) = .FALSE.
620       
621        ENDIF
622
623        DT(:) = 0.
624        n_i(:)=0
625        Tbef(:)=zt(:)
626        qlbef(:)=0.
627
628        ! Treatment of non-boundary layer clouds (lognormale)
629        ! We iterate here to take into account the change in qsat(T) and ice fraction
630        ! during the condensation process
631        ! the increment in temperature is calculated using the first principle of
632        ! thermodynamics (enthalpy conservation equation in a mesh composed of a cloud fraction
633        ! and a clear sky fraction)
634        ! note that we assume that the vapor in the cloud is at saturation for this calculation     
635
636        DO iter=1,iflag_fisrtilp_qsat+1
637
638                keepgoing(:) = .FALSE.
639
640                DO i=1,klon
641
642                    ! keepgoing = .true. while convergence is not satisfied
643
644                    IF (((ABS(DT(i)).GT.DDT0) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN
645
646                        ! if not convergence:
647                        ! we calculate a new iteration
648                        keepgoing(i) = .TRUE.
649
650                        ! P2.2.1> cloud fraction and condensed water mass calculation
651                        ! Calculated variables:
652                        ! rneb : cloud fraction
653                        ! zqn : total water within the cloud
654                        ! zcond: mean condensed water within the mesh
655                        ! rhcl: clear-sky relative humidity
656                        !---------------------------------------------------------------
657
658                        ! new temperature that only serves in the iteration process:
659                        Tbef(i)=Tbef(i)+DT(i)
660
661                        ! Rneb, qzn and zcond for lognormal PDFs
662                        qtot(i)=zq(i)+zmqc(i)
663
664                      ENDIF
665
666                  ENDDO
667
668                  ! Calculation of saturation specific humidity and ice fraction
669                  CALL calc_qsat_ecmwf(klon,Tbef,qtot,pplay(:,k),RTT,0,.false.,zqs,zdqs)
670                  CALL calc_gammasat(klon,Tbef,qtot,pplay(:,k),gammasat,dgammasatdt)
671                  ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
672                  zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:)
673                  ! cloud phase determination
674                  IF (iflag_t_glace.GE.4) THEN
675                  ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
676                       CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
677                  ENDIF
678
679                  CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop,ztemp_cltop,zfice,dzfice)
680
681                  !--AB Activates a condensation scheme that allows for
682                  !--ice supersaturation and contrails evolution from aviation
683                  IF (ok_ice_supersat) THEN
684
685                    !--Calculate the shear value (input for condensation and ice supersat)
686                    DO i = 1, klon
687                      !--Cell thickness [m]
688                      delta_z = ( paprs(i,k) - paprs(i,k+1) ) / RG / pplay(i,k) * Tbef(i) * RD
689                      IF ( iftop ) THEN
690                        ! top
691                        shear(i) = SQRT( ( (u_seri(i,k) - u_seri(i,k-1)) / delta_z )**2. &
692                                       + ( (v_seri(i,k) - v_seri(i,k-1)) / delta_z )**2. )
693                      ELSEIF ( k .EQ. 1 ) THEN
694                        ! surface
695                        shear(i) = SQRT( ( (u_seri(i,k+1) - u_seri(i,k)) / delta_z )**2. &
696                                       + ( (v_seri(i,k+1) - v_seri(i,k)) / delta_z )**2. )
697                      ELSE
698                        ! other layers
699                        shear(i) = SQRT( ( ( (u_seri(i,k+1) + u_seri(i,k)) / 2. &
700                                           - (u_seri(i,k) + u_seri(i,k-1)) / 2. ) / delta_z )**2. &
701                                       + ( ( (v_seri(i,k+1) + v_seri(i,k)) / 2. &
702                                           - (v_seri(i,k) + v_seri(i,k-1)) / 2. ) / delta_z )**2. )
703                      ENDIF
704                    ENDDO
705
706                    !---------------------------------------------
707                    !--   CONDENSATION AND ICE SUPERSATURATION  --
708                    !---------------------------------------------
709
710                    CALL condensation_ice_supersat( &
711                        klon, dtime, missing_val, &
712                        pplay(:,k), paprs(:,k), paprs(:,k+1), &
713                        cf_seri(:,k), rvc_seri(:,k), ql_seri(:,k), qi_seri(:,k), &
714                        shear, tke_dissip(:,k), cell_area, &
715                        Tbef, zq, zqs, gammasat, ratqs(:,k), keepgoing, &
716                        rneb(:,k), zqn, qvc, issrfra(:,k), qissr(:,k), &
717                        dcf_sub(:,k), dcf_con(:,k), dcf_mix(:,k), &
718                        dqi_adj(:,k), dqi_sub(:,k), dqi_con(:,k), dqi_mix(:,k), &
719                        dqvc_adj(:,k), dqvc_sub(:,k), dqvc_con(:,k), dqvc_mix(:,k), &
720                        Tcontr(:,k), qcontr(:,k), qcontr2(:,k), fcontrN(:,k), fcontrP(:,k), &
721                        flight_dist(:,k), flight_h2o(:,k), &
722                        dcf_avi(:,k), dqi_avi(:,k), dqvc_avi(:,k))
723
724
725                  ELSE
726                  !--generalised lognormal condensation scheme (Bony and Emanuel 2001)
727
728                   CALL condensation_lognormal( &
729                       klon, Tbef, zq, zqs, gammasat, ratqs(:,k), &
730                       keepgoing, rneb(:,k), zqn, qvc)
731
732
733                  ENDIF ! .NOT. ok_ice_supersat
734
735                  DO i=1,klon
736                      IF (keepgoing(i)) THEN
737
738                        ! If vertical heterogeneity, change fraction by volume as well
739                        IF (iflag_cloudth_vert.GE.3) THEN
740                            ctot_vol(i,k)=rneb(i,k)
741                            rneblsvol(i,k)=ctot_vol(i,k)
742                        ENDIF
743
744
745                       ! P2.2.2> Approximative calculation of temperature variation DT
746                       !           due to condensation.
747                       ! Calculated variables:
748                       ! dT : temperature change due to condensation
749                       !---------------------------------------------------------------
750
751               
752                        IF (zfice(i).LT.1) THEN
753                            cste=RLVTT
754                        ELSE
755                            cste=RLSTT
756                        ENDIF
757                       
758                        ! LEA_R : check formule
759                        IF ( ok_unadjusted_clouds ) THEN
760                          !--AB We relax the saturation adjustment assumption
761                          !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
762                          IF ( rneb(i,k) .GT. eps ) THEN
763                            qlbef(i) = MAX(0., zqn(i) - qvc(i) / rneb(i,k))
764                          ELSE
765                            qlbef(i) = 0.
766                          ENDIF
767                        ELSE
768                          qlbef(i)=max(0.,zqn(i)-zqs(i))
769                        ENDIF
770
771                        num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT &
772                              +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)
773                        denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) &
774                              -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k)      &
775                              *qlbef(i)*dzfice(i)
776                        ! here we update a provisory temperature variable that only serves in the iteration
777                        ! process
778                        DT(i)=num/denom
779                        n_i(i)=n_i(i)+1
780
781                    ENDIF ! end keepgoing
782
783                ENDDO     ! end loop on i
784
785        ENDDO         ! iter=1,iflag_fisrtilp_qsat+1
786
787        ! P2.2> Final quantities calculation
788        ! Calculated variables:
789        !   rneb : cloud fraction
790        !   zcond: mean condensed water in the mesh
791        !   zqn  : mean water vapor in the mesh
792        !   zfice: ice fraction in clouds
793        !   zt   : temperature
794        !   rhcl : clear-sky relative humidity
795        ! ----------------------------------------------------------------
796
797
798        ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
799        IF (iflag_t_glace.GE.4) THEN
800           CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
801           distcltop(:,k)=zdistcltop(:)
802           temp_cltop(:,k)=ztemp_cltop(:)
803        ENDIF
804
805        ! Partition function depending on temperature
806        CALL icefrac_lscp(klon, zt, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
807
808        ! Partition function depending on tke for non shallow-convective clouds
809        IF (iflag_icefrac .GE. 1) THEN
810           CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), omega(:,k), qice_ini, ziflcld, zqn, &
811           rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
812        ENDIF
813
814        ! Water vapor update, Phase determination and subsequent latent heat exchange
815        DO i=1, klon
816            ! Overwrite phase partitioning in boundary layer mixed phase clouds when the
817            ! iflag_cloudth_vert=7 and specific param is activated
818            IF (mpc_bl_points(i,k) .GT. 0) THEN
819                zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
820                ! following line is very strange and probably wrong
821                rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
822                ! water vapor update and partition function if thermals
823                zq(i) = zq(i) - zcond(i)       
824                zfice(i)=zfice_th(i)
825            ELSE
826                ! Checks on rneb, rhcl and zqn
827                IF (rneb(i,k) .LE. 0.0) THEN
828                    zqn(i) = 0.0
829                    rneb(i,k) = 0.0
830                    zcond(i) = 0.0
831                    rhcl(i,k)=zq(i)/zqs(i)
832                ELSE IF (rneb(i,k) .GE. 1.0) THEN
833                    zqn(i) = zq(i)
834                    rneb(i,k) = 1.0
835                    IF ( ok_unadjusted_clouds ) THEN
836                      !--AB We relax the saturation adjustment assumption
837                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
838                      zcond(i) = MAX(0., zqn(i) - qvc(i))
839                    ELSE
840                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))
841                    ENDIF
842                    rhcl(i,k)=1.0
843                ELSE
844                    IF ( ok_unadjusted_clouds ) THEN
845                      !--AB We relax the saturation adjustment assumption
846                      !-- qvc (grid-mean vapor in cloud) is calculated by the condensation scheme
847                      zcond(i) = MAX(0., zqn(i) * rneb(i,k) - qvc(i))
848                    ELSE
849                      zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
850                    ENDIF
851                    ! following line is very strange and probably wrong:
852                    rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
853                    ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
854                    IF (iflag_icefrac .GE. 1) THEN
855                        IF (lognormale(i)) THEN 
856                           zcond(i)  = zqliq(i) + zqice(i)
857                           zfice(i)=zfice_turb(i)
858                           rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
859                        ENDIF
860                    ENDIF
861                ENDIF
862
863                ! water vapor update
864                zq(i) = zq(i) - zcond(i)
865
866            ENDIF
867
868            ! c_iso : routine that computes in-cloud supersaturation
869            ! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
870
871            ! temperature update due to phase change
872            zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
873            &     * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
874                  +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
875        ENDDO
876
877        ! If vertical heterogeneity, change volume fraction
878        IF (iflag_cloudth_vert .GE. 3) THEN
879          ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
880          rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
881        ENDIF
882
883
884    ! ----------------------------------------------------------------
885    ! P3> Precipitation processes, after cloud formation
886    !     - precipitation formation, melting/freezing
887    ! ----------------------------------------------------------------
888
889    ! Initiate the quantity of liquid and solid condensates
890    ! Note that in the following, zcond is the total amount of condensates
891    ! including newly formed precipitations (i.e., condensates formed by the
892    ! condensation process), while zoliq is the total amount of condensates in
893    ! the cloud (i.e., on which precipitation processes have applied)
894    DO i = 1, klon
895      zoliq(i)  = zcond(i)
896      zoliql(i) = zcond(i) * ( 1. - zfice(i) )
897      zoliqi(i) = zcond(i) * zfice(i)
898      ! c_iso : initialisation of zoliq* also for isotopes
899    ENDDO
900
901    !================================================================
902    ! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
903    IF (ok_poprecip) THEN
904
905      CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
906                            ctot_vol(:,k), ptconv(:,k), &
907                            zt, zq, zoliql, zoliqi, zfice, &
908                            rneb(:,k), znebprecipclr, znebprecipcld, &
909                            zrfl, zrflclr, zrflcld, &
910                            zifl, ziflclr, ziflcld, &
911                            qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
912                            dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
913                            dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
914                            dqsmelt(:,k), dqsfreez(:,k) &
915                            )
916      DO i = 1, klon
917          zoliq(i) = zoliql(i) + zoliqi(i)
918      ENDDO
919
920    !================================================================
921    ELSE
922
923      CALL histprecip_postcld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
924                            ctot_vol(:,k), ptconv(:,k), zdqsdT_raw, &
925                            zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, &
926                            rneb(:,k), znebprecipclr, znebprecipcld, &
927                            zneb, tot_zneb, zrho_up, zvelo_up, &
928                            zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, &
929                            zradocond, zradoice, dqrauto(:,k), dqsauto(:,k) &
930                            )
931
932    ENDIF ! ok_poprecip
933
934    ! End of precipitation processes after cloud formation
935    ! ----------------------------------------------------
936
937    !----------------------------------------------------------------------
938    ! P4> Calculation of cloud condensates amount seen by radiative scheme
939    !----------------------------------------------------------------------
940
941    DO i=1,klon
942
943      IF (ok_poprecip) THEN
944        IF (ok_radocond_snow) THEN
945          radocond(i,k) = zoliq(i)
946          zradoice(i) = zoliqi(i) + qsnowdiag(i,k)
947        ELSE
948          radocond(i,k) = zoliq(i)
949          zradoice(i) = zoliqi(i)
950        ENDIF
951      ELSE
952        radocond(i,k) = zradocond(i)
953      ENDIF
954
955      ! calculate the percentage of ice in "radocond" so cloud+precip seen
956      ! by the radiation scheme
957      IF (radocond(i,k) .GT. 0.) THEN
958        radicefrac(i,k)=MIN(MAX(zradoice(i)/radocond(i,k),0.),1.)
959      ENDIF
960    ENDDO
961
962    ! ----------------------------------------------------------------
963    ! P5> Wet scavenging
964    ! ----------------------------------------------------------------
965
966    !Scavenging through nucleation in the layer
967
968    DO i = 1,klon
969       
970        IF(zcond(i).GT.zoliq(i)+1.e-10) THEN
971            beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
972        ELSE
973            beta(i,k) = 0.
974        ENDIF
975
976        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
977
978        IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
979
980            IF (temp(i,k) .GE. t_glace_min) THEN
981                zalpha_tr = a_tr_sca(3)
982            ELSE
983                zalpha_tr = a_tr_sca(4)
984            ENDIF
985
986            zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
987            frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
988
989            ! Nucleation with a factor of -1 instead of -0.5
990            zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
991
992        ENDIF
993
994    ENDDO
995
996    ! Scavenging through impaction in the underlying layer
997
998    DO kk = k-1, 1, -1
999
1000        DO i = 1, klon
1001
1002            IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN
1003
1004                IF (temp(i,kk) .GE. t_glace_min) THEN
1005                    zalpha_tr = a_tr_sca(1)
1006                ELSE
1007                    zalpha_tr = a_tr_sca(2)
1008                ENDIF
1009
1010              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
1011              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
1012
1013            ENDIF
1014
1015        ENDDO
1016
1017    ENDDO
1018   
1019    !------------------------------------------------------------
1020    ! P6 > write diagnostics and outputs
1021    !------------------------------------------------------------
1022   
1023    !--AB Write diagnostics and tracers for ice supersaturation
1024    IF ( ok_ice_supersat ) THEN
1025      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs)
1026      CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:,k),RTT,2,.false.,qsati(:,k),zdqs)
1027
1028      DO i = 1, klon
1029
1030        IF ( zoliq(i) .LE. 0. ) THEN
1031          !--If everything was precipitated, the remaining empty cloud is dissipated
1032          !--and everything is transfered to the subsaturated clear sky region
1033          rneb(i,k) = 0.
1034          qvc(i) = 0.
1035        ENDIF
1036
1037        cf_seri(i,k) = rneb(i,k)
1038
1039        IF ( .NOT. ok_unadjusted_clouds ) THEN
1040          qvc(i) = zqs(i) * rneb(i,k)
1041        ENDIF
1042        IF ( zq(i) .GT. min_qParent ) THEN
1043          rvc_seri(i,k) = qvc(i) / zq(i)
1044        ELSE
1045          rvc_seri(i,k) = min_ratio
1046        ENDIF
1047        !--The MIN barrier is NEEDED because of:
1048        !-- 1) very rare pathological cases of the lsc scheme (rvc = 1. + 1e-16 sometimes)
1049        !-- 2) the thermal scheme does NOT guarantee that qvc <= qvap (or even qincld <= qtot)
1050        !--The MAX barrier is a safeguard that should not be activated
1051        rvc_seri(i,k) = MIN(MAX(rvc_seri(i,k), 0.), 1.)
1052
1053        !--Diagnostics
1054        gamma_cond(i,k) = gammasat(i)
1055        subfra(i,k) = 1. - cf_seri(i,k) - issrfra(i,k)
1056        qsub(i,k) = zq(i) - qvc(i) - qissr(i,k)
1057        qcld(i,k) = qvc(i) + zoliq(i)
1058      ENDDO
1059    ENDIF
1060
1061    ! Outputs:
1062    !-------------------------------
1063    ! Precipitation fluxes at layer interfaces
1064    ! + precipitation fractions +
1065    ! temperature and water species tendencies
1066    DO i = 1, klon
1067        psfl(i,k)=zifl(i)
1068        prfl(i,k)=zrfl(i)
1069        pfraclr(i,k)=znebprecipclr(i)
1070        pfracld(i,k)=znebprecipcld(i)
1071        d_q(i,k) = zq(i) - qt(i,k)
1072        d_t(i,k) = zt(i) - temp(i,k)
1073
1074        IF (ok_bug_phase_lscp) THEN
1075           d_ql(i,k) = (1-zfice(i))*zoliq(i)
1076           d_qi(i,k) = zfice(i)*zoliq(i)
1077        ELSE
1078           d_ql(i,k) = zoliql(i)
1079           d_qi(i,k) = zoliqi(i)   
1080        ENDIF
1081
1082    ENDDO
1083
1084
1085ENDDO ! loop on k from top to bottom
1086
1087
1088  ! Rain or snow at the surface (depending on the first layer temperature)
1089  DO i = 1, klon
1090      snow(i) = zifl(i)
1091      rain(i) = zrfl(i)
1092      ! c_iso final output
1093  ENDDO
1094
1095  IF (ncoreczq>0) THEN
1096      WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
1097  ENDIF
1098
1099
1100RETURN
1101
1102END SUBROUTINE lscp
1103!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1104
1105END MODULE lmdz_lscp
Note: See TracBrowser for help on using the repository browser.